Improving the Advection Discretization in D-FLOW Flexible Mesh

(1)

Improving the Advection Discretization in D-FLOW

Flexible Mesh

Jason Houtekamer

Supervisors:

Dr. Jeroen Gerrits, VORtech Dr. Mart Borsboom, Deltares Dr. Valeria Krzhizhanovskaya, UvA

A Thesis presented for the degree of Master of Science

Computational Science Lab University of Amsterdam

The Netherlands November 14, 2017

(2)

Declaration of Authorship

I, Jason Houtekamer, confirm that this thesis titled ”Improving the Advection Discretization in D-FLOW Flexible Mesh” is my own work, and I have documented clearly all sources used. I take full responsibility for the contents of this thesis, and the work done. I have acknowledged all main sources of help.

(3)

Abstract

The Netherlands has a rich history of water management. With advances in technology, this water management has become increasingly reliant on computer simulation. Daily predictions of water levels and currents are done by Rijkswaterstaat, using hydrodynamical simulation packages based on the shallow-water equations (SWE). D-FLOW Flexible Mesh (FM) is an engine for hydrodynamical simulations which has been in development at Deltares for the past decade, as a successor to the currently operational packages. D-FLOW FM was built to have the efficiency of the previous generation hydrodynamics models like Delft3D, combined with the modelling flexibility offered by a flexible mesh. The applications that D-FLOW was optimized for are large scale flows where water level differences and Coriolis force are the most important factors. Unfortunately it turned out that this optimization, combined with the highly unstructured meshes that are sometimes required to fit the domain, are causing errors in the advection terms, that are too large to be accepted. A new discretization of the advection term in the 2DH depth-integrated shallow-water equations is presented. The discretization is based on a numerical scheme that has full velocity vectors as unknowns on cell faces instead of only the normal velocity components, as is currently the case. This allows for a compact finite volume discretization. Two methods to reconstruct the velocity within the cells are investigated. A very simple method turned out to make large errors and did not converge. A velocity reconstruction method that is based on conservation of mass provides better results. Errors are reasonably small, and the convergence order is first order for the L∞ norm, and second order for the L1 norm. While the initial results are promising,

extensive further study is required before considering to implement this method into the D-FLOW FM code.

(4)

Acknowledgements

This thesis would not have been possible without the help of others. I would like to express my sincere gratitude to all the people who provided any help. In no particular order I would like to thank:

• my supervisor Jeroen Gerrits at VORtech, for all of the helpful advice, educational discussions and inspiring encouragement.

• my supervisor Mart Borsboom at Deltares, for being a patient teacher, providing honest criticism and sharing his excitement for numerical mathematics.

• my supervisor Valeria Krzhizhanovskaya, for allowing me to pursue this thesis. • everyone at VORtech, for providing advice and ideas.

• Sander van der Pijl and Frank Platzek at Deltares, for helping with MATLAB and D-FLOW related things.

(5)

1 Introduction

The Netherlands has a rich history of water management. With advances in technology, this water management has become increasingly reliant on computer simulation. Daily predictions of water levels and currents are done by Rijkswaterstaat, using hydrodynamical simulation packages based on the shallow-water equations (SWE). The shallow-water equations are a simplification of the Navier-Stokes equations that can be applied when the horizontal length scales are much larger than the vertical length scales. The term ‘shallow’ in the name is deceiving, since these equations can also be used for modeling ‘deep’ waters, such as oceans. D-FLOW Flexible Mesh (FM) is an engine for hydrodynamical simulations which has been in development at Deltares for the past decade, as a successor to packages such as Delft3D (one of the most commonly used packages) and SIMONA (the software used by Rijkswaterstaat). The operational software uses structured meshes, for which accurate and fast methods have been developed over the past decades. The downside of packages using structured meshes is that they are difficult to apply in areas with irregular geometry, such as estuaries. D-FLOW FM works on unstructured meshes, which can fit these irregular environments with far greater ease. The downside of finite difference and finite volume schemes on unstructured meshes is that they are not able to compute the advection terms very accurately if the mesh is irregular. Like many SWE packages, D-FLOW FM uses a staggered scheme to solve the equations. This means that the unknowns, water level and velocity, are computed at different locations on the mesh. The water levels are computed in the cell circumcenters, and the normal components of the flow velocities are computed on the faces of the cells. D-FLOW FM’s scheme, originally developed by Kramer and Stelling [21], based on previous work by Leendertse [22], Stelling [28], and Casulli [7]. It uses a finite volume discretization of the advection based on the scheme developed for the Navier-Stokes equation by Perot [24]. The downside of Perot’s scheme is that it only allows for a very inaccurate calculation of the advection term on irregular meshes, which causes D-FLOW FM to make relatively large errors where advection is dominant. Since the size of these errors is often larger than acceptable, alternative schemes should be investigated.

An alternative approach will be investigated, that is similar to D-FLOW, but has the full velocity vector as the unknown on the cell faces instead of only the normal component as is currently the case. Using the finite volume method, this seems to allow for a much more compact discretization with significantly less averaging and approximating of values. In applications where advection can be neglected, as is often the case in large open bodies of water, the properties of the scheme should be almost identical to the scheme currently used in D-FLOW. In the cases where convection plays an important role a notable improvement is required. This thesis will focus on developing details of the scheme and testing it using an implementation that focuses purely on the accuracy of the spatial discretization of the advection term employed in the scheme.

(8)

2 Problem Statement

Many numerical methods for solving the (2DH) shallow-water equations have already been studied in academic literature. Many more methods have been used to solve various other forms of the Navier-Stokes equations. These methods generally are not strictly better or worse than other methods but have their unique advantages and disadvantages.

D-FLOW FM was built to have the efficiency of the previous generation hydrodynamics models like Delft3D, combined with the modelling flexibility offered by the flexible meshes that had traditionally mostly been used in finite element methods. The applications that D-FLOW was optimized for are large scale flows where water level differences and Coriolis force are the most important factors. Unfortunately it turned out that this optimization, combined with the highly unstructured meshes that are sometimes required to fit the domain, are causing errors in the advection terms, that are too large to be accepted.

The objective of this thesis is to study a possible extension to the numerical scheme of D-FLOW FM, that would improve the accuracy of the advection discretization without compromising on the accuracy on any of the other parts. At the same time this should ideally be an extension in the sense that additions are made so that other components of D-FLOW can remain largely unchanged.

3 Background and Literature Study

3.1 History and Context

The focus of the work in this thesis will not be on applications, but rather on a small part of the mathematics underlying hydrodynamic simulation packages, namely the advection term. Nonetheless, it is useful to consider the context in which these packages are actually applied. D-FLOW FM and various other products made at Deltares are commercial products, and generally not developed with research in mind. As such the amount of academic literature created using the software is rather limited. Nonetheless there are some conference proceed-ings available. Jagers et al. [16] used D-FLOW FM to create a global 2D hydrodynamic model. Typically models are created for smaller regions, with open boundaries that allow for external forcing terms to simulate the influence of global tidal effects. The versatility of the D-FLOW’s flexible mesh allowed for a much larger global model while retaining more detail in coastal areas. The output of water level and flow data from hydrodynamics packages is also used as input for other models. Achete et al. [1] used D-FLOW FM simulations as input for a sediment transport model. A case study was performed on the Sacramento-San Joaquin Delta in the USA, where the model was 90% accurate in predicting the sediment transport behavior. The difference in scale between both of these applications illustrates very well why having a flexible mesh is very attractive.

(9)

Figure 1: The layout of the mesh used by Harlow and Welch, where h is the pressure and u and v are the velocity components normal to the cell faces in x and y direction respectively.

parts of the country are situated several meters below sea-level, and their safety not only depends on levees, but also depends on various barriers that close during storm surges. Using these barriers effectively requires accurate predictions of water levels at sea, but also in the various rivers that run through the country. For this reason, the Netherlands has a long-standing tradition of hydrodynamic modelling. Early in the 20th century, this modeling was done using physical scale models, but as computers became more powerful and more readily available, software modelling became more important than experimentation. These models generally used finite difference methods on Cartesian and curvilinear structured meshes, and the foundations were initially laid by Leendertse [22] in his PhD dissertation. Leendertse describes what would be, by current standards, a very simple finite difference discretization of the depth-integrated 2DH shallow-water equations. He used a staggered mesh, which had been pioneered only a few years earlier by Harlow and Welch [13] for the Navier-Stokes equations.

The method developed by Harlow and Welch worked on a mesh that computed normal flow velocity data on cell edges, and pressures in cell centers, thus staggering the variables as seen in figure 1. The motivation behind this idea was simple: The discretization of the pressure gradient term between 2 pressure data point is most accurate exactly in the middle of these 2 pressure data points. Since this gradient is part of the momentum equation it makes sense to place the velocity data exactly in this point. Similarly, the velocity derivatives in the continuity equation can be computed most accurately if the point at which they are approximated is exactly in the middle of the points used to approximate it. When finite difference methods are applied, staggered meshes also solve the odd-even decoupling that occurs on regular meshes.

Development of the finite difference models used in the Netherlands was continued by Stelling [28]. Stelling introduced upwind schemes into the methods by Leendertse, and refined the way

(10)

flooding and drying of cells was handled. His work provided the basis of the systems that are currently operational in the Netherlands, namely WAQUA [32] and Delft3D [12]. Both these packages use extensions of Stelling’s method for curvilinear and/or spherical meshes. The downside of these methods, and the meshes they require, is that their current implementation does deal well with very irregular boundaries.

Of course it is possible to add local refinement to these staggered meshes, but it makes the discretizations much more complex. Recently van der Plas [30] completed a PhD thesis in which he explored in great detail a method to do local refinement and to deal with boundaries cutting through cells. While his work focused on expanding the free-surface flow Navier-Stokes model in ComFLOW, a CFD package based on finite-volume methods, it seems very likely the method could be adapted to the 2DH shallow-water equations. He concludes that his methods should in theory contain first- and zeroth-order errors around the areas in which the refinement occurs. His numerical experiments, however, show the method converging must faster than expected based on analysis of local truncation error. This is known as supraconvergent behavior. Compared to D-FLOW FM, where these errors are present both in theory and in practice, this is a good sign. If finding a satisfactory solution to the issues in D-FLOW FM turns out to be impossible, applying the methods by van der Plas to the shallow-water equations could be worth exploring.

Many researchers, however, went in another direction and focused on developing methods using unstructured meshes. Casulli and others developed UnTRIM [10], a package that uses a hybrid of finite-difference and finite-volume methods to solve the Navier-Stokes equations on unstructured staggered meshes. This is one of the methods that D-FLOW FM is based on. Another example is FVCOM [9], a finite-volume based shallow-water model that also contains temperature and salinity as part of the core computational method. FVCOM is somewhat similar to the method detailed in this thesis, in that it uses the complete velocity component. The control volumes and layout of the variables is very different from D-FLOW FM however.

Luettich, Westerink and Scheffner developed ADCIRC [23], a finite-element based shallow-water circulation model. ADCIRC was developed for the US government, with the goal of long (a year) numerical simulations on large domains (the east coast of the US). ADCIRC is currently still being used, in a computational model that integrates it with SWAN [11], to predict hurricane waves and storm surges, and the results can be interacted with online. Some other examples of flexible/unstructured mesh packages for the shallow-water equations are TELEMAC [14], ELCIRC [3] and SLIM [6]. These models all have their own unique and interesting approaches to solving the shallow-water equations, but they are all quite different from the approach taken in D-FLOW FM. While this literature is interesting, it does not provide much help when it comes to improving the advection discretization in D-FLOW FM because the improvements to the advection should be made while maintaining as much of the current scheme as possible.

(11)

3.2 Numerical Methods in D-FLOW FM

D-FLOW FM was presented in 2011 by Kernkamp, van Dam, Stelling and de Goede [17], in a paper that roughly outlines the numerical methods used but mostly focuses on its applications and performance. The main reason that the numerical method is not discussed in great detail, is that it is largely based on previous work, as will be outlined in the rest of this section. According to the paper the main novelties are the combining of 1D, 2D and 3D modelling concepts, allowing arbitrary degree polygon shaped cells and using a very efficient matrix solver tailored to the method. The mathematical details of the numerical scheme used in D-FLOW FM will be explored in section 4. This section of the literature review will discuss the various papers that provided the basis for the numerical methods in D-FLOW, and their contributions.

The mesh used in D-FLOW FM is an unstructured orthogonal Arakawa-C mesh [2]. This type of mesh, which has normal velocity data on cell faces and water level data in cell centers, was applied to the Navier-Stokes equations by Casulli and Zanolli [8]. They use a simplified version of the Navier-Stokes equations with assumptions very similar to those made in the derivation of the shallow-water equations. The reason that orthogonal meshes are appealing, is that the pressure gradient term in the momentum equation only requires information from the 2 cell centers adjacent to the face for which the equation is solved. Another novelty was their application of a very efficient semi-implicit time integration scheme to the shallow-water equations. Some of the variables are computed implicitly while others are computed explicitly. The variables that are computed implicitly are chosen such that conditional stability is guaranteed. The resulting matrix is positive definite which allows for very fast and efficient solution methods.

One of the main weaknesses of staggered schemes, however, is that while they are very efficient for smooth flows, they do not perform very well for rapidly varying or critical flows on unstructured meshes. When phenomena like hydraulic jumps occur in these rapidly varying flows, the hydrostatic pressure assumption inherent to the shallow-water equations is locally invalid. Stelling and Duinmeijer [29] created some extensions to the classical staggered mesh implicit schemes, that can deal with these flows by locally applying a scheme that is based purely on conservation laws.

Another issue that occurs when using staggered meshes, is that the full velocity vector is known nowhere on the mesh, while it is required to compute the Coriolis term. The Coriolis term is one of the terms driving the tidal forces which shallow-water models are often used to simulate. The challenge of reconstructing the velocity vector on an unstructured mesh is a topic on its own in the literature. The full velocity vector at cell centers (or circumcenters) is reconstructed from the normal velocities at the surrounding faces. Then if, required, the reconstructed full velocity vector at the cell center is projected back to the cell faces tangential direction. This results in two different projections at every cell face, and the final tangential velocity will be some linear combination of these two projections. Klepstova, Pietrzak and Stelling [19] compare several methods, and conclude that the method by Perot [24] is the most suitable. Perot’s reconstruction, and accompanying proofs of various conservation properties,

(12)

however, were done for an application to the Navier-Stokes equations.

The full velocity vector is also needed for computing the advection term. Kramer and Stelling [21] investigated the properties of an advection scheme based on the Perot reconstruction. They claim that such a scheme can be made to conserve momentum if the water level data is chosen properly. This is the scheme that the approach to advection in D-FLOW FM is based on.

In an internal Deltares memo by Borsboom [4], a very detailed analysis of the numerical scheme used in D-FLOW FM showed that momentum is not conserved. Borsboom also notes that conservation of momentum would be possible with minor changes, and implies that not being perfectly conservative is a conscious choice. This on itself is not a huge issue; a discretization can be accurate while not fully conserving momentum and the sacrifices required to achieve momentum conservation might outweigh the benefits. Unfortunately it is hard to judge exactly what these sacrifices and benefits are, because they depend on the typical application, mesh, computing time and other such factors. Unfortunately, asides from conservation properties, the Perot method is found to be only first-order accurate in general. Borsboom and van der Pijl [5] found, in another internal Deltares memo, that in certain test cases the current advection scheme in D-FLOW FM is very inaccurate. The test case they considered is a Poiseuille flow, a steady laminar viscous flow through a straight, constant depth, channel. While the results are acceptable when using a uniform rectangular mesh, the results when using a mixed rectangular-unstructured mesh are very poor. These reports motivated an effort to improve the accuracy of the advection discretization. Visee [31] explored several techniques to improve on the Perot velocity vector reconstruction per cell. These techniques, however, mostly increased the number of mesh cells involved in the discretization of the equations on each cell face, while providing only a small increase in accuracy. This brings us to the goal of this thesis, which is to explore alternative approaches to improve the discretization of the advection in D-FLOW FM.

3.3 Improving the Advection Discretization

Most of the problems with the advection discretization in D-FLOW FM are caused by the fact that there is insufficient velocity information to accurately compute all the required terms accurately on a highly irregular mesh. A potential solution could be to compute the full velocity vector at the locations where previously only the normal velocity was computed. This would provide much more velocity information for the computation of the advection terms. A method that uses the full velocity vectors at the faces was proposed by Hwang [15] for the incompressible Navier-Stokes equations. Hwang uses finite volumes to discretize both the momentum and continuity equations on different volumes. Figure 2 shows the control volumes used for the continuity and momentum volumes respectively. Note that figure 2 shows circumcenters, while Hwang uses centroids.

Ride et al. [25] tested 2 different staggered schemes with full velocity vectors. One similar to the scheme used by Hwang with the full velocity vector stored in the centers of the faces as

(13)

Figure 2: On the left the control volume for the continuity equation (the mesh cell), on the right the control volume for the momentum equation (momentum volume).

seen in figure 2, and one with the full velocity vectors stored at the vertices of the continuity volume. Rida et al. also use an advecting velocity field u that is reconstructed in a way that ensures conservation of mass.

Shala et al. [26] similarly use a finite volume based scheme that stores the full velocity vector on the faces of the continuity volumes. They, however, opt to use the two adjacent continuity volumes as the momentum volume around every face. The advantage of this method is that no interpolation or reconstruction of velocities is required at all. In the methods used by Hwang and Rida, the velocity information is needed within the continuity volumes at the locations xk as seen in figure 2, while it is not known there. Thus some interpolation is

required. The method used by Shala on the other hand requires velocity information at the midpoints of the cell faces, which happens to be exactly where the velocity information is stored. So the full velocity vectors are known at every point where they are needed. The downside of this approach is that momentum conservation cannot be guaranteed because the control volumes overlap.

This thesis contributes by applying a scheme similar to those applied by Hwang [15] and Rida [25] to the 2DH depth-integrated shallow water equations. Such a scheme could improve a weakness of D-FLOW FM, the accuracy of the advection discretization on highly irregular meshes, while maintaining its strengths.

(14)

4 Numerical Methods

4.1 Discretization in D-FLOW FM

The continuity equation is discretized by a finite volume method, with the control volume being a single cell in the mesh as seen in the left part of figure 2. The index i is used to indicate variables located at circumcenters, j indicates variables located at the midpoints of cell faces and k indicates the center of a momentum face. The total water depth hi is

computed at the circumcenter xi, and can be split as h = ζ + b. Since b, the bathymetry, is

constant in this study, the variable that is actually solved for is ζi. The discretization of the

continuity equation at cell i is given by: dζ

dtVi = − X

j

ljhjuj, (1)

where j are the indices of the cell faces that make up cells i. lj is the length of the face, hj

and uj are the total water depth and normal velocity respectively, and Vi is the volume of

cell i. Note that hj has to be taken from either cell i, or from the upwind cell, depending on

the direction of flow, because the water depth h is not stored at the faces j.

Because the normal velocity uj is one of the unknowns, the continuity equation is easily and

compactly discretized. In order to solve the momentum equation for uj, a discretization of

the water-level gradient ∇ζ in the direction normal to the face is required. This becomes very simple when ζi is stored in the circumcenter, because the line connecting the circumcenters

of two adjacent triangles is normal to their shared face. This requires the circumcenter to exist, which it always does for triangles, but there are some restrictions on higher order polygons. For triangles, it is still preferable for the circumcenter to be within the triangle, and as close to the centroid as possible, since this minimizes the discretization error. After all, in the continuity equation it is assumed that the water level hi is representative for the

water level in all of cell i. Orthogonal meshes fullfill these requirements, although they do not guarantee that the circumcenter is close to the centroid, it is always within the polygon. The discretization of the momentum equation is now also relatively compact. The discretization used in D-FLOW is duj dt + Advec(u) + g ζi,R− ζi,L ||xi,R− xi,L|| + g||˜u||uj hjC2 = 0, (2)

where Advec(un) is an explicit discretization of the advection roughly following the Perot method [24] and ˜u is a reconstruction of the full velocity vector at the center of face j. The subscripts R and L indicate the cells to the right and left of face j.

When discretizing this equation in time, a compact, symmetric and diagonally-dominant system of equations is obtained that can be easily solved. Many of the choices in the numerical scheme of D-FLOW FM seem to have been made with the goal of achieveing such a system

(15)

Figure 3: The control volume for the advection term.

of equations. Unfortunately this seems to be detrimental to discretization of the advection term.

4.2 A Different Approach to the Advection Term

The approach to the advection that we are exploring in this thesis discretizes the advection term using the finite volume method (FVM) on the marked volumes that are shown in Figure 3. The two triangles are continuity volumes. Figure 4 shows the locations of, and notation used to refer to, several important points on the mesh. The tangential and normal velocity data is stored at the midpoints xk of the faces fj of the continuity volumes. Normal

vectors are defined outward for the normal component, and in counterclockwise direction for tangential components. The reason to compute values at the midpoints of faces is because that it is convenient for the purposes of a one-point Gauss approximations used to compute the fluxes of mass and momentum across a cell face. The marked quadrilateral in figure 3 is a momentum volume. The momentum volume is created by connecting the two mesh nodes at the ends of the face to the circumcenters xi of the two continuity volumes adjacent to the

face. This method is similar to those used by Rida et al. [25] and Hwang [15]. One of the main reasons to use FVM, is that it guarantees conservation of momentum, both globally and locally, because any momentum leaving a momentum cell through a face fk is added to

the adjacent momentum cell that shares face fk. While conservation can be an important

property, it does not automatically imply accuracy. In order to determine the accuracy and convergence behaviour of the method, mathematical analysis and numerical testing is required.

While the placement of variables on the mesh for this method is different from the current D-FLOW mesh, in that it has not only the normal but also the tangential velocity components at the faces, this change should bring significant advantages with it. Not only the advection discretization can be improved, but diffusion and Coriolis also require tangential velocity information and can potentially be improved.

(16)

Figure 4: Two continuity cells and a momentum cell.

4.2.1 Discretizing the Advection Term

Recall that the shallow-water momentum equation is given by ∂q

∂t + ∇ · (qu) + gh∇ζ = τττ , (3)

where q = hu and τ is a broad term representing various stresses like viscosity and external forces like Coriolis and bottom friction. Integration of the advection term over a momentum volume yields:

Z

V

∇ · qudV. (4)

The main reason to write q instead of hu here is to make clear that these represent different things physically. q is the transported quantity while u is the transporting quantity. In an actual implementation it might be beneficial to use an upwind approximation for q in order to prevent wiggles and increase numerical stability.

Part of the FVM approach is to transform a volume integral to a contour integral using the divergence theorem. Applying this theorem to the advection integral yields

Z

S

(u · n)qdA, (5)

where A is the border of the volume and n is the normal to S. Because the volume V is a quadrilateral consisting of four straight line segments, this integral can be written as a sum over those four line segments. The individual integrals over the faces of the momentum volume can be numerically approximated using a one-point Gaussian quadrature rule. The reason to use a one-point quadrature is because it is computationally very efficient and still second-order accurate. The velocities that are used in the quadrature will also be at most second-order accurate, so the one-point approximation should not be the bottleneck for the accuracy. Using the one-point Gauss quadrature yields:

(17)

Figure 5: Two continuity cells and a momentum cell, with faces and important points named.

Z

k

(u · n)qdfk= (qk· nk)uklk, (6)

where qm = hum and lm is the length of the face. Note that the face water level hk is

determined by the continuity cell in which the momentum face is located. The entire integral can now be computed as the sum over the four momentum faces surrounding the cell face for which we are solving the momentum equation:

Z V ∇ · qudV = 4 X k=1 (qk· nk)uklk. (7)

lk and nk are both properties of the mesh geometry and h is assumed to be 1 for simplicity.

These quantities are all known. uk however, are not known, but uj are known instead. How

to compute uk given these uj will be discussed in the section 4.2.2.

4.2.2 Reconstructing the Velocity Field

In order to find velocities at the faces of the momentum volume, a very simple approximationis an average of the velocities at the two adjacent cell faces. For example, in figure 5 the velocity on the N-E face could be computed as the average of the velocities on the NE and C faces. There are some subtle concerns. Namely, do the fluxes computed using these velocities adhere to conservation of mass as it is enforced by the continuity equation?

Equation 7 shows that qk and uk are required in the centers of the momentum faces. Since

the velocity unknowns are not available there, they have to be constructed from the velocities at the cellfaces in some fashion. Because the continuity equation 46 has to hold everywhere, the discrete version eq. 1 should also hold for the momentum volumes. With the assumption h = 1 that means that

(18)

X

k

lkqk· nk= 0 (8)

should hold. This provides a good constraint on the computation of the advecting uk. The

advected qk can then be taken in some upwind fashion, but that is outside the scope of this

study.

The problem of reconstructing velocities within cells in a such a conservative manner is also encountered when trying to track particles using streamlines. In 2016 Ketefian, Gross and Stelling [18] proposed a method to reconstruct the velocity field within a cell such that continuity is not just guaranteed per cell, but also across cell faces. The method they propose is second-order accurate and linearly exact. This is done by not only using the fluxes through the cell faces, but also using the linear variation of the flux across the face. The continuity volumes have three faces, so the flux and linear variation of the flux provide enough information to compute the 6 unknowns of a linear velocity field. A note that should be made here is that this conservative velocity reconstruction does some averaging of velocity information from neighboring cells that could be detrimental to the accuracy.

For any cell i a linear velocity field ui = (ui, vi) has the form

ui(x) = Aix + bi, (9)

where x = (x, y) and Ai and bi are given by

Ai = a_i bi ci di , bi = b_x,i by,i . (10)

In order to compute such a velocity field for, for example, the cell with circumcenter E in figure 5, the 6 required quantities would be the momentum fluxes through the edges QNE, QSE

and QC, and the linear variation of the momentum along every edge dqNE/dsNE, dqSE/dsSE

and dqC/dsC. The capitalization is used to indicate that the quantity is integrated over a

face. The 6 × 6 matrix that relates these quantities to the velocity field contains geometric information about the cell. Finding the velocity fields requires this 6×6 matrix to be inverted for every cell. Since the matrix only depends on geometry it does not change after the mesh has been created, so this step only has to be completed once upon mesh generation.

In order to apply this velocity reconstruction, the linear variation of the momentum in the di-rection tangential to every continuity cell face has to be determined. Because the full velocity vectors are known at the continuity volume faces this can be done relatively straightforwardly. Once again referring to figure 5, dqC/dsC can be approximated using finite difference from

the velocity information at NE, SE, NW and SW. Assuming that the mesh is reasonably regular, eg. the circumcenters are close to the centroids, the middle between NE and NW will be roughly on face C at 3/4th of its length, and the middle between SE and SW will be roughly on face C at 1/4th of its length. So we can write

(19)

dqc ds = qN W + qN E 2 − qSW + qSE 2 /lC 2 = qN W + qN E− qSW − qSE lC . (11)

In this expression q is the full momentum vector. Since only the normal momentum compo-nent is required, a dot product with nC finally gives dqC/dsC as it is needed to reconstruct

the velocity field.

This method is exact for linear velocity fields, but how it performs for more complex velocity fields will be examined in the results section.

4.2.3 Analytical Solution for the Advection Integral

To find an analytical solution to the advection integral, it is convenient to start with the form created by applying the divergence theorem

Z

S

q · nudS, (12)

where instead of using a one-point Gaussian quadrature rule, a line integral can be solved for every face. For convenience, assume h = 1, so that the integral can be written as

Z

S

u(x, y)(u(x, y)nx+ v(x, y)ny)dS, (13)

where nx and ny are constant because the faces are straight lines. Using the fact that faces

are straight makes parameterizing them straightforward. Let c1 and c2 be the begin and end

points of any face of the momentum volume, in a counterclockwise manner, as seen in Fig 3. The parameterization of these faces is then given by

c(s) = c1(s) + s(c2− c1) (14)

where c1 = (x1, y1) and c2 = (x2, y2) are the coordinates of the points c1 and c2. Given a

velocity field u(x, y) we can substitute the parameterization into the integral Z

S

u(c(s))(u(c(s))nx+ v(c(s))ny)dS. (15)

Using dS = |c2− c1|dt, the integral over s becomes

Z 1

0

u(c(s))(u(c(s))nx+ v(c(s))ny)|c2 − c1|dt. (16)

Given a specific velocity field u(x, y), this integral can be easily evaluated using some com-puter algebra system to give an expression for in terms of c1 = (x1, y1) and c2 = (x2, y2).

(20)

corresponding integrals for each momentum volume will provide the value of the advection integrated over that volume.

4.2.4 Error Analysis and Convergence Behaviour

It is useful to do some analysis of the expected errors and convergence behavior. This will give some insights in why the errors are the way they are.

Figure 6 shows a momentum volume with the size of some relevant local mesh-related dis-tances indicated. Given a velocity field u(x, y), the series expansion of its x component u around the center of this momentum volume, which we assume to be (0, 0) is:

˜ u(x, y) =u0+ xux+ yuy+ 1 2x 2_u xx+ xyuxy + 1 2y 2_u yy+ 1 2xy 2 uxyy + 1 2x 2 y2uxxyy+ 1 4x 2 yuxxy+ O(x, y)3, (17)

where u0 and all the derivatives are evaluated at the center of the momentum volume. Given

this approximation, we can examine the numerical integration error in the advection integral made by using the one-point Gauss quadrature. Consider the point value of a single term of the advection in the center of the momentum volume. The exact value of the integrated ∂uu/∂x term of the advection using the series expansion is

1 V Z V ∂ ˜u2_{(x, y)} ∂x dV. (18)

The numerical approximation is a sum of 4 one-point quadratures over the faces, which we can evaluate, so the error can then be written as

1 V Z V ∂ ˜u2(x, y) ∂x dV − 1 V ˜ u ∆x 2 , (∆y + δy) 2 (∆y + δy) ˜ u ∆x 2 , (∆y + δy) 2 + ˜ u ∆x 2 , − ∆y 2 ∆y ˜ u ∆x 2 , − ∆y 2 + ˜ u −∆x 2 , (∆y + δy) 2 (−∆y − δy) ˜ u −∆x 2 , (∆y + δy) 2 + ˜ u −∆x 2 , − ∆y 2 (−∆y) ˜ u −∆x 2 , − ∆y 2 . (19) Evaluating both of these expressions, and subtracting them, will give us an idea of the kind of terms that will occur in our error. Mathematica was used to do this, and the result it

(21)

Figure 6: 2 Continuity volumes and the momentum volume around their shared edge. produces has many terms. The interesting terms are mainly the lowest order terms, because they are leading, so we will focus on those.

The first interesting case to consider is δy = 0, this occurs when the mesh is regular, but it also occurs when an arbitrary triangular mesh is refined. Figure 7 shows the same set of 2 mesh cells used in figure 6, and their shared momentum volume. When the mesh is refined, every mesh cell is refined into 4 smaller cells, these 4 cells will all have the same shape as the original cell. The same does not apply to the momentum volumes however. The new momentum volumes that are on the original edge between the 2 continuity volumes retain the same shape as the original momentum volume, but are 4 times smaller. In these cases not only do ∆x and ∆y halve their sizes, δy also goes down in size by a factor 2. This will be important later on when looking at the convergence of error terms that contain δy. In the middle of the old unrefined mesh cell there is now an edge connecting 2 new mesh cells that have the same shape. The momentum volume around this edge will be perfectly regular, so δy = 0. As the mesh is refined further, a larger and larger part of the mesh will consists of these regular volumes. The contribution of these irregular volumes to the total area of the domain goes down by a factor 2 at every refinement step because the number of irregular volumes doubles, but their area goes down by a factor of 4.

The error terms that remain in the δy = 0 case are: 1 12∆y 2_u xuyy+ 1 6∆y 2_u yuxy + 1 12∆y 2_uu xyy − 1 4∆x 2_u xuyy+ O(∆x, ∆y)4 (20)

These are all second-order terms or higher, as expected. Interesting to note is that even among the higher order error terms, there are no odd-powered error terms (x3_{, y}5 _{etc.). The}

(22)

Figure 7: On the left 2 continuity volumes and the momentum volume around their shared edge. On the right these same 2 continuity volumes have been refined into 8 volumes. Momentum volumes are drawn in to indicate the increasing regularity in their shapes.

1 3δyuxuy+ 1 3δyuuxy+ 1 6 δy2 ∆yuxuy+ 1 6 δy2 ∆yuuxy + 1 8δy 2_u xuyy+ 1 4δy 2_u xuxy + 1 8δy 2_uu xyy+ 1 24 δy3 ∆yuxuy + 1 12 δy3 ∆yuyuxy + 1 24 δy3 ∆yuxyy+ 1 8δy∆yuxuyy+ 1 4δy∆yuxuxy + 1 8δy∆yuuxyy+ 1 12∆y 2 uxuyy+ 1 6∆y 2 uyuxy + 1 12∆y 2 uuxyy− 1 4∆x 2 uxuyy+ O(∆x, ∆y)3 (21) Notice that the final row of terms in Eq. 21 are actually the exact same error terms we saw in the δy = 0 case (20). All the other terms however are additional error terms caused by the irregularity in the shape of the momentum volume. Most notably there are now also 2 first-order terms in δy, 2 first order δy2_{/∆y terms, and second-order terms. These terms}

occur here because the irregularity of the cells means they no longer cancel out. What we can conclude from this is that the L∞ norm should converge at a first-order rate, but the L1

norm should show second-order convergence behaviour. For definitions of these norms see section 4.2.5. One of the underlying assumptions here is that the reconstructed velocity used for the numerical approximation is also at least second-order accurate, otherwise the error in the velocity could be dominant. An nth-order error in velocity leads to a local n − 1th-order error in advection. So a first-1th-order accurate velocity reconstruction, can lead to local zeroth-order advection errors on an unstructured mesh.

In the case of the very simple velocity reconstruction that was tested, the velocity on the momentum face is the average of the velocities at the 2 adjacent cell faces, there is some

(23)

added complication. Figure 6 shows that the points where the velocities are approximated are not the same points that are used in the one-point Gauss integration. Instead these points are separated by some distance that scales with ∆x and ∆y, and a constant ξ that depends on the geometry of the continuity volumes (and the momentum volumes, but their geometry in turn is determined by the continuity volumes). For simplicity, and to keep the number of terms somewhat limited, it is assumed that there is only some error in the x direction. The error made in the integration when using this velocity reconstruction can then be written as

1 V Z V ∂ ˜u2_{(x, y)} ∂x dV − 1 V ˜ u ∆x 2 + ξ1∆x, ∆y 2 ∆y ˜ u ∆x 2 + ξ1∆x, ∆y 2 + ˜ u ∆x 2 + ξ2∆x, − ∆y 2 ∆y ˜ u ∆x 2 + ξ2∆x, − ∆y 2 + ˜ u −∆x 2 + ξ3∆x, ∆y 2 (−∆y) ˜ u −∆x 2 + ξ3∆x, ∆y 2 + ˜ u −∆x 2 + ξ4∆x, − ∆y 2 (−∆y) ˜ u −∆x 2 + ξ4∆x, − ∆y 2 . (22)

There are 4 distinct ξi constants to indicate that they are not necessarily all the same. They

would be the same (non-zero) size if the triangles were equilateral, but if the triangles are skewed at all this symmetry is lost. Because the shape of the mesh cells does not change when refining the mesh, any error terms that depends on ξ and not on ∆x or ∆y are zeroth-order constant errors. This is very similar to what happens in the current D-FLOW FM discretization. A first-order error in the velocity used in the advection discretization leads to a local zeroth-order error in the advection. The resulting error terms are shown in appendix D, because the list is too long to show here.

Some of the terms are constant, as they do not contain ∆x and ∆y. Useful to note is that all these constant zeroth-order terms cancel out when ξ1 = ξ2 = ξ3 = ξ4 = ξ, which

is the case when dealing with equilateral triangles. So in the case of a mesh consisting of purely equilateral triangles this reconstruction will be first order accurate. In the case of the unstructured meshes this method will have local zeroth-order errors. The remaining terms in when all ξ are the same are

1 12∆y 2_u yyux+ 1 6∆y 2_u yuxy + 1 12∆y 2_uu xyy − 1 4∆x 2_u xuxx

+ 2∆xξu2_x+ 2∆xξuuxx+ 3∆x2ξ2uxuxx+ O(∆x, ∆y)3

(23)

While this analysis provides clear grounds not to use the very simple velocity reconstruction of averaging velocities from adjacent cells, a very similar analysis would also apply when

(24)

considering an upwind scheme. Upwind schemes are often used to prevent wiggles and add dissipation to the method. In upwind schemes the transporting quantity is computed by methods such as the ones outlined above, but the transported quantity is taken often the ”upwind” continuity edge (or in case of higher order methods it is taken from several continuity edges as some weighted average). When using an upwind scheme with the methods outlined above, zeroth-order constant errors would be present unless the upwind velocities are chosen from points such that these errors cancel. Detailed numerical and analytic analysis of upwind schemes has not been performed because it is out of the scope of this thesis. 4.2.5 Error Norms

Useful measures for comparing errors are the L1- and L∞- norms. The L1 norm is useful

because it shows the average behaviour of the method, and this norm should thus not be dominated by the effects of very irregular cells on the accuracy of the method. The L∞ norm

is useful because it shows exactly this effect of irregularity on the accuracy of the method. The norms are computed as

L1 = 1 V Z ∂ ˜u2 ∂x − ∂u2 ∂x dV =⇒ 1 n X i   R i ∂ ˜u2 ∂xdVi− R i ∂u2 ∂xdVi Vi  , (24) L∞= max i   R i ∂ ˜u2 ∂xdVi− R i ∂u2 ∂xdVi Vi  , (25)

where V represents the area of entire domain, ∂ ˜u2/∂x and ∂u2/∂x respectively are the numerical approximation, and true value,of the advection term and n is the total number of cells in the mesh.

(25)

5 Implementation

All the implementation is done in MATLAB. The main goal of the thesis is to explore the accuracy of advection schemes, not necessarily to create a highly optimized large-scale code. Therefore performance is not as important as the ability to quickly prototype. Nonetheless, it is preferable for the code to run relatively quickly so experiments can be repeated easily. Because MATLAB is highly optimized for doing linear algebra, vector and matrix multipli-cations are used instead of for-loops whenever possible. For symbolic analysis like in section 4.2.4 Mathematica was used. Mesh generation was done by Sander van der Pijl at Deltares using the D-FLOW FM Mesh generation tool.

Some of the used MATLAB code had previously been created by Sander van der Pijl for various experiments related to D-FLOW FM, most notably functionality to read the NetCDF mesh files already existed. The way the data is structured is straightforward: All variables are stored as vectors. The vectors with velocities, normal vectors, x, y-coordinates etc. contain the information about cell face j at position j. Similarly there are vectors momentum faces, where the information about moment face k is always at position k. Finally, information about cell i, and it’s circumcenter, will be at position i in the respective vectors. Most of the geometry information is stored in matrices that relate the various vectors to each other. For example, there is a Ni× 3 matrix, where the 3 entries of row i contain the indices j of the

3 faces of cell i. In many other cases, sparse matrices with 1s and 0s are used to represent some operator. For example when the total advection on a moment face is computed, the contributions from each of the momentum faces are stored in a vector. There are then 2 k × j matrices which contain a 1 at location k, j if momentum face k contributes to the advection for cell face j, and a 0 otherwise. Construction of these matrices is generally nontrivial, and one has to be very careful to enforce consistent conventions about minus signs and directions. Once these matrices are constructed, writing the actual methods is generally just a matter of doing a series of matrix vector multiplications.

The analytical values for the advections were computed using MATLABs symbolic features. The velocity field was integrated over every face analytically. The resulting expression only depended on the x and y coordinates of the end points of the face. These symbolic integrals are by far the most time consuming part of the code, but luckily they only have to be ran once for every velocity field/mesh combination.

(26)

6 Results

The general procedure followed in all the performance tests is to prescribe some divergence free velocity field u(x, y), and then use that to compute the analytical advection value. Next the velocity field is used to set the velocities at the center of each cell face j using the coordinates xj. A small issue with this approach is that for velocity fields that are quadratic

or higher order, the discrete continuity equation no longer holds. From there the velocities at the momentum faces k are computed by first computing dqj/dsj at every cell face, then using

the velocity field reconstruction and finally inserting the coordinates xkinto the reconstructed

velocity field. The mesh used for the test is an unstructured triangular mesh with 0 ≤ x ≤ 50 and 0 ≤ y ≤ 100. Figure 8 shows the mesh. Norms are computed over the x- and y-terms of the advection, but plots are made for all the separate terms to be able to more easily isolate errors.

Whenever the error in the x-advection term or y-advection term is mentioned, this refers to the error made in the entire x or y advection term. So for the x-error those terms are ∂uu/∂x + ∂uv/∂x and for the y-error they are ∂uv/∂y + ∂vv/∂y. To compare errors made on various velocity fields relative errors can be useful. However in some of these velocity fields there will be locations with 0 advection. In order to have some alternative means of creating relative errors to compare, a sort of typical value is constructed that will be non-zero for all velocity fields with non-constant components. The typical value used here is:

TypicalVal = Z 100

0

Z 50

0

|u(x, y)| + |v(x, y)| 100 · 50 dxdy· Z 100

0

Z 50

0

|ux(x, y)| + |uy(x, y)| + |vx(x, y)| + |vy(x, y)|

100 · 50 dxdy.

(26)

The reason this particular value was chosen, is that it is somewhat similar to ∇·(qu) integrated over the entire domain in a way that guarantees a non-zero value for non-zero u. The relative error in cell i δAdvi can then be computed using the same errors that were used in the norms

in eqs 24 and 25 as follows

δAdvi = 1 TypicalVal   R i ∂ ˜u2 ∂xdVi− R i ∂u2 ∂xdVi Vi  . (27)

For the simplest test case, a constant velocity field, the method should be exact. As seen in equation 20 and 21 every term of the error made in the one-point Gauss integration contains some derivatives of u. These terms are all zero for a constant velocity field, so the error should also be zero. Numerical results confirm this, with all terms being accurate up to machine precision. When the mesh is refined, the error norms do get larger, but this is simply because more machine precision rounding errors are being added up. Even the simple

(27)

0 20 40 x 0 20 40 60 80 100 y

Figure 8: The mesh used in the various tests

velocity reconstruction method is accurate in this case, since the fact that the velocity is the same everywhere means it is actually estimated correctly.

6.1 Linear Velocity Field

The next step is a linear velocity field. The linear velocity field that was studied is u(x, y) = (0.01x + 2, −0.01y + 2), figure 9 shows a vector plot of this velocity field. This field is chosen because it is divergence free, so ∇ · u = 0. The velocity within the entire domain is between 1 and 3, which are values one could expect in a real application. Another linear velocity field was also studied, but because both velocity fields showed identical behavior analyzing the second field would not add any value.

First we examine the convergence behavior. Figure 10 shows the convergence behavior of the x and y terms of the advection for the linear velocity field u(x, y) = (0.01x + 2, −0.01y + 2). The convergence is exactly as predicted by the numerical analysis. The velocity reconstruc-tion is exact in this case, so the only error should be the integrareconstruc-tion error from the one-point Gauss approximation. The error norms remain roughly constant when using the simple ve-locity reconstruction, marked in the figures as ”x simple” and ”y simple”. When using the conservative velocity reconstruction, marked in the figure just as ”x” and ”y” we see first-order convergence of the L∞ norm and second-order convergence of the L1 norm. The plot

actually shows a convergence slightly faster than second order for the conservative velocity reconstruction. This is caused by the fact that every time the mesh is refined, the total volume of the domain that is considered in the norms changes slightly because of the way

(28)

0 20 40 x 0 20 40 60 80 100 y

Figure 9: A vector plot of the linear velocity field u(x, y) = (0.01x + 2, −0.01y + 2)

boundaries are handled.

1st refinement 2nd refinement 3rd refinement 0 1 1.5 2 Convergence order L 1 x L 1 x simple L ∞ x L ∞ x simple L 1 y L 1 y simple L ∞ y L ∞ y simple

Figure 10: The convergence of the L1 and L∞ norms for u(x, y) = (0.01x + 2, −0.01y + 2).

(29)

convergence nicely, but if the error is initially very large than that does not matter as much. Figure 11 shows the error in 1 of the 4 components of the advection for both the simple and conservative velocity reconstructions.

Focusing on the ∂uu/∂x error term for the conservative reconstruction, it is clear that the biggest errors occur in cells that are irregular in x direction. In the notation used in sec-tion 4.2.4 this would mean that δx is large. In equasec-tion 21, a variasec-tion in y direcsec-tion was considered, but the linear terms containing δy also contained derivatives of u(x, y) that are 0 in the linear velocity field we are currently considering. If a variation in the x direction is considered, this linear term would be 1₆δxu2

x, which does not contain any zero terms, and

causes these larger errors.

When looking at the error in the ∂uu/∂x term for the simple velocity reconstruction, it is clear that the largest errors occur in very different cells. Firstly, it is worth noting that the errors are much larger in general. Secondly, from looking at figure 11 it is clear that in the case of the simple velocity reconstruction, cells with similar shapes show similar errors. The errors are largest in the cells where ∆x_∆y is either very small or very large. This is exactly what was predicted based on the mathematical analysis in section 4.2.4.

0 20 40 x 0 20 40 60 80 100 y 0.02 0.04 0.06 0.08 0.1 0 20 40 x 0 20 40 60 80 100 y 0.5 1 1.5 2 2.5 3 ×10-5

Figure 11: The absolute error in the ∂uu/∂x term for the simple (left) and conservative (right) velocity reconstructions with a linear velocity field.

Figure 12 shows the relative errors made in the x-advection term as a whole. It is clear that the conservative velocity reconstruction method is much more accurate.

Figure 13 shows the same ∂uu/∂x term that was shown in figure 11, after the mesh has been refined once. The behavior is as expected. The conservative method shows smaller errors

(30)

0 20 40 x 0 20 40 60 80 100 y 0.2 0.4 0.6 0.8 1 1.2 1.4 0 20 40 x 0 20 40 60 80 100 y 1 2 3 4 5 6 7 8 9 ×10-4

Figure 12: The relative error in the x-advection term for the simple (left) and conservative (right) velocity reconstructions with a linear velocity field.

on cells that are created within what used to be the continuity volumes, and slightly larger errors on the faces that previously featured larger errors. The simple velocity reconstruction still shows similar errors in similarly shaped cells, and while the errors did get smaller, they also increased in number at the same rate, thus not converging.

(31)

0 20 40 x 0 20 40 60 80 100 y 0.05 0.1 0.15 0.2 0.25 0.3 0 20 40 x 0 20 40 60 80 100 y 0 0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 ×10-5

Figure 13: The absolute error in the ∂uu/∂x term for the simple and conservative velocity reconstructions, on a mesh that has been refined once, with a linear velocity field.

6.2 Quadratic Velocity Field

The next step is to consider a quadratic velocity field. Here the conservative velocity recon-struction will no longer be exact. In the final stages of the study, we found out that because of the way the velocity was discretized it will also no longer conservative. It is unknown exactly how exactly this effects the result, but it should not change the order of convergence. It is likely that correcting the test would only result in improvements to the results. The velocity field that was chosen is u(x, y) = (.001x2 _{− .002xy + 8, .001y}2 _{− .002xy + 8), and}

figure 14 shows a vector plot of this velocity field. This specific field is divergence free, and its values are roughly between 0 and 10. These are properties one could expect to find in an actual water flow, and thus this should provide a good reference for the accuracy in practice. Figure 15 shows the convergence behavior. As expected, the simple velocity reconstruction once again shows constant norms. The conservative velocity reconstruction also behaves as expected. Initially the errors in the most irregular cells play a bigger part in the L1 norm,

but as the mesh is refined more the second-order behavior starts to become dominant. Figure 16 shows the ∂uu/∂x error term for both velocity reconstructions. Once again it is clear that the simple reconstruction makes larger errors than the conservative method. The variation in velocity over the domain is much larger for the quadratic velocity field, and the second derivatives are no longer zero, both of which also introduce a slightly larger variation in errors. Still, when looking at the errors made by the simple velocity field, the shape of

(32)

0 20 40 x 0 20 40 60 80 100 y

Figure 14: A vector plot of the quadratic velocity field u(x, y) = (.001x2−.002xy +8, .001y2₋

.002xy + 8)

the cell does still seem to be a good predictor. The errors also becomes larger when the difference between x and y becomes larger in the top-left and bottom right corner. This is not surprising, since that is the areas in which the advection terms are largest as well. The figure showing the errors for the conservative reconstruction are slightly harder to inter-pret because there is a very large error in the bottom-right cell making it harder to distinguish differences in the rest of the mesh. The very large error unsurprisingly occur in a very irreg-ularly shaped cell.

After refining the mesh once there are no real surprises. When zooming in on the bottom right of the conservative velocity field picture, the errors in the irregular cells are larger in the 2 rightmost ones. This suggests that it is indeed the shape of the cell that is largely contributing to the size of the error there. Finally figure 18 shows the relative errors of the entire x-advection term on the unrefined mesh. The relative errors here are significantly larger than in the case of the linear velocity field. The difference between the simple and conservative velocity reconstruction is also once again much larger when looking at relative velocities.

Because almost all of the large errors are made near the boundaries of the domain, where the mesh is irregular and the triangles are not as nicely shaped, it would be interesting to see the result on a triangular domain. The errors would be much smaller because the mesh could be regular, with every cell in the mesh being identical.

(33)

Figure 15: The convergence of the L1 and L∞ norms for u(x, y) = (.001x2 − .002xy +

8, .001y2− .002xy + 8). 0 20 40 x 0 20 40 60 80 100 y 0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 1.8 2 0 20 40 x 0 20 40 60 80 100 y 0.1 0.2 0.3 0.4 0.5

Figure 16: The absolute error in the ∂uu/∂x term for the simple (left) and conservative (right) velocity reconstructions with a quadratic velocity field.

(34)

0 20 40 x 0 20 40 60 80 100 y 2 4 6 8 10 12 0 20 40 x 0 20 40 60 80 100 y 0.05 0.1 0.15 0.2 0.25 0.3 0.35 0.4

Figure 17: The absolute error in the ∂uu/∂x term for the simple (left) and conservative (right) velocity reconstructions on a mesh that has been refined once with a quadratic velocity field.

0 20 40 x 0 20 40 60 80 100 y 0.2 0.4 0.6 0.8 1 1.2 0 20 40 x 0 20 40 60 80 100 y 0.02 0.04 0.06 0.08 0.1

Figure 18: The relative error in the x-advection term for the simple (left) and conservative (right) velocity reconstructions with a quadratic velocity field.

(35)

6.3 Sinusoidal Velocity Field

Another type of velocity field that might occur in a real application, is some sort of periodic field. Therefore a good testcase would be some form of sinusoidal velocity field. For this test the field u(x, y) = (sin(y(π/50)), sin(x(π/25))) is chosen, and figure 19 shows a vector plot of this velocity field. The coefficients are chosen such that there is exactly one period of the sine in both directions of the domain. This continuous velocity field is also divergence free, and the values of the sine, between -1 and 1, are within the range of what could be expected in a water flow.

Figure 20 shows the convergence behavior for this velocity field. This is arguably the most interesting convergence plot so far, as the L∞ norm for conservative velocity reconstruction

actually increases on the first refinement. On the second and third refinement however, the L∞ norm quickly converges to first-order behavior. The L1 norm shows the expected

behavior.

The cell that causes the strange behavior in the L∞norm is located in the top left, roughly at

(10,90). This can be see in figure 22, which shows the errors after one refinement, after some zooming in, since the cell is very small. When looking at eq. 20, all the positive term depend on ∆y, and the only negative term depends on ∆x. If either ∆x >>> ∆y, or ∆y >>> ∆x, this error could be large. Due to the way the refinement is done, it is possible that during the first refinement cells are created in which the relative difference between ∆x and ∆y is larger than in any of the cells before refinement. This could be contributing to this increase in the L∞ norm. Another factor that could be contributing to the error here is the fact that

the velocity field does not satisfy the discrete continuity equation.

When looking at figure 21, we see relatively large errors, but still see the conservative velocity reconstruction doing better than the simple velocity reconstruction.

Looking at the simple velocity reconstruction here, both in figure 21 and 22, there seems to be some periodicity to the error. The errors are smallest when y = 0, y = 25, y = 50, y = 75 and y = 100. These correspond to the points where either u(x, y) = 0, namely y = 0, y = 50 and y = 100, or the points where uy(x, y) = 0, namely y = 25 and y = 75.

(36)

0 20 40 x 0 20 40 60 80 100 y

Figure 19: A vector plot of the sinusoidal velocity field u(x, y) = (sin(y(π/50)), sin(x(π/25)))

(37)

0 20 40 x 0 20 40 60 80 100 y 0.01 0.02 0.03 0.04 0.05 0.06 0.07 0.08 0.09 0 20 40 x 0 20 40 60 80 100 y 2 4 6 8 10 12 14 16 18 ×10-3

Figure 21: The absolute error in the ∂uu/∂x term for the simple (left) and conservative (right) velocity reconstructions with a sinusoidal velocity field.

0 20 40 x 0 20 40 60 80 100 y 0 0.05 0.1 0.15 0.2 0.25 0 20 40 x 0 20 40 60 80 100 y 0.005 0.01 0.015 0.02 0.025 0.03

Figure 22: The absolute error in the ∂uu/∂x term for the simple (left) and conservative (right) velocity reconstructions, on a mesh that has been refined once, with a sinusoidal velocity field.

(38)

0 20 40 x 0 20 40 60 80 100 y 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0 20 40 x 0 20 40 60 80 100 y 0.02 0.04 0.06 0.08 0.1 0.12

Figure 23: The relative error in the x-advection term for the simple (left) and conservative (right) velocity reconstructions with a sinusoidal velocity field.

(39)

6.4 Poiseuille Velocity Field

The Poiseuille velocity field is interesting because it shows a weakness of the conservative method relative to the simple method. Figure 24 shows the Poiseuille velocity field. The conservative method makes an error in the ∂uu/∂x term because the velocity reconstruction is not exact due to the quadratic term in v. The simple reconstruction is exact in ∂uu/∂x-direction, because there is no y information in that computation. Figure 25 shows this difference nicely. Other than that, the behavior is what would be expected. Second-order L1 and first-order L∞ convergence with the conservative reconstruction, and no convergence

(40)

0 20 40 x 0 20 40 60 80 100 y

Figure 24: A vector plot of the Poiseuille velocity field u(x, y) = (0, 0.002x(x − 50))

0 20 40 x 0 20 40 60 80 100 y -1 -0.5 0 0.5 1 0 20 40 x 0 20 40 60 80 100 y 0.5 1 1.5 2 2.5 3 3.5 4 4.5 ×10-4

Figure 25: The absolute error in the ∂uu/∂x term for the simple (left) and conservative (right) velocity reconstructions respectively with a Poiseuille velocity field.

(41)

6.5 Comparison to D-FLOW FM

A direct quantitative comparison of the current results to D-FLOW is unfortunately not possible, because the required results are not present. However Visee [31] performed extensive experimentation on the error in the advection terms that D-FLOW makes. She also looked at this same error when higher-order discretizations were applied to the advection, while trying to stay very close to the current D-FLOW method. Her results show linear convergence, or even no convergence at all, depending on the method, for both linear and quadratic velocity fields. So while the magnitude of the errors cannot be directly compared, we can certainly conclude that the method detailed in this thesis has better convergence behavior than the previously tried improvements.

(42)

7 Conclusion

D-FLOW Flexible Mesh is a state-of-the-art hydrodynamics simulation package developed by Deltares. It is used all over the world, and also intended to be the successor to the currently operational hydrodynamics packages used by Rijkswaterstaat in the Netherlands. Traditionally Cartesian and curvilinear meshes were used, but they could not provide the required modelling flexibility. In order to achieve the flexibility their users desired, Deltares started to work on D-FLOW FM about a decade ago.

The discretization in D-FLOW FM is optimized for situations in which the pressure gradient terms are dominant, like tidal wave modelling. Unfortunately, as a result of this, when advective terms are dominant the advection discretization tends to cause errors on highly irregular meshes that can be too large to be acceptable. It is possible to reduce these errors by means of mesh redesign, but this brings about undesirable additional work. A discretization of the advection that is less sensitive to mesh irregularities would be preferable. Small attempts at improving the advection discretization while making minimal changes to the rest of the numerical scheme have so far had limited success. This thesis details an alternative and compact discretization of the advection term. The method relies on having the full velocity vector as the unknown per cell faces, instead of just the normal component as is presently the case. The finite volume method is used. The control volumes for the momentum equation are centered around the faces on which the velocity vectors are computed. The control volume is created by connection the two circumcenters adjacent to the face, to the two mesh nodes at the ends of the face.

Two methods to reconstruct the velocity within the cell are tested. A very simple linear interpolation proves to be too inaccurate. The second approach is that velocities within the cells are reconstructed using a method that ensures conservation of mass over any arbitrary volume. This conservative velocity reconstruction method is exact for linear velocity fields, and second-order accurate for higher order velocity fields.

Algebraic analysis predicts first-order errors that only depend on the irregular aspects of the mesh, and second-order errors on regular meshes. Computed values for the advection are compared to analytic values, and the results look promising, but are based on a central discretization. This cannot be applied in practice, because the dissipation from upwind discretizations is required to ensure stability. On the coarsest meshes that were tested, relative errors of the advection vary between around 1% for a typical cell, and 10% for the most irregular cell. The convergence behavior of the method is good. The L1 norm

has second-order convergence and the L∞ norm has first-order convergence. The errors are

larger in cells that are irregular. There are two types of irregularity in the diamond-shaped momentum volumes. The first type is the difference in length between its diagonals ∆x and ∆y, which leads to very long and thin control volumes. The second type of irregularity occurs when one of the diagonals does not intersect the other diagonal in its center. This creates a kite-shaped volume rather than a diamond-shaped volume. The difference in length between the two segments of the diagonal is called δx or δy. When the velocity field has a non-zero

(43)

second derivative, the dominant error terms depend on ∆y/∆x, which represent how narrow a cell is. If the velocity field is linear, the dominant error terms scale with δx or δy, which represent how far from an ideal diamond shape the cell is.

7.1 Future Research

The results were acquired omitting any variation in water depth, and without any kind of upwind discretization. Before any decisions on the implementation of this method into the D-FLOW FM code are made, more extensive testing is required. The experiments conducted in this thesis did not consider the discretization near the boundaries of the domain. When-ever staggered schemes are used, the handling of boundaries is not trivial, so the boundary discretization of the advection scheme and velocity reconstructions in this thesis should be studied. This study only looked at a velocity field reconstruction within the continuity cells that is based on the assumption that the velocity reconstruction has to fit the constraints of the continuity equation. This is by no means a proven fact, and it would be relevant to see if a more straightforward bi-linear interpolation could provide similar or perhaps better results.

Another limitation of this study is that it focused purely on the computation of the advection term, without considering the complete discretization of the shallow-water equations in both space and time. The premise was that most of the discretization as it is currently used in D-FLOW FM should be maintained. However, the method detailed in this thesis also requires the tangential velocity component on every cell face to be known. The extra velocity information that is available by having the full velocity vector at the cell faces, can be used to compute the diffusion term in a manner similar to the advection term. The main challenge will be the discretization of the pressure gradient term for the tangential velocity component, which is yet to be investigated.

Improving the Advection Discretization in D-FLOW Flexible Mesh