вход по аккаунту


Icosahedral-hexagonal grids on a sphere for CFD applications.

код для вставкиСкачать
Asia-Pac. J. Chem. Eng. 2011; 6: 110–119
Published online 9 July 2010 in Wiley Online Library
( DOI:10.1002/apj.479
Special Theme Research Article
Icosahedral-hexagonal grids on a sphere for CFD
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
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:
 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
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.
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
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
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
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
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
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
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
(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
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
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
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
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
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
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:
Parent node = Child node1, if DCHILD 1 <
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
 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
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
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
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
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
+ ∇.(ρq V ) = 0
Asia-Pac. J. Chem. Eng. 2011; 6: 110–119
DOI: 10.1002/apj
Asia-Pacific Journal of Chemical Engineering
∂ 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
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 =
ρi qi li (vRi n) and Un =
i =1
ρ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 + cos(πr/R) if r < R
q o (λ, ϕ) = 2
if r ≥ R
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
Figure 9.
Asia-Pac. J. Chem. Eng. 2011; 6: 110–119
DOI: 10.1002/apj
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
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.
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
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.
[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
[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,
[7] T. Ringler, L. Ju, M. Gunzburger. Ocean Dynamics, 2008; 58,
[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:
[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,
[16] Q. Du, V. Faber, M. Gunzburger. SIAM Rev., 1999; 41,
 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,
[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,
[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,
[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
Без категории
Размер файла
395 Кб
cfd, application, hexagonal, grid, icosahedral, sphere
Пожаловаться на содержимое документа