ASIA-PACIFIC JOURNAL OF CHEMICAL ENGINEERING Asia-Pac. J. Chem. Eng. 2011; 6: 110–119 Published online 9 July 2010 in Wiley Online Library (wileyonlinelibrary.com) DOI:10.1002/apj.479 Special Theme Research Article Icosahedral-hexagonal grids on a sphere for CFD applications H. C. Upadhyaya, Om P. Sharma,* R. Mittal and H. Fatima Centre for Atmospheric Sciences, Indian Institute of Technology Delhi, New Delhi 110016, India Received 13 August 2009; Revised 17 May 2010; Accepted 18 May 2010 ABSTRACT: An efficient method for achieving the non-uniform grids on a spherical geometry is presented here. It uses the icosahedral–hexagonal grid with successive refinement to arrive at the target grid suitable for computational fluid dynamics (CFD) applications. It reduces the complexity of search from O(n 2 ) to O(n) on icosahedral–hexagonal grids where n refers to the refinement level. The numerical solution of the transport equation is performed using the initial conditions of a well-known problem (solid-body rotation); the numerical scheme is second-order accurate. Since advection of chemically active tracers in the atmosphere and their modelling are becoming a major area of concern within the climate change scenario, this study assures the efficiency and accuracy of the numerical scheme as tested here for tracer transport. 2010 Curtin University of Technology and John Wiley & Sons, Ltd. KEYWORDS: CFD; Delaunay-Voronoi association; tracer transport INTRODUCTION A large number of computational fluid dynamics (CFD) problems require numerical solutions of partial differential equations on spheres or in spherical geometries. In chemical engineering problems too, spherical shapes are one among the typical geometrical configurations on which numerical solutions are especially sought for the convection-diffusion equation.[1] The present day computing resources allow very fine grid implementations of local approximations using finite difference, finite element or finite volume methods to maintain the required order of accuracy of the numerical solutions. It is thus desirable to use a grid system which is sufficiently adaptive in nature to exploit the potential of parallel processing and to resolve accurately the developing steep gradients in evolution-type CFD problems in science and engineering. A geometrically ideal grid system for a sphere is the one which divides it into elements of equal area and equal shape, retains homogeneity and preserves isotropy of the sphere around any grid point as much as possible. The regular surface coordinate grid is most commonly used for solving geophysical flows (referred to as lat–lon grid). But it suffers from computational complications associated with the convergence of meridians at *Correspondence to: Om P. Sharma, Centre for Atmospheric Sciences, Indian Institute of Technology Delhi, New Delhi 110016, India. E-mail: opsharma@cas.iitd.ac.in 2010 Curtin University of Technology and John Wiley & Sons, Ltd. Curtin University is a trademark of Curtin University of Technology the poles.[2] As an alternative to lat–lon meshes, quasiuniform grids on the earth’s surface use regular polyhedrons which are topologically equivalent to a sphere.[3] In this class, icosahedral–hexagonal grids[4] and cubed sphere[5] deserve special mention, as they both constitute promising grid structures. This paper presents the details of a grid system designed from icosahedron geometry and aims to present an efficient computational framework for designing numerical schemes on spherical surfaces. Since these meshes are capable of creating non-uniform grids, often this grid system finds some interest in calculating accurately simultaneous heat and mass transfer from a single sphere in a turbulent stream.[6] In recent atmospheric modelling studies, highresolution coupled general circulation models (GCMs) are being designed using the icosahedral–hexagonal meshes for solving atmospheric, oceanic, and chemical and aerosol processes on a single platform.[7 – 9] The modelling of the advection of chemically active tracers in the atmosphere has become a major area of concern with the increasing awareness of climate change. Studies have shown that ad hoc techniques of tracer convection are retro-fitted into atmospheric models. This, however, removes consistency during advection of tracers between mass (and momentum) and the chemical tracer concentration, causing significant errors. Most of the atmospheric constituents undergo gas phase and aqueous phase reactions in the atmosphere and their Asia-Pacific Journal of Chemical Engineering ICOSAHEDRAL–HEXAGONAL GRIDS FOR CFD APPLICATIONS distribution depends upon horizontal and vertical transport of industrial emissions released from the ground. In order to deal with this important chemical engineering component in the evolution of various constituents of the atmosphere, a consistent scheme, like the one proposed by Miura,[10] has been tested in this paper for atmospheric tracer transport. The success of this grid system in atmospheric circulation modelling may prove it equally a promising candidate for numerical simulations in other areas of modelling irrespective of the science discipline. This paper has been organised as follows: first, a complete description of the mathematical modelling framework based on icosahedral–hexagonal grid is presented. Next, the main features of indexing, storage search, and static or dynamic adaptability on this geometry are given, and finally this computational framework is tested for a transport problem that arises in several chemical and environmental engineering applications besides weather and climate predictions. ICOSAHEDRAL–HEXAGONAL GRID GENERATION There are several ways to construct a discretization on the surface of a sphere from the icosahedron geometry. The approach first given by Sadourny et al .[4] has been followed here. The starting grid configuration consists of 12 grid points and 30 great circle arcs. On connecting the nearest neighbours, the spherical surface subdivides into 20 main spherical equilateral triangles. The angles between two sides of the triangles are 2π/5 or 72◦ and the length (ω) of any side of these triangles on a unit sphere follows from Fig. 1 as cos(ω/2) = cos(π/5)/ sin(2π/5) = 1/2 sin(π/5) ⇒ ω ≈ 1.10715 (1) Beginning from this grid consisting of 20 icosahedral triangles (referred to as the zero-level triangulation in the text), a newer level of refined grid of triangles is generated in the manner as described by Sadourny et al .[4] The procedure generates a level-n resolution mesh when each side of triangle in the zero-level triangulation is subdivided into n parts. Figure 2 gives the icosahedral–hexagonal mesh at different resolutions (level-6 and level-10), illustrating the maintenance of isotropy as resolutions becomes finer and finer. At level-n resolution there are n − 1 node points on each side of the triangle and {(n − 1)(n − 2)/2} points within each of 20 triangles. This gives a total number of {12 + 30 (n − 1) + 20 (n − 1)(n − 2)/2} = {10n 2 + 2} points (referred to as np points in the text). The total of n 2 smaller triangles (nearly equilateral) are generated within each level-0 triangle. This gives a total of 20 n 2 triangles (referred to as nt triangles in the text) on the sphere and there are 30 n 2 smaller great circles’ arcs (referred to as ne edges in the text) constituting the sides of these small triangles. Length (L) of a side of zero-level discretization on the spherical earth is equal to ωRE ≈ 7053.589 km, where RE is the radius of the earth. Thus, the arc length or the grid size at a resolution level-n is ≈ ωRE /n, which implies that as n increases, becomes smaller and smaller. Thus, desired levels of fine discretizations on the sphere can be accomplished by choosing an appropriate value of n. With regard to numerical accuracy on the sphere, there exist a number of variants to icosahedral– hexagonal grid system. For example, one may adopt a recursive grid generation method to obtain higher level refinement by dyadic divisions of the original zero-level grid.[11] It has serious restrictions on the flexibility of target grid resolution as pointed out by Giraldo.[12] Another derived grid system is constituted of ‘twisted icosahedral–hexagonal grids’,[13] wherein symmetric grid points across the equator are obtained by simply rotating all the points in the southern hemisphere by 36◦ . They further modified the grid system after refinement in order to further minimize error of numerical differential operators. Tomita et al .[14,15] proposed another sort of grid modification by applying the spring dynamics. They could generate a more homogeneous grid system in comparison to the standard icosahedral grid system by tuning the natural spring length. By setting the natural spring length appropriately, one can obtain a grid system which is most efficient in terms of the Courant–Friedrichs–Levy (CFL) condition. Here, the modelling framework based on the original method of grid generation of Sadourny et al .[4] is further developed. The comparative statistics available in the literature on lat–lon and geodesic grids has been utilized to check the accuracy of the grid generation procedure. GRID STRUCTURE: GEOMETRY, INDEXING, STORAGE AND LOCAL REFINEMENTS Figure 1. One main spherical triangle at the North Pole.[4] . 2010 Curtin University of Technology and John Wiley & Sons, Ltd. The icosahedral–hexagonal grid modelling framework in spite of being very promising remained dormant for a long time, and the foremost reason being utterly Asia-Pac. J. Chem. Eng. 2011; 6: 110–119 DOI: 10.1002/apj 111 112 H. C. UPADHYAYA ET AL. Asia-Pacific Journal of Chemical Engineering Figure 2. Icosahedral–hexagonal grid at resolutions: (a) level-6 and (b) level-10. This figure is available in colour online at www.apjChemEng.com. insufficient computer power (processing speed, memory, and storage) which is mandatory for numerical algorithms and methods for resolving the partial differential equations (PDEs) on this grid. However, it has generated new interest in the researchers due to the availability of vast computing resources at relatively affordable cost than in the past. Grid geometry Geodesic grids use triangles constructed from gridpoints to obtain the Delaunay mesh, whereas Voronoi polygons that surround the Delaunay mesh points are generated by joining the centroids of the neighbouring triangles.[16] These polygons are hexagons or pentagons but Delaunay mesh cells are all nearly equilateral triangles (Fig. 3). From the duality of Voronoi polygons and the Delaunay triangulation, grid refining by successive triangulation on the sphere could be accomplished in a precise manner. Thus, there are nt number of Voronoi points (VPs) and np number of Delaunay points (DPs). From the illustrations given in Figs 3 and 4, it may be noted that each DP is surrounded by 6 neighbouring VPs, except for the 12 fixed DPs that are surrounded by 5 VPs, while all the VPs are surrounded by only 3 DPs. These associations form the key to design adaptive grids for producing non-uniform resolutions on the sphere for achieving greater accuracy of numerical solutions in regions whenever the situation demands. The procedure of laying an adaptive grid over any region of interest on spherical surfaces could be made automatic following VP–DP associations which make this grid system very attractive in industrial applications for 2010 Curtin University of Technology and John Wiley & Sons, Ltd. Figure 3. Schematic diagram of Delaunay-Voronoi association: dashed line represents the Delaunay triangulation and solid lines represent Voronoi polygons. The solid circles represent the basic gridpoints (DPs) and solid hexagons around DPs represent their area of influence. DPs and VPs are labelled as Ps and Qs, respectively. processes which undergo on a single sphere or inside an annular region of concentric spheres. It is quite likely that hanging nodes shall invariably appear in the automatic generation of grids but it can be easily handled as explained in the later part of this paper. Asia-Pac. J. Chem. Eng. 2011; 6: 110–119 DOI: 10.1002/apj Asia-Pacific Journal of Chemical Engineering ICOSAHEDRAL–HEXAGONAL GRIDS FOR CFD APPLICATIONS Figure 4. Labelling of level-0 triangulation: North Pole (NP) at node 1 and South Pole (SP) at node 12. Triangle and node numbering Various conventions for this process have been suggested earlier because much of the efficiency of numerical computations on this grid is linked to numbering of nodes. Further, node numbers are subsequently used for generating and numbering the triangles. Faster retrieval of information during the computing process moreover depends a good deal on how the indices of nodes and triangles have been stored. At the outset, the lowest-order object of a grid is a vertex; next comes in an edge that conveys the connectivity of mesh points along with indices of its starting and terminating vertices. A face represents the interface area between adjacent volumetric cells and could be identified from edges that bound it, or from a sequence of its corner vertices. This kind of systematic indexing of various objects associated with this grid system avoids undesired latency in numerical computations. As a matter of fact, information on the following basic quantities is required to be stored: 2010 Curtin University of Technology and John Wiley & Sons, Ltd. (i) Nodal Points: Neighbouring nodal points; Common vertices of triangles; Edges connecting them. (ii) Triangle Centroids: Neighbouring centroids; Vertices of triangles; Constituting edges of triangles. (iii) Edges: Edge neighbours; Edge extremity nodes (vertices) etc. Node numbering The starting configuration of zero-level icosahedron has 12 vertices, 20 triangles and 30 edges, which are illustrated in Fig. 4 with appropriate labelling. At a resolution level-n, there are a total number of {10n 2 + 2} nodes and their indexing is done using the following procedure: (a) first, index the 12 nodes of the icosahedron (level0); (b) then create indices for nodes lying on 30 edges; and (c) finally, number the nodes located inside the 20 triangles. The numbering of nodes in triangle T1 that one will arrive at, in particular, at a resolution level-4 has been depicted in Fig. 5. Asia-Pac. J. Chem. Eng. 2011; 6: 110–119 DOI: 10.1002/apj 113 114 H. C. UPADHYAYA ET AL. Asia-Pacific Journal of Chemical Engineering Figure 5. Node numbering of Triangle 1 at level-4 resolution. Thus, one may notice that the number of nodes in each icosahedron triangle at a resolution level-n will be (n + 1)(n + 2) . Now for each of the main icosahedral 2 triangles, the node indices can be stored in a lower triangular part of a [(n + 1) × (n + 1)] matrix (abbreviated as LTM) in the following manner. To fix the idea, let us consider the Triangle T1 in Fig. 4 whose nodes have been numbered with the index of each node as shown on this figure. The procedure to fill the LTM entries with node numbers is next presented for the triangle T1 of Fig. 4. It begins by filling location (1, 1) in LTM with the index of the starting node (Node: 1). The location (2, 1) will be filled by the index of node at boundary B1 (Node: 13), while the node on B2 (Node: 16) will occupy the location (2, 2). Further, the location (3, 1) will be occupied by the index of node on B1 (Node: 14), while the location (3, 2) will store the interior node index (Node: 103) and the location (3, 3) in the matrix will be occupied by node on B2 (Node: 17). Naturally, the last row of the LTM will be occupied by nodes of B6 (Node: 2, Node: 28, Node: 29, Node: 30, Node: 3). Note that the nodes on boundary B1 occupy the locations in the first column of LTM, whereas all nodes on boundary B2 will occupy the diagonal locations of LTM and clearly the nodes on the third boundary B6 (in this example) will occupy the locations in the last row, i.e. (n + 1)th row of LTM. We are thus required to create 20 such LTMs to look up for the required information during the integration of governing equations. Since the nodes of icosahedral–hexagonal discretization belong to different categories, viz., Delaunay, Voronoi, boundary and interior nodes, the entire procedure of index storage could be described in a stepwise manner as follows: Step 1: Define Delaunay mesh: Join the two lower triangular matrices to compose a full [(n + 1) × (n + 1)] matrix represented by 10 rhombuses in Fig. 4. This constitutes the Delaunay mesh and it is stored in an array of size {10 × (n + 1) × (n + 1)}. 2010 Curtin University of Technology and John Wiley & Sons, Ltd. Figure 6. Search algorithm. Step 2: Define Voronoi mesh: From 20 LTMs which store information about grid nodes, the n 2 smaller equilateral triangles are identified in each one of them; it gives an array of size {(nt × 3)}. The centroids of these smaller triangles give the coordinates of points on Voronoi polygons. The 20 n 2 VPs are subsequently numbered from 1 to n 2 for the first main triangle; from (n 2 + 1) to 2n 2 for smaller triangles of the second main triangle and so on till the smaller triangles of the 20th main triangle have been numbered. The numbering of the smaller triangles has been shown in Fig. 6. Thus, the array that stores the Voronoi nodes of 20 triangles has the dimension {20n 2 }. Step 3: Assign Voronoi polygons to Delaunay nodes: Voronoi polygon around a Delaunay node is composed of smaller triangles which have this particular Delaunay node as their common vertex (Fig. 3). Now, the Voronoi polygons of pentagonal points can be identified from the (1, 1), (n + 1, 1) and (n + 1, n + 1) elements of each LTM. To identify the polygons for a node lying on any boundary of a triangle, one needs to restrict search to three boundaries [vertical (1, :), horizontal (n, :), slant (diagonal)] of two LTMs adjacent with respect to that boundary. Finally, for internal points of triangle, polygons are composed of smaller triangles lying within the corresponding LTM at their specified locations. This gives an array of size {np × 6}. Step 4: Assign Delaunay neighbours to Delaunay nodes: This information can be derived from the corresponding Voronoi polygons obtained in the previous step. That is, each 6-node (5-node for first 12 nodes) triangle whose centroids form a Voronoi polygon (which encloses the grid point in question) is considered and its vertices give the Delaunay neighbours of this particular node. Again this gives the array of size {np × 6}. Asia-Pac. J. Chem. Eng. 2011; 6: 110–119 DOI: 10.1002/apj Asia-Pacific Journal of Chemical Engineering ICOSAHEDRAL–HEXAGONAL GRIDS FOR CFD APPLICATIONS Step 5: Assign edge neighbours to Delaunay mesh and to Voronoi polygon edges: Delaunay neighbours to each Voronoi polygon edge and Voronoi neighbours to each Delaunay mesh edge are stored. The information like the rth edge of nth node is the sth edge of mth node is also stored. This information is very useful in designing discrete operators especially when Green’s theorem is applied. The array size in this case will be {nt × 6 × 2}. The first entry in the three-dimensional array stores the neighbouring node indices while second entry gives the local number of the common edges of the neighbouring polygons. Search procedure The search for a grid cell on the disctretization is one of the prime necessities during interpolation, adaptive grid refinements or while dealing with the transport of a scalar within a predefined computational domain. Therefore, efficiency of search algorithms warrants efficient computations at higher resolutions. As a matter of fact, grid-point search can be made highly efficacious, and the search algorithm used here is, in effect, of O(n) complexity at level-n resolution. This is a very significant step towards removing latency of computations in simulating transport of a scalar, which essentially results from the efficiency of the search algorithm presented here. This search algorithm, termed as the ‘nearest nodes navigator (NNN)’, exploits the parent–child node associations in finding the nearest grid-point (nearest node) to the target location (X ) in the following manner: Step 1: Identify which of the 20 main triangles of icosahedron contains the location X . Step 2: Identify in this icosahedron, the node nearest to X (say node A in Fig. 6). Step 3: Begin search from node A of this icosahedron triangle to identify its node nearest to point X : Set parent node = A and compute DPARENT = DIST (X , parent node) where DIST represents the spherical distance between X and parent node. Step 4: Identify the child nodes Child node1 and Child node2 for a parent node: For example, if parent node = A then from Fig. 6, Child node1 = B; Child node2 = C Step 5: Compute distances of X from the identified child nodes as DCHILD 1 = DIST (X , Child node1) DCHILD 2 = DIST (X , Child node2) Step 6: Next, find the new parent node using the following decision process: IF MIN(DCHILD 1, DCHILD 2) < DPARENT then DPARENT = MIN(DCHILD 1, DCHILD 2) and Parent node = Child node1, if DCHILD 1 < DCHILD 2 Parent node = Child node2, otherwise. Repeat from Step 4 ELSE Nearest node = Parent node; Stop The maximum number of floating point operations required in this process is 12 (Step 1) + 3 (Step 2) + 2 × (n + 1/2) (Steps 4–6), where n indicates resolution level-n. This node-search algorithm has been tested on regular icosahedral–hexagonal grid by Mittal et al .[17] for transport on a sphere free of unwanted numerical diffusion. However, this search algorithm is equally applicable on adaptive grid discretizations. Similar kind of approach has been earlier formulated by Giraldo[12,18] for solving advection and shallow water equations using Lagrange-Galerkin Methods. Local refinement and adaptive meshing Figure 7. Local refinement in icosahedral–hexagonal grid. This figure is available in colour online at www.apj ChemEng.com. 2010 Curtin University of Technology and John Wiley & Sons, Ltd. The quest for a unified atmospheric model[19] for simultaneous global- and regional-scale modelling is one of the preferred goals of designing future generation models.[7] In particular, increasing local resolution to better capture topography influences or important features of atmospheric circulation that result from the interplay of diabatic processes, adaptive grid techniques have become common and are rapidly evolving. Grid refinement on a sphere has been addressed in both grid-point and spectral techniques. Earlier Asia-Pac. J. Chem. Eng. 2011; 6: 110–119 DOI: 10.1002/apj 115 116 H. C. UPADHYAYA ET AL. Asia-Pacific Journal of Chemical Engineering (a) Base Grid (b) Resolved Grid Figure 8. Adaptive meshing on icosahedral–hexagonal grid. studies by Schmidt[20] using spectral technique, FoxRabinovitz et al .[21] using finite difference technique, Bacon et al .,[22] Giraldo,[12] Behrens et al .[23] and Jablonowski et al .[24] using finite element/finite volume method, and Fournier et al .[25] using a spectral element technique are a few of them. But adaptive gridding in the icosahedral–hexagonal framework has yet not been studied much, which indeed has great potential to downscale and upscale numerical models. This is a very important feature of this grid system. For this reason, we present here one such mesh refinement strategy that could be very attractive as shown in Fig. 7. The global resolution of the base grid in this figure is at level-8, on which further refinements with refine-4 and refine-8 are imposed on selected triangles. The strategy refine-k as described below signifies the local refinement on a Delaunay mesh of some base resolution. It is important to note that refine-k should be of the form 22k (k ≥ 0). This means that with k equals to 1, 2, and 3, the target triangle will be further subdivided into 4, 16, and 64 triangles, respectively. However, mesh refinements may result in hanging nodes which have to be dealt with during computations. Indeed this situation can be handled easily as any variable at the hanging nodes can be interpolated from its base grid values. This problem arises frequently in refined limited area modelling within coarser grid resolutions of a global model which provides the boundary conditions. But if adaptive mesh is chosen, then quite likely, numerical computations at the hanging nodes are also required, which may be performed in an efficient manner following a method described below. Let us consider the level-8 base grid resolution as shown in Fig. 8a and assume triangle T0 is to be refined further with refine-3. Then the systematic refinement of T0 and its neighbouring triangle steps (Fig. 8b) are as follows. Apply refine-3 (64 triangles) to T0 and refine-2 (16 triangles) to its Voronoi neighbours, i.e. triangles T1 , T2 , and T3 . The Voronoi neighbours T4 , T5 , T6 , T7 , T8 , and T9 of triangles T1 , T2 , and T3 will then be divided into four triangles each. That is, refinement 2010 Curtin University of Technology and John Wiley & Sons, Ltd. of these triangles is at refine-1 level. Once refine-1 level is reached then no further subdivision is required corresponding to target resolution of the central triangle (T0 in this case). It is important to notice that for all the hanging nodes generated in this process there always exists a node exactly opposite to it as shown in Fig. 8b. This can be made a part of the target grid by joining the central node with the opposite node (dashed lines in Fig. 8b) and one can always construct a pentagon (though not regular) around it. An important situation is encountered when two adjacent triangles like T0 and T3 have to have same refinement (say refine-k ); then application of this procedure to T0 will assign k − 1 level refinement on T3 , which is not desirable in view of earlier remark which suggests that the refinement level should be refine-k on T3 . In order to alleviate this problem following steps are required: (a) Start refinements of triangles in decreasing order, i.e. first the procedure is applied to triangles which are required to have highest refinement followed by triangles with the requirement of lower refinement. (b) While assigning refinement to Voronoi neighbours at any stage of the application of the abovementioned procedure, check which of the two, viz., existing refinement level or the refinement level, for which newer value has to be assigned, is greater. It is then required that the greater of these two values should be assigned to the Voronoi neighbour. Testing of computational framework: Tracer transport The conservation laws describe tracer transport of a scalar variable such as aerosols, hydrometeors and other chemical species in chemistry-climate models. On simplifying the equations of passive tracer advection (without source/sink), one obtains the following: ∂ρ + ∇.(ρ V ) = 0 ∂t ∂ρq + ∇.(ρq V ) = 0 ∂t (2) (3) Asia-Pac. J. Chem. Eng. 2011; 6: 110–119 DOI: 10.1002/apj Asia-Pacific Journal of Chemical Engineering ICOSAHEDRAL–HEXAGONAL GRIDS FOR CFD APPLICATIONS ∂ is the local time derivative; V is velocity field where ∂t and q represents the tracer concentration. Given the velocity field V (t) at time t, the distribution of tracer as time progresses from t to t + t can be computed with a stable numerical advection scheme in a uniform flow, warranting the initial shape of the tracer field being preserved in the absence of numerical diffusion. The integration of the flux forms of the continuity equation (2) and tracer conservation equation (3) converts the surface integrals to line integrals. Consequently, the discrete form of these equations may be written in terms of fluxes through the surface of a polyhedron with N faces. That is ∂m = Un ∂t n=1 (4) ∂qm = Fn ∂t n=1 (5) N N where m is the total air mass inside the polygon, q is the mean tracer mixing ratio in the same volume, and Un and Fn are the mass and tracer inward fluxes for face n, respectively. The tracer distribution within an element may be assumed constant, linear or higher order polynomial depending upon the required accuracy, time complexity and efficiency of the algorithm. The average tracer concentration crossing the interface boundary is then calculated in this procedure. This provides contributions from all boundary walls. The most common methods to estimate fluxes are upstream and centred schemes. But these methods have some typical drawbacks; upstream scheme is highly diffusive while centred scheme causes oscillations in the tracer field. Nevertheless, a typical second-order accurate scheme may be formulated in the following manner. Consider Fig. 5 where solid and dashed lines represent the Delaunay triangulation and the Voronoi polygons (or the control volumes), respectively. The control volume (CVo ) centred at O with Qi as its vertices has neighbouring control volumes with Pi as their centres. The mass and tracer fluxes at the boundaries of control volume CVo are computed as described by Miura[10] and the model is integrated using Matsuno-leapfrog time integration. The use of Matsuno step, after every five leapfrog steps during time integration, inhibits the growth of computational mode associated with pure leapfrog time integration.[26] A linear tracer distribution in each of the triangles OPi Pi +1 in surface coordinate (geographical) directions are computed and are then averaged over triangles having O as common vertex.[27] This gives the average tracer reconstruction for the cells centred at O. Finally, following Miura,[10] the incoming and outgoing fluxes from each edge of the control volume are 2010 Curtin University of Technology and John Wiley & Sons, Ltd. estimated as Fn = 6 ρi qi li (vRi n) and Un = i =1 6 ρi li (vRi n) i =1 In the above expressions, vRi is the resultant velocity at the midpoint (Ri ) of the corresponding edge and qi represents the average tracer mixing ratio for the outgoing flux. The method is mass-conserving secondorder accurate in time and space, but does not preserve monotonicity, and therefore one may use a slope or flux limiter to impose monotonicity. The consistency of the numerical method on the framework presented in previous sections has been examined for the most commonly used test case for any advection scheme.[28] In this test case, a structure of compact support is advected by a specified velocity field corresponding to solid-body rotation whose axis is not necessarily coincident with that of the rotation of the earth. The velocity components of the wind field are given by u = u0 (cos α cos ϕ + sin α cos λ sin ϕ) v = −u0 sin α sin λ The initial tracer distribution is prescribed as 1 (1 + cos(πr/R) if r < R q o (λ, ϕ) = 2 0 if r ≥ R (6) (7) where r represents the great circle distance from the centre of the cosine bell, (λ0 , ϕo ) = (3π/2, 0), and the bell radius R is 7π/64. The transport scheme should preserve the shape of the initial tracer distribution at Global error norms: solid-body rotation of cosine bell. This figure is available in colour online at www.apjChemEng.com. Figure 9. Asia-Pac. J. Chem. Eng. 2011; 6: 110–119 DOI: 10.1002/apj 117 118 H. C. UPADHYAYA ET AL. Asia-Pacific Journal of Chemical Engineering Numerical (dashed) and analytical (solid) cosine bell after one complete rotation. This figure is available in colour online at www.apjChemEng.com. Figure 10. all resolutions. Also, the scheme should not produce diffusion which may otherwise enlarge tracer distribution shape outward in the course of time integration. This phenomenon is tested rigorously in this numerical advection experiment. This numerical experiment has been performed for three resolutions: Level-16, Level-32, and Level-64 corresponding to 2562, 10 242, and 40 962 points on the sphere. The simulations are performed with Courant number ∼0.65 corresponding to each resolution in order to satisfy the CFL condition. The computed L1 , L2 , L∞ global error norms are shown in Fig. 9 showing the second-order accuracy of the scheme as found by Miura.[10] The analytical and simulated fields after one complete resolution are shown in Fig. 10, which represent that the numerical scheme described here preserves the tracer profile shape as desired, but since the limiters are not implemented here, the undershoots tracer profile develops. Reproduction of the Miura’s[10] results for this test case proves the consistency of the numerical framework presented here. CONCLUSIONS An efficient modelling framework has been developed for performing numerical computations on icosahedral–hexagonal grids capable of fast search and adaptive meshing. A novel method for indexing the nodes of the icosahedral–hexagonal grid has been presented in this work, which could be useful in CFD applications in chemical engineering, especially for computing 2010 Curtin University of Technology and John Wiley & Sons, Ltd. multiphase flows in the annular regions of rotating or stationary concentric spheres. One significant achievement of this development work is that an efficient search algorithm has been devised, which reduces the complexity of search from O(n 2 ) to O(n) on icosahedral–hexagonal grids where n refers to the refinement level. Therefore, the index search of nodes appears optimal because it can be easily done on the elements of lower triangular matrix of a [n + 1 × n + 1] matrix. Moreover, this framework has also been successfully tested for second-order accurate transport scheme, which is widely used in computational mathematics. Acknowledgements The authors would like to thank KDP Nigam for constant encouragement and the anonymous reviewers for their comments which improved the presentation of this paper. Rashmi Mittal and Hashmi Fatima would like to thank the Indian Institute of Technology Delhi for providing financial help in the form of a Ph.D fellowship to carry out this study. REFERENCES [1] M. Fernandino, C.A. Dorao, H.A. Jacobsen. Chem. Eng. Sci., 2007; 62, 6777–6783. [2] D.L. Williamson. J. Meteorol. Soc. Jpn., 2007; 85, 241–269. Asia-Pac. J. Chem. Eng. 2011; 6: 110–119 DOI: 10.1002/apj Asia-Pacific Journal of Chemical Engineering ICOSAHEDRAL–HEXAGONAL GRIDS FOR CFD APPLICATIONS [3] K. Sahr, D. White, A.J. Kimerling. Cartography Geographic Information Sci., 2003; 30, 121–134. [4] R. Sadourny, A. Arakawa, Y. Mintz. Mon. Wea. Rev., 1968; 96, 351–356. [5] R. Sadourny. Mon. Wea. Rev., 1972; 130, 1227–1245. [6] M.I. Boulos, D.C.T. Pei. Can. J. Chem. Eng., 1969; 47, 30–34. [7] T. Ringler, L. Ju, M. Gunzburger. Ocean Dynamics, 2008; 58, 475–498. [8] M. Satoh, T. Matsuno, H. Tomita, H. Miura, T. Nasuno, S. Iga. J. Comput. Phys., 2008; 227, 3486–3514. [9] D. Majewski, D. Liermann, P. Prohl, B. Ritter, M. Buchhold, T. Hanisch, G. Paul, W. Wergen, J. Baumgardner. Mon. Wea. Rev., 2002; 130, 319–338. [10] H. Miura. Mon. Wea. Rev., 2007; 135, 4038–4044. [11] J.R. Baumgardner, P.O. Frederickson. SIAM J. Numer. Anal., 1985; 22, 1107–1115. [12] F.X. Giraldo. J. Comput. Phys., 2000; 160, 336–368. [13] R. Heikes, D. Randall. Mon. Wea. Rev., 1995; 123: 1862–1880. [14] H. Tomita, M. Tsugawa, M. Satoh, K. Goto. J. Comput. Phys., 2001; 174, 579–613. [15] H. Tomita, M. Satoh, K. Goto. J. Comput. Phys., 2002; 183, 307–331. [16] Q. Du, V. Faber, M. Gunzburger. SIAM Rev., 1999; 41, 637–676. 2010 Curtin University of Technology and John Wiley & Sons, Ltd. [17] R. Mittal, H.C. Upadhyaya, O.P. Sharma. Mon. Wea. Rev., 2007; 135, 4214–4215. [18] F.X. Giraldo. J. Comput. Phys., 1997; 136, 197–213. [19] M.J.P. Cullen, T. Davies. Quart. J. Roy. Meteor. Soc., 1991; 117, 993–1002. [20] F. Schmidt. Beitr. Phys. Atmos., 1977; 50, 211–217. [21] M.S. Fox-Rabinovitz, L.L. Takacs, G.L. Stenchikov, M.J. Suarez, R.C. Govindaraju. Mon. Wea. Rev., 2001; 129, 453–469. [22] D.P. Bacon, N.N. Ahmed, Z. Boybeyi, T.J. Dunn, M.S. Hall, P.C.S. Lee, R.A. Sarma, M.D. Turner, K.T. Waight III, S.H. Young, J.W. Zack. Mon. Wea. Rev., 2000; 128, 2044–2076. [23] J. Behrens, N. Rakowsky, W. Hiller, D. Handorf, M. Lauter, J. Papke, K. Dethloff. Ocean Modelling, 2005; 10, 171–183. [24] C. Jablonowski, M. Herzog, J.E. Penner, R.C. Oehmke, Q.F. Stout, B. van Leer, K.G. Powell. Mon. Wea. Rev., 2006; 134, 3691–3713. [25] A. Fournier, M.A. Taylor, J.J. Tribbia. Mon. Wea. Rev., 2004; 132, 726–748. [26] F. Messinger, A. Arakawa. GARP Publication Series, Vol I, World Meteorological Organization: Geneva, 1976; p.64. [27] K.S. Yeh. J. Comput. Phys., 2007; 225, 1632–1652. [28] D.L. Williamson, J.B. Drake, J.J. Hack, R. Jakob, P.N. Swarztrauber. J. Comput. Phys., 1992; 102, 211–224. Asia-Pac. J. Chem. Eng. 2011; 6: 110–119 DOI: 10.1002/apj 119

1/--страниц