Mesh Generation Concept. Structured Grids V. Selmin

1. Mesh Generation Concept. Structured Grids V. Selmin Multidisciplinary Computation and Numerical Simulation

2. Mesh Generation Concept The computational grid on which a discretisation of an equation is a fundamental component of the numerical solution procedure. The general requirement for a grid are: 1- The computational grid must discretise the domain and all the boundaries . This closely links the grid to the geometry of the domain. The distribution of points in the domain must adequately represent the geometry of the region. 2- The nodes and elements of the grid should be inter-related or connected, so that appropriate approximations to mathematical operators, i.e. the discrete approximation, can be computed in a efficient manner. This determines the inherent structure within the grid. 3- The spatial distribution of nodes should provide a well-shaped elements with no large discontinuities in nodal an element spacing. The accuracy of a numerical discretisation is dependent of the properties of the underlying spatial discretisation. For example, the classical second order accurate finite difference representation of a second spatial derivative is only second order if the grid point spacing is uniform. Furthermore, it can be shown that other spatial derivatives only achieve second order accuracy if the grid lines are approximately orthogonal. Hence the spatial discretisation should be achieved without discontinuous grid point spacing and without the introduction of higly skewed cells or elements. 4- The distribution of nodes in the grid should allow for an appropriate resolution of features in the solution. For accurate simulations, it is necessary that there are adequate points in regions where there is high activity, i.e. where the unknown variables change rapidtly. 5- Ideally, the minimum number of nodes should be used to achieve a given level of accuracy of the solution of the governing equations. Mesh Generation Concept

3. Mesh Generation Concept Structured Grids • Mesh structure: • Domain divided into a structured assembly of • quadrilateral cells • Each interior nodal points is surrounded by • exactly the same number of mesh cells (or • elements) • Directions within the mesh can be immediately • identify by associating a curvilinear co-ordinates • system • it is possible to immediately identify the nearest • neigthbours of any node j on the mesh • Advantages: • Large number of algorithms for discretisation are • available • The algorithms can be normally implemented in • a computationallly efficient manner • Disadvantages: • Difficulty to generate grids of regions of general • shapes  multi-block grids • Very high elapsed time necessary to produce a • grid for domains of extremely complex shape Finite difference, finite volume and finite element discretisations Mesh Generation Concept

4. Mesh Generation Concept Unstructured Grids • Mesh structure: • Computational domain divided into an unstructured • assembly of computational cells • The number of cells surrounding a typical interior • node is not necessarily constant • The nodes and the elements has to be numbered • To get the necessary information on the neightbours • the numerotation of the nodes wich belong to each • element has to be stored • The concept of directionality does not exist anymore • Advantages: • Powerful tool for discretising domains of complex • shapes • Unstructured mesh methods naturally offer the • possibility of incorporating adaptivity • Disadvantages: • Alternative solution algorithms are more limited • Computational implementation places large demands • on both computer memory and CPU Finite volume and finite element discretisations Mesh Generation Concept

5. Geometry Modelling The boundary of the domain to be discretised needs to be represented in a suitable manner before the generation can start. If the automatic discretisation of an arbitrary domain is to be achieved, the mathematical description of the domain topology has to possess the greatest possible generality. The computer implementation of this description must provide means for automatically computing any geometrical quantity relevant to the generation procedure. Planar two-dimensional case In the planar two dimensional case, the boundary is represented by closed loops of orientated composite curves. For simple connected domains these curves are orientated in a counter-clockwise sense while for multiple connected regions the exterior boundary curves are given a counter-clockwise orientation and all the interior boundary curves are orientated in a clockwise sense Mathematical description of the geometry → CAD model Geometry modelling

6. Geometry Modelling Three dimensional case In three dimensions, the domain to be discretised is viewed as a region bounded by surfaces which intersect along curves. The portions of these curves and surfaces needed to define the three dimensional domain of interest are called curve and surface components, respectively. Decomposition of the boundary of a three dimensional domain into its surfaces and curve components Geometry modellingt

7. Geometry Modelling In addition, boundary curves and surfaces are oriented. This is important in the generating process as it is used to define the location of the region to be discretised. The orientation of a boundary surface is defined by the direction of the inward normal. The orientation of the boundary curves is defined with respect to the boundary surfaces that contain them. Each boundary curve will be common to two boundary surfaces and will have opposite orientations with respect to each of them. An example of the approximated geometry of the surfaces component and its curve components is here presented. Geometry modelling

8. Curve Representation A support curve may be described by a piecewise parametric representation. In those representation, the curve is subdivided into arcs, and the position vector r of a generic point on each arc is expressed as a function of a single real parameter u, which by convention varies into the interval [0,1]. In general, this function is represented by a polynomial whose rank may change from arc to arc and which can be expressed as where each is a vector formed by three coefficients which represents the components of with respect to a cartesian reference system (x,y,z) , whereas n is the rank of the polynomial. The position vector can be rewritten as follows where and A point on the curve can then be identified by the number of the arc on which it is lying and the value of the parametric coordinates u, which is usually called local parametric coordinate. Geometry modelling

9. Surface Representation The surfaces are represented by using a similar scheme. They are subdivided into patches which forms a regular grid on the so-called parametric plane. On this plane, we can define with respect to each patch the local parametric coordinates u and v with varies on the interval [0,1]. The position r of a node on the surface can then be expressed as a polynomial expansion in u and v on each patch : The total rank of the polynomial is nxm, whereas n and m represent the rank of the polynomial with respect to the cartesian reference system (x,y,z): The position vector can be rewritten in the following form where The quantities U and V are expressed as Geometry modelling

10. Structured Grid Generation and Mapping Concept Consider in a one-dimensional domain , with boundaries and . To generate a set of n+1 equally spaced nodes within the domain, it is convenient to introduce a new variable such that This representation may be generalised to give, for an arbitrary distribution of nodes, where is the domain within which the nodes are uniformly spaced, x is the domain of and the function defines the nodal distribution. In this sense, the generation of a grid is equivalent to finding a function which maps nodes distributed evently in one space (or domain) to the required distribution of nodes in another space. More specifically, the equally spaced nodes in the logical (computational) space are mapped through the grid generation function, to the physical space. The basic concept also applies in two and three dimensions. An arbitrary shape domain is mapped by into a square in two dimensions and by into a cube in three dimensions. The grid in the logical space is uniform and when mapped to the physical space provides the required grid. 0 1 0 L Structured grid generation

11. Mappings and Jacobians Consider the continuous mapping of the form . If a function f exists for which then For any differential function f(x) , it is clear that the change of f with respect to is equal to the change of f with respect to x multiplied by a scaling parameter J which is the change of x with respect to the new variable . Using this equation to solve for gives The term is called the jacobian of the mapping and can be seen to be an effective scaling introduced through the change of variables. The quantity is the ratio of the length in the physical space to the length in the logical space. If the Jacobian vanishes within the domain, then the transformation is invalid In two dimensions, let consider a continuous transformation of the form If we consider a function f with then Structured grid generation

12. Mappings and Jacobians that may be rewritten in the more concise form Solving this equation for the derivative of f with respect to the vector x gives If vanishes the the transformation is invalid. Thus, if the Jacobian is never zero, the inverse function guarantees the existence of the local inverse map: In a similar manner, in three dimensions, we consider a continuous transformation of the form Again consider a function f with then Structured grid generation

13. Mappings and Jacobians There is clearly an infinite number of maps between the logical and physical space. The issue in grid generation is to find a map which has the appropriate properties. The properties which are required to a mapping include • Positive Jacobian determinant, • Smoothness • Orthogonality or near othogonality • Uniformity of cells The following table provides some useful quantities used in grid generation Structured grid generation

14. 0 1 Algebraic and Transfinite Interpolation Interpolation The root of algebraic methods for grid generation are in simple interpolation. In one dimension, the boundaries of a domain are the two end points and . It is natural and intuitive that a grid is formed by explicit interpolation between these boundaries. Linear interpolation in one dimension can be written This represents a linear interpolation between the end points with and providing the linear variation. The more general formulation is with The functions and are called blending functions since they blend the boundary data and into the domain . This interpolation may be generalised by choosing other blending functions, for example functions that allows to control slopes of the grid lines at the boundaries It is clear that such simple techniques cannot be in general be applied to domains in two dimensions which have four general boundary curves. In general, interpolation between curves would lead to grid cross-over. Structured grid generation

15. Algebraic and Transfinite Interpolation Multi-directional Interpolation – Two Dimensions A two-dimensional domain can be represented by four bounding curves, two for which and two with . The boundaries are defined as and the four curves intersect at the four corners of the domain: Structured grid generation

16. Algebraic and Transfinite Interpolation As a starting position, the one-dimensional interpolation scheme can be written for a two-dimensional domain in the form where represents the interpolation derived between the boundaries . Similarly, for the boundary curve , the interpolant is Finally, the following operator is formed The operator is called the Boolean sum of and , and can be written as or This type of interpolation is called transfinite interpolations and the functions are called blending functions. This interpolation may be generalised by choosing other blending functions, for example functions that allows to control slopes of the grid lines at the boundaries Structured grid generation

17. Without boundary slope control With boundary slope control Algebraic and Transfinite Interpolation 2D Example Structured grid generation

18. Algebraic and Transfinite Interpolation Multi-directional Interpolation – Three Dimensions Transfinite interpolation schemes can be extended into three dimensions where a domain is assume to be bounded by six faces and naturaly mapped into a cube with local coordinates . The resulting grid will consist of a structured set of hexahedra. To form the interpolant , three unidirectional interpolation schemes are formed between opposite faces of the domain, namely: Finally the interpolant is formed from Structured grid generation

19. Grid Generation by PDE PDE: Partial Differential Equation The mathematical formulation of may problems in science and engineering involve partial differential equations. Special cases of the two-dimensional second order equation: where a, b, c, d, e, f and g may be functions of the independent variables x and y and of the dependent variable , occur frequently because they often take the mathematical form of a conservation principles in physics. This general equation can be classified according to the coefficients. The equation is Elliptic equations are generally associated with equilibrium or steady-state problems. For example, steady flow of incompressible irrotational fluid, or heat conduction in a domain. The solution of an elliptic equation is a function of x and y which satisfies the partial differential equation at every point in a domain. Problems involving changes of variables with time usually lead to parabolic or hyperbolic equations. A common example of a parabolic equation derives from the theory of heat conduction, when the temperature is a function of time. Hyperbolic problems generally originate from vibration problems, or from problems where discontinuities can persist in time, such as shock waves. Perhaps the simplest hyperbolic equation is the one-dimensional wave equation. Although these different types of equations exhibit quiet different properties, they have a role to play in grid generation. Structured grid generation

20. Grid Generation by PDE Grid Generation using Laplace and Poisson equations An appropriate form of grid generation equation in one dimension can be expressed as subject to the generalised boundary conditions where P is an arbitrary function used to control the grid point distribution. Thus, the grid generation problem involves 1- the specification of the two end nodes, and , 2- the generation of the one-dimensional grid by solving the equation In two dimensions, the equivalent equations are where (x,y) are the coordinates in the physical space, the coordinates in the logical space, P and Q are arbitrary functions of the independent variables . The grid generation problem then involves 1- the specification of the four corner nodes, 2- the generation of one-dimensional grids on the four bounding curves, 3- the generation of the two-dimensional grid by solving the elliptic system of equations Structured grid generation

21. Grid Generation by PDE The elliptic operator naturally extend to three dimensions where the form is The equations are solved for the grid coordinates (x,y,z) as a function of the independent coordinates of the logical space . The functions P, Q and R can be used to control the grid node spacing. The grid generation problem then involves 1- the specification of the eight corner nodes, 2- the generation of one-dimensional grids on the twelve bounding curves, 3- the generation of the two-dimensional grids on the six bounding faces, 4- the generation of the three-dimensional grid by solving the elliptic system of equations For the numerical solution procedure, second order accurate discretrisation of the operators are usually used combined with iterative solution procedure and over-relaxing techniques in order to improve convergence. The iterative solution process requires that initial conditions be prescribed to enable the initial solution to be computed. It is possible to begin the iteration with any suitable conditions. However, an effective procedure is to use transfinite interpolation to form the initial guess of the solution. Structured grid generation

22. Grid Generation by PDE A function that satisfies Laplace’s equation is said to be harmonic which has the following properties: 1- Smoothness It is expected that the solution at any given point to be an average of the solution in surrounding region and, as such, be smoothly varying over the domain. In formals terms, it is possible to prove the mean value theorem for harmonic functions. This states that the value of a harmonic function at a point, is equal to the average of its values over the area of any circle centred on those point. In this case, the average would be defined as 2- Extremum principle It is possible to demonstrate to prove by using the mean value theorem, that if is harmonic, then the maximum or minimum of must occur on the boundary. An immedfiate consequence of the maximum principle is that the solution to the problem with specidied on every boundary, is unique and is indefinitely differentiable, i.e. . For these properties, it is possible to infer the behaviour to be expected of grids generated using the system of Lapacian equations: 1- The nodal distribution in the grids will be smooth and any discontinuities in the boundary data will not be propagated into the domain. 2- The solution will be unique 3- The extremum principle may guarantee a one-to-one mapping between spaces. This is equivalent to ensuring a mapping exists for which the Jacobian J of the mapping is positive everywhere. However, this is not true for the grid generation equations written in the previous Laplacian forms: the solution admits negative Jacobians, which means the grid lines fold, or overlap, creating cells that have a negative area or volume. Structured grid generation

23. Grid cross-over Grid Generation by PDE To overcome this practical limitation of the use of the Laplacian operator for the solution of the grid coordinates, it is necessary to consider the inverse mapping, using x and y as the independent variables and and as the independent variables. In two dimensions, the Laplacians for and with control functions and are where Inverting these equations leads with It is noted that the previous equations are nonlinear, elliptic and coupled in x and y. Unlike the original formulation, which could be solved for x and y , these equations must be solved together. Structured grid generation

24. Grid Generation by PDE The major motivation of the use of the inverse equations is to obtain a mapping which possesses the property of the Jacobian being positive everywhere. Rado’s theorem shows that by requiring the inverse map to be harmonic instead of the map itself, then the mapping is one-to-one and onto and hence exhibits the property of a positive Jacobian everywhere. Certainly, in the exact differential form, the inverse equations provide a mean of generating a grid which is smooth, and does not produce a folding of the grid lines. It is emphatised that this property is exhibited by the exact differential form of the equations. This property does not necessarily hold for the numerical solution. Moreover, the advantageous property of no grid folding exhibited by the inverse equations does not in general hold true for the inhomogeneous form of the equations. Rados’ theorem no longer applies. The property of non-folding of grid lines does hold true provided the source term do not change sign within the domain. Grid control functions The inherent smoothing properties of Laplace equations ensures, that, in absence of boundary curvature, the mesh points are evently spaced. However, near convex boundaries the grid points will become more closely spaced, whilst near concave boundaries the mesh spacing will be more sparse. Structured grid generation

25. Grid Generation by PDE The properties are not always desirable for mesh generation where it is essential that the mesh near a boundary reflects the shape of that boundary. Control of mesh point spacing can be achieved by introducing source terms into Laplace equation thus converting the system of equations into Poisson type. The source terms, often call grid control functions, can be so designed to provide the flexibility required to control mesh points spacing. Negative values of Q tend to cause the -coordinates lines to move in the direction of decreasing , while negative values of P cause the -lines to move in the direction of decreasing . Thomson at al have reviewed the many forms of functions which can be used to attract or repel the mesh to lines or points. A typical function is and similar for Q. Here the effect of the amplitude of the is to attract -lines towards the lines while the effect of the amplitude of is to attract -lines toward the single point . Structured grid generation

26. Grid Generation by PDE Grid Generation using Hyperbolic equations Computational grids can also be generated by solving hyperbolic partial differential equations which are marched outwards from a specified distribution of points on an inner boundary. Constraints on mesh size and angle are generally used to determine the governing grid generation equations. Because the outer boundary of a grid can not specified in hyperbolic mesh generation, the method is best suited to applications which have asymptotic like outer boundaries. In external aerodynamics, for example, the outer boundary is in general a far-field condition of undisturbed flow. However, domain decomposition methods that use overset composite meshes for complex field simulations can also effectively employ hyperbolic grid generation. In two dimensions, the idea is to march outward from a user specified distribution of nodes on an inner boundary. The marching is performed by using a hyperbolic systems of equations, in which one of the independent variables corresponds to the marching direction. It is assumed that the boundary is defined at Following Steger, two equations are chosen which represent constraints on mesh orthogonality and mesh cell area. The orthogonality constraint is or and the mesh area or where A is the area of any one mesh cell and must be prescribed Structured grid generation

27. Grid Generation by PDE These two equations form a system of non-linear partial differential equations with initial data specified on the boundary . Steger proposed to solve these equations using a linearisation procedure, leading to with If is the marching direction, then Using the fact that the inverse of a matrix H is where is the determinant of H , and G is the transpose of the matrix of co-factors of H , then it can be shown that The matrix exists if , and the matrix is symmetric. A symmetric matrix has real eigenvalues, so the system is hyperbolic with used as the marching direction. Structured grid generation

28. Grid Generation by PDE The numerical solution is achieved by using backward differences in the -direction and central differences in the -direction, i.e. This gives the marching equations These equations can be solved by either an explicit or implicit finite difference schemes. The implicit scheme can lead to marching the equations for any arbitrary step interval in the -direction. The grid is controlled by the initial nodal spacing on the boundary and by the cell area A , which is clearly a function of x and y, A=A(x,y). The cell areas, which provide the grid control, must be specified at every node in the mesh and as such it is necessary to devise automatic techniques to specify this information. Two approaches may be identified: 1- Use a reference (background) grid where the area distribution is specified 2- Determine the cell area as the grid is generated (extrapolation). An example of a mesh around an aerofoil is here shown. The discontinuity propagated from the trailing edge is clearly visible The advantage of a hyperbolic mesh generation procedure are that orthogonal grids which can offer good user control can be generated very rapidly. The disadvantages are that boundaries in the field, including outer boundary can not be specified and discontinuities in the boundary data can be propagated into the field. Structured grid generation

29. Grid Generation by PDE Grid Generation using Parabolic equations Grid Generation based upon partial differential equations can also be performed. In general, such equations are constructed by modifying elliptic generation systems so that second derivatives, in one co-ordinate direction, do not appear. The solution of the parabolic system can then be marched away from a boundary in much the same way as that for hyperbolic grid generation. Here, however, some influence of the other boundaries is retained in the equations. Although some development of such schemes has been made, it has proved to be limited in its flexibility and applicability. Structured grid generation

30. Multiblock Mesh Generation It is clear that structured meshes result in a set of curvilinear coordinates. One way to seeing this is that, in two dimensional space, the physical domain is mapped to a rectangle in the logical (computational) domain. That means that a cartesian mesh in this logical space is mapped through a transformation to the physical space to provide a mesh consisting of quadrilaterals which are not necessary regular. The global property of curvilinear coordinates makes structured mesh generation somewhat limited in its application to complex geometries. The idea behind multiblock mesh generation is that, instead of using one global curvilinear coordinate system, several local curvilinear systems are constructed and connected together. The domain is subdivided into blocks and within each block a curvilinear system is derived. The block subdivision provides the necessary flexibility to construct structured meshes for complex geometrical shapes. The multiblock concept is not restricted to any particular mesh point generation technique. The generation of points can be performed using the algebraic or the elliptic partial differential approach. Structured grid generation

31. Multiblock Mesh Generation a b Block structure for C-mesh topology Topology Definition: Central to the multiblock concept is the block subdivision of a domain. Here a block is defined as a geometrical domain which is topologically equivalent to a rectangle in two dimensions or a cube in three dimensions. A trivial example of a block subdivision of a domain around an airfoil is given in the following figure. Here the mesh topology is a so-called C-mesh. The block subdivision for the mesh topology is not unique. Two equally valid subdivisions are denoted by a and b. In the case a, the aerofoil boundary is restricted to one block side whilst in the second case the aerofoil boundary extends over one complete side and partially over two others. Either are possible but for discussion purposes it will be assumed that only one type of boundary can exist on any block side. This assumption does not in practice restrict the applicability of the multiblock approach, but does provide some simplification in the implementation. Structured grid generation

32. Multiblock Mesh Generation Block structure for O-mesh topology Block structure for H-mesh topology Another two examples of block subdivisions are shown for a O and a H type grid. The block structure for a polar or O-topology can be reduced to a single block in which two opposite sides map onto each other. The cartesian or H-topology mesh around an aerofoil is given by six blocks suitable connected. Structured grid generation

33. Multiblock Mesh Generation Let consider now the more complex geometry of two aerofoils, with rounded leading edges, arranged in a tandem configuration. It is possible in this case to generate a mesh which effectively has one global coordinate system. A cartesian mesh could be constructed local to each airfoil or alternatively the forward aerofoil could be favoured with a C-mesh but the aft-element would then have a cartesian mesh. On the grounds of mesh quality and flow sensitivity in the vicinity of the leading edges a cartesian mesh local to leading edge may be found to be unsuitable. Ideally, it can be argued that a C-mesh or a polar mesh is required to provide the necessary grid quality for accurate flow simulations. If this is the case, it is no longer possible to build such mesh topologies and conserve a global curvilinear coordinate system. However, it is topologically possible to construct such meshes and the concept of a multiblock subdivision of the domain provides the appropriate mechanism. C-mesh around the leading airfoil and H-mesh around the second airfoil Cartesian mesh around both airfoils Structured grid generation

34. Mesh Generation within a Block Each block has its own curvilinear coordinate system. This system is independent of those in the adjacent blocks. To generate a mesh within a block, it is necessary to construct an algorithm which operates on a rectangular computational domain. With a suitable treatment of boundary points, the generation of the global mesh then involves a loop over all blocks. The mesh within each block is in general build after discretisation of its boundaries. Thus, assuming that the boundary points are known, the interior points may be determined by using a system of elliptic partial differential equations - in this case, the boundary points location corrisponds to a Dirichlet condition – or by using transfinite interpolations. The basic idea is that the boundary points are either known or compute by some other techniques. Structured grid generation

35. Point and Slope Continuity at Interfaces Clearly, a computational mesh can not contain holes or void regions. For this condition to be satisfied, it is necessary that given two blocks which are adjacent and have a common edge, the grid points in physical space, which correspond to points on the edge, must be continuous. With fixed data points on block boundaries, it is possible to form a well defined problem for the generation of points within a block. But, unless special treatment to be applied, the grid lines extending from such boundaries will not possess special angular relationships with boundary. If a block boundary is adjacent to another block boundary, and the boundary points are fixed, the resulting mesh will in general be discontinuous in gradients (and higher derivatives) accross the boundary. It is possible to achieve slope continuity at block interfaces, by providing to the internal nodes grid generation technique, grid lines discretisation and grid lines intersection angles on both block edges. Structured grid generation

36. Singular Points Treatment Flexibility for structured multi-block grid generation introduces singular points in the mesh. Usually in two dimensions two families of lines pass through a given point. Singular points are those points in the field which more or less than two families of curvilinear lines intersect. Classical discretisations, like finite difference formulations cannot be defined or have to be adapted at singular points. Alternative means of expressive derivatives at those points must be found. → Impact on the solver Structured grid generation

37. Example of Multiblock Structured Grid Structured grid generation

38. Summary • Introduction of the mesh generation concept. • Mathematical representation of the geometry. • Structured grid generation • Algebraic and Transfinite Interpolations • Grid generation by using PDE (elliptic, hyperbolic, parabolic) • Multi-block structured grid generation concept • Structured grid generation is efficient for simple geometry, but lack of efficiency for • complex geometry. • Strong impact on the solver that has to be adapted at particular features of the grid • like singular points. • New developments (Chimera techniques, ….) transfer the complexity from the grid generation process • to the solver: easier(?) grid generation process, but more complex solvers. • The technologies based on the unstructured grid generation approach seems to be more adapted for application • of CFD (Computational Fluid Dynamics) around complex shapes Structured grid generation