Lecture Notes Numerical Hydraulics
5-1
5 Numerical Methods
This chapter is a general introduction to numerical methods. It explains the main features of the different approaches and defines the important terms and concepts. This should also give a basis for the correct choice and the proper use of commercial software.
5.1 Introduction
Due to the rapid developments in the field of computer technology, computing capacities are available to the user today, which just a few years ago were possible only on a super computer. This computing capacity allows us to use numerical solutions, where in the past estimates obtained via simplified approaches had to be used. In the field of flow calculations, the term "Computational Fluid Dynamics" (CFD) is often used. CFD deals with the numerical solution of partial differential equations (PDG) in the area of fluid mechanics, i.e. with the Navier Stokes-, Reynolds- or shallow water equations, among other things. For the solution of these differential equations, mainly three different methods are used: the Finite Differences method (FD), the Finite Volumes method (FV) and the Finite Elements method (FE).
5.2 Discretization
To allow for the numerical solution of the PDGs, in a first step the computational domain in space and time has to be divided into smaller subdomains. This procedure is called discretization. It is necessary to reduce the infinite number of degrees of freedom of continuous equations to a finite number, which can be handled by the computer. The finer the discretization into smaller domains, the smaller is the approximation error. There are different ways to discretize the computational domain. Not necessarily, but often the choice of the discretization method depends on the computational method. The aim of the following paragraphs is not so much to explain the discretization methods in detail, but rather to present the different approaches and outline their advantages and disadvantages. We start with the spatial discretization. We distinguish between structured and unstructured grids. 5.2.1 Structured grids
Structured grids are the simplest grids that can be used for the discretization of a rectangular domain. They can be characterized by the fact that the number of subdivisions in one direction, independently of the number of subdivision in another direction, is always constant, i.e. a 2D grid has a certain amount of cells nx and ny in x and y direction. Therefore, each cell in the computational grid or every calculation node can be identified unambiguously using two indices (for a 2D grid). Structured grids can be further subdivided into: Regular grids, where the cell length in each direction is always constant Orthogonal grids, where the grid lines in the different directions are normal to each other and Curvilinear grids, where the grid lines can be arbitrary curves in space
5-2
Lecture Notes Numerical Hydraulics
Fig. 5-1: Different types of structured grids, regular and orthogonal grids (A and B), curvilinear grid (C) and orthogonal, curvilinear grid (D).
Structured, and especially regular structured grids are very suitable for the FD-method. Thanks to the properties of this type of grid, the derivatives in all spatial directions can be determined well, without a need for coordinate transformations. A great advantage of the structured grids lies in the small storage requirements, since the nodes of a cell and its neighbours are straightforwardly determined using their indices. This is also of advantage for computation times. The laborious calculation of contributions to a node from neighbouring cells is not necessary and the structure of the solution matrix is known in advance. The main disadvantage of structured grids is that they are not flexible enough to discretize complex domains. 5.2.2 Multiblock grids
Multiblock grids are a combination of many blocks, where each individual block is structured. Hence they combine the advantages of structured grids and avoid their rigid flexibility constraints. On each individual block of a multiblock grid the computation is done exactly as it is done on a structured grid. However, the passage from one block to the next requires special treatment.
Fig. 5-2: Multiblock grid consisting of 2 blocks.
5.2.3
Unstructured grids
Unstructured grids are the most complex to generate, but on the other hand they offer the highest flexibility. On an unstructured grid the number of neighbours of a cell is variable. Therefore it is no more possible to identify a cell using two (2D) or three (3D) indices only. Moreover, the nodes
Lecture Notes Numerical Hydraulics
5-3
forming a cell are not known a priori, but they have to be identified. An unstructured grid must contain at least two types of information: All nodes and their coordinates, and all cells and the related information, which nodes constitute each individual cell. Unstructured grids are mostly used in conjunction with the finite elements method. Partly they are also employed in the FV method. Unstructured grids have the disadvantage of requiring significantly higher computation time. They have, however, the advantage that they can be refined locally.
Fig. 5-3: Unstructured grid.
5.2.4
Grid generation
For simple domains, which require only a small number of cells, the grid generation can eventually also be done by hand. For more complex domains with a large number of cells this would be too time consuming. Therefore tools are used, which, for a given domain contour automatically subdivide the domain into cells. Such a tool is called grid generator. Algebraic grid generators create cells by using algebraic relations. Hereby they take different criteria into account, e.g., maximize the cell area, angle criteria, etc. Alternatively there are grid generators, which solve PDGs to subdivide the domain. Depending on the type of PDG used, the grid generators are denominated accordingly, e.g. elliptic grid generators use an elliptic PDF like the Poisson equation. The procedure is that a simple initial grid is transformed over the complex domain by solving a PDG, where the boundary coordinates of the new domain are the boundary conditions of the PDG.
5.3 Boundary and initial conditions
To allow for a numerical treatment, first of all the domain has to be discretized, as explained above. After that, a PDG is solved over the domain with the use of a numerical method. An unambiguous solution of the PDG is only possible provided that boundary conditions are specified. For an unsteady problem, in addition to the boundary conditions, also boundary conditions in time, i.e. the initial conditions have to be specified. The numerical solution of the problem is significantly influenced by boundary and initial conditions, respectively. Even the best numerical algorithm will not help much, when boundary conditions are not specified appropriately, because they influence the solution at the boundary and
5-4
Lecture Notes Numerical Hydraulics
often even deep into the computational domain. In so-called blind tests, which are contests where numerical modellers calculate flow problems without knowledge of the measured data, using the same program, typically very different correspondence with reality was achieved, depending on the ability of the modeller to set the boundary conditions adequately. As for the initial conditions, the situation is somewhat different. Obviously, they have to be specified correctly, especially when the temporal evolution of the flow has to be correct from the beginning. However, the flow forgets the initial condition after some time. Often, the initial condition is of secondary importance, because only a final steady situation is of interest. In this case it is important to set physically meaningful initial conditions, so that the solution does converge, especially for the case of nonlinear flow problems. Often it is then meaningful to start from a known initial condition, which can for example be a flow at rest.
5.4 Numerical methods
5.4.1 Basic principles
Numerical methods for the calculation of flows can be compared to experimental investigations. In both cases only a limited number of flow quantities can be measured or calculated, respectively. However, their number can be chosen arbitrarily high, so that this constraint does not imply any limitation as far as practical questions are concerned. Different from experimental investigations, where primary unknowns are measured via probes in a limited number of points, numerical techniques prescribe a discrete number of grid points for the unknowns and solve the algebraic equations using a numerical algorithm. Since the focus is on variables to be determined at a finite number of nodes, a discontinuous distribution of variables replaces the continuous one. This approach is known as discretization. The algebraic expressions are derived assuming a distribution between the unknown node-variables. The type of distribution depends on the discretization method. Since the assumed distributions are typically simple expressions, their validity is obviously restricted to a small subdomain. Such a subdomain is termed element and the subdivision that leads to such elements is called space discretization, as explained in Section 5.2 above. Up to now we have treated the discretization method in a rather general way. As mentioned, the type of equations used for the approximation of the distribution of the variables within the elements depends on the discretization method. Without outlining this for now in detail, it can be mentioned that for the finite differences method the differential form of the PDGs is used, whereas the finite volumes method uses the integral form and the finite elements method uses the weak integral form. What this means more precisely will be explained in the following on the basis of a concrete example. 5.4.2 Example
We use a simple, one-dimensional problem, namely the Poisson Equation (e.g. groundwater flows), to illustrate the application of the different discretization methods and related equations. The Poisson equation has the following form: 2 = f x 2
(5-1)
where represents, e.g. the velocity potential. In addition we assume homogeneous boundary conditions on both sides, i.e. (x=0) = (x=1) = 0.
Lecture Notes Numerical Hydraulics 5.4.3 Finite Differenzes method
5-5
The finite differences method uses the differential form of the equations, as mentioned in (5-1). For the FD these equations are approximated by a Taylor expansion. Depending on the accuracy requirements, higher order terms are neglected. Starting from the discretization in Fig. 5-4, where point two lies between points 1 and 3, and assuming that all the points are separated by the same distance, we can apply the following Taylor expansion:
Fig. 5-4: Discretization in the finite differences method.
2 1 2 1 = 2 x + x 2 .... x x 2 2 2 2 1 2 3 = 2 + x + x 2 + .... x x 2 2 2
(5-2)
(5-3)
Truncating after the third term and subtracting/adding Eq. (5-2) from/to Eq. (5-3), we get:
3 1 and = x 2 2x 2 1 + 3 2 2 2 = x x 2 2
Hence, in the finite-differences formulation the Poisson-Equation becomes:
+ 3 2 2 f 2 = 1 x 2
(5-4)
(5-5)
(5-6)
Since the Taylor series was expanded only up to the quadratic terms, one can infer that the estimate of the error arising from the space discretization to be on the order of O(x2), or, in other words, the accuracy of the applied discretization is second order in space. Generalizing the formulation for a generic node i, Eq. (5-6) is written as: + i +1 2 i f i = i 1 for i = 2,...n-1 x 2
(5-7)
This relation can be formulated for all nodes of the computational domain. By doing so, we obtain a matrix, which has a tridiagonal form, i.e. it has nonzero elements only in the main diagonal, the first diagonal below this, and the first diagonal above the main diagonal:
5-6
Lecture Notes Numerical Hydraulics
0 0 0 0 2 1 0 1 2 1 0 . 0 0 . 0 0 1 2 1 . 1 A= 0 0 1 2 1 . 0 2 x 0 . . 1 2 1 0 0 . . 1 2 1 0 0 0 0 0 0 1 2
(5-8)
If we denote the solution vector as xT = [1, 2, 3, ...n] and the right-hand-side of the equation as bT = [f1, f2, f3, ... fn], we obtain a matrix equation:
A x = b fr i = 2,...n-1
(5-9)
Including the two boundary conditions at x=0 and x=1, the matrix equation (5-9) can be solved using a standard equation solver. The example showed that the finite differences method transforms the PDG into a system of algebraic equations. In addition, it can be shown that the approximate solution converges towards the real one, when the number of discretization points goes to infinity. 5.4.4 Finite Volumes method
For the finite volume method, Eq. (5-1) is recast in a somewhat different form. We write:
= f x x
(5-10)
Taking the definitions of Fig. 5-5 into account, Eq. (5-10) is integrated over the control volume, which is x long and extends from the western point w to the eastern point e. The integration yields:
= f dx = 0 + x e x w w
e
(5-11)
Similar to the finite differences method before, in the finite volumes method the derivatives are approximated by difference quotients. Note that the derivatives at the nodes w (west) and e (east) lie between the discretization nodes. Using the values at the nodes, Eq. (5-11) is discretized as:
P P W + E x x e w
= fx
(5-12)
where f is the average of f over the control volume. Two important advantages of the finite volumes method have not been mentioned so far: With the use of the integral form (5-11) the conservation of mass, momentum or energy is guaranteed in the respective transport equations, because an explicit balance between inflow, outflow and storage is required for all elements. The balance of these fluxes is in principle independent of the form of the elements, i.e. the finite volume method is much more flexible compared to the finite differences method.
Lecture Notes Numerical Hydraulics
5-7
Fig. 5-5: Finite Volumes discretization of the model problem.
5.4.5
Finite Elements method
We start once more from Eq. (5-1) and use the same spatial discretization employed previously for the other two methods (FD and FV). However, different from the other two methods, we will not use either differential or integral form, but instead we use the so-called weak integral form. The procedure is to apply the method of weighted residuals to the initial Eq. (5-1). For each discretization an error is introduced, i.e. the solution of the discrete system is different from the solution of the initial equation. The method of weighted residuals requires that the error associated with the initial equation multiplied with a weighting function vanish, when integrating over the domain, i.e.:
2 W 2 + f dx = 0 for i=2,..n x xi 1
xi
(5-13)
Moreover, in the finite elements (FE) method it is assumed that the variation of a variable over a space element can be approximated using a polynomial function. In the simplest case this will be a linear interpolation funcitn, i.e. for the variable we assume:
= 1 + 2 x
(5-14)
Now, with the use of this assumption, the values 1 and 2 can be expressed by using the values of the variables at the nodes, i.e. they can be written as functions of 1 and 2. We can write:
1 x1 1 e = 1 = = K 2 1 x 2 2
With Eq. (5-15) written in matrix form and inserting for we obtain:
= [1 x ]K 1 e
(5-15)
(5-16)
In Eq. (5-16), the variation of a variable over an element is related to the discrete values of the variables at the nodes. The quantity:
N i = [1 x ]K 1
(5-17)
is termed basis function. The basis function at a given node of an element must be equal to one at the node itself and must vanish at all remaining nodes of the element. With reference to the definition on Fig. 5-6, we note that the value of the variable at the location x in element 1 can be obtained by summing the values at the nodes at both nodes of this element, multiplied by the respective basis function, i.e.
5-8
Lecture Notes Numerical Hydraulics
x = N 1x 1 + N 2 x 2
(5-18)
Fig. 5-6: Finite Elements diskretisation of the model problem.
The tentative solution is inserted into the integral (5-1). A good choice for the weighting function W is the basis function itself (Galerkin method). Since the basis function chosen here is not differentiable twice, in the differentiation one derivative is removed through partial integration and carried over to the weighting function W. The resulting system of linear equations for the unknown functions at the nodes must be solved using an equation solver. At the nodes at the boundaries, the boundary conditions of the problem have to be specified.
5.5 Time discretization
The discretizations in time and space are typically well separated from each other. Exceptions are space/time elements in the FE-method. Such elements are successfully used for problems, where the boundary of the computational domain varies in time (free water surface, movable boundaries, etc.). For now we will not discuss the different time discretization in major detail, but rather outline the basic possibilities at hand. Fist of all, one distinguishes between explicit and implicit methods for time discretization. In the explicit time discretization, variables at the new time instance depend only on values at the old time step. This means that the values of each variable at the nodes at the new time depend only on known values at the previous time step. In this way no equation solver is necessary and the solution is straightforward and easy to program. However, for too large time intervals the solution becomes instable. In contrast, the fully implicit method is strictly stable. In this method, the values at the new time depend on its neighbouring values at the new time, i.e. a system of coupled equations must be set up and solved. This is also the case for so-called semiimplicit methods, where the variable at the new time depends both on the old and new time step. One can improve the accuracy of the time discretization using semi-implicit methods, however, the big advantage of explicit methods, which is that equation solvers are not required, is lost completely.