fully anisotropic unstructured grid generation with ... - CiteSeerX

European Conference on Computational Fluid Dynamics ECCOMAS CFD 2006 P. Wesseling, E. O˜ nate and J. P´eriaux (Eds) c TU Delft, The Netherlands, 2006

FULLY ANISOTROPIC UNSTRUCTURED GRID GENERATION WITH APPLICATION TO AIRCRAFT DESIGN V. Selmin∗ , E. Pelizzari† and A. Ghidoni†† ∗ Alenia

Aeronautica S.p.A. Corso Marche 41, 10146 Torino, Italy e-mail: [email protected] † Alenia Aeronautica S.p.A. Corso Marche 41, 10146 Torino, Italy e-mail: [email protected] †† University of Brescia, Mechanical Engineering Dep., Via Brianze 38, 25123 Brescia, Italy e-mail:[email protected]

Key words: surface meshing; anisotropic mesh; background grid; unstructured and hybrid grids Abstract. This paper describes a fully automatic 3D anisotropic mesh generation method for domains of arbitrary shape. The spacing of the boundary mesh is computed by the analysis of the principal curvatures and directions of the boundary surfaces. The spacing in the domain is obtained by interpolation of the spacing at the boundaries on a suitably constructed background mesh. Examples which illustrate the performance of the proposed methodology are presented.

1

INTRODUCTION

The availability of ever increasing computational resources and of efficient and robust CFD solution algorithms allows the numerical simulation of very complex flow fields. The generation of the grid is probably the bottleneck in terms of human resources required to complete a CFD simulation for domains of complex shape such as, for example, an aircraft configuration. The construction of a grid which minimizes the number of nodes required to compute a flow field within a prescribed tolerance level, or even of a grid which approximates the geometry of the domain with sufficient accuracy to allow the computation of an initial coarse solution in an adaptive solution process, remains in fact a difficult and time consuming task with the commonly available software. Among the many tasks which require substantial user input in grid generation, that of prescribing the mesh spacing for a complex domain is probably the most demanding in terms of human resources required. We here describe a three-dimensional anisotropic 1

V. Selmin, E. Pelizzari and A. Ghidoni

unstructured grid generation algorithm which instead requires a very limited user input even for very complex geometries since the mesh spacing is automatically defined as a function of the curvature of the boundary of the domain. We follow the widely adopted approach [3, 5] which consists in augmenting the domain to be meshed with a metric field M(x) (i.e. to define a Riemannian structure on the domain) and in requiring that all the mesh edges have the same length with respect to M(x) (typically unit length edges lengths in M(x) are chosen). This leads to an anisotropic control of the mesh spacing which means to control not only the mesh element size but also their aspect ratio. This approach decreases significantly the nodes number used when the flow field or the geometry of the domain are characterized by anisotropic features Following the ideas introduced in [4], the metric field on the boundary is computed according to the principal curvatures and directions of the boundary surface. This approach is, however, not appropriate if the surface curvature is very small, since it leads to excessive mesh spacings. The generally adopted approach to deal with “nearly-flat” surfaces is to limit the spacing with a threshold value prescribed in an isotropic manner. As a consequence all the informations regarding the anisotropy of the mesh are completely lost in regions of small surface curvature and much of the potential savings allowed by anisotropic unstructured meshes are not fully realized in practice. This situation is typical in many applications such as, for example, high aspect ratio wings and turbomachinery bladings. An original anisotropic treatment for “nearly-flat” surfaces has been developed in this work. After establishing the metric field distribution on the boundary, we proceed with the surface mesh generation using an approach, sometimes referred to as the “indirect approach”, which consists in generating the mesh in the parametric domain associated to the surface and in mapping the resulting mesh onto the real surface. Using the surface mesh as input data it is then possible to build a background grid, which is used to interpolate the value of the boundary metric field at any other points. The computational domain is then discretized by using a 3D anisotropic mesh generator based on Front Advancing method. In order to build hybrid grids suited to solve Navier-Stokes equations, layers of prisms are introduced at the skin boundary, by using a spring analogy based node movement algorithm. Starting with layers of null thickness, a force proportional to the required displacement is applied to the skin nodes which result to push on the grid formed by tetrahedra, leaving then room to introduce prisms. The anisotropic grid generation approach has been applied with success to the computation of inviscid and viscous flows around complex aeronautical configuration configurations. In addition to provide unstructured/hybrid grids of good quality formed by a relatively small number of elements, it allows to control some desirable properties of the grid. As an example, it is possible by using this method to maintain the chordwise number of nodes along the wing span almost constant, which is hard to obtain by using an isotropic unstructured mesh generation process.

2


2

RIEMANNIAN METRIC

An anisotropic control of the mesh spacing has been chosen for this work, which means that we will be able to control not only the mesh elements size but also their aspect ratio. This is realized by requiring that all the edges of the mesh are of unit length with respect to a suitable defined metric field M(x), where x is a point of Rd (d = 2, 3). M(x) is a (d × d) symmetric positive definite matrix. For example, in two dimensions, we consider a(x) b(x) , (1) b(x) c(x) such that a > 0, c > 0 and ac − b2 > 0, for a, b, c ∈ R. If the field of tensors thus defined is known, it induces a Riemannian structure over Rd . The length with respect to M of an arbitrary mesh edge v is approximated as p ||v||M = v T M(b x)v , (2)

b is the midpoint of edge v. When M(x) is independent of the position x, we where x again find the classical Euclidean case, otherwise we are in the Riemannian context. The matrix M(x) can be factorized as M(x) = R(x)Λ(x)R(x)T ,

(3)

R(x) = [r 1 (x), . . . , rn (x)]

(4)

where is the orthogonal matrix of the normalized (right) eigenvectors of M(x) and Λ(x) = diag(γ1 (x), . . . , γn (x)) is the diagonal matrix of real eigenvalues γi > 0, i = 1, . . . , n. We can associate to M(x) a transformation matrix T(x) defined as p p T(x) = M(x) = R(x) Λ(x)R(x)T ,

(5)

(6)

and we therefore can express M(x) in terms of T(x) as

M(x) = T(x)T T(x) = T(x)T(x)T .

(7)

The transformation matrix T maps an edge v of unit length with respect to M into an edge w = Tv of unit length in the usual Euclidean norm. We have in fact that √ √ √ ||w|| = wT w = v T TT Tv = v T Mv = ||v||M .

3


3

GEOMETRICAL CHARACTERISTICS OF AN ELEMENT

The geometrical characteristics of an element can be defined in terms of a set of n orthogonal directions αi and n associated element sizes δi , where n is the number of dimensions. The transformation matrix T can then be regarded as the result of combining n scaling operations, with factor 1/δi , in each direction αi , i.e. n X 1 T= αi αTi . δ i i=1

(8)

By rewriting Eq. (6) as T=

n X √

γi r i r Ti

(9)

i=1

and by comparing with Eq. (8) we therefore have that 1 δi = √ , γi

αi = r i .

(10)

The spacing s(β) in the generic direction defined by the unit vector β is given by s(β) = ||Tβ||−1 . 4

MESH SPACING CONTROL ON THE BOUNDARY

The objective of surface mesh generation is to construct a triangulation in which the maximum distance from the real surface is below a prescribed threshold value — this is a useful and practical way to arrive at meshes which approximate the real surface as accurately as required. In practice, we have however to construct a metric field M 3 on the real surface Σ so that the above criterion is satisfied by unitary (in M3 ) mesh edges, and a corresponding metric M2 in the parametric space Ω which will be used as input data by the surface mesh generator. The manner in which M3 and M2 are constructed is described in the following. 4.1

Mesh spacing based on local curvature

Let Σ be a surface defined by the parametrization r(u, v) ∈ R3 , P a point of Σ of coordinates r P (uP , vP ), n the normal to Σ at P , and t the unit tangent vector along some curves on Σ passing at P . By locally approximating in the neighborhood of P the curve by its osculating circle, the ratio ε between δ, the maximal gap between the osculating circle and its chord of arc, and the osculating circle radius ρ, is given by q δ ε ≡ = 1 − 1 − (λ/2)2 , (11) ρ 4


where λ = ht /ρ, and ht = λρ is the chord length. The parameter λ is therefore given as a function of ε according to p ht = 2 ε(2 − ε). λ≡ (12) ρ The above equation is employed to compute the chord length ht for a prescribed value of ε — a parameter which measures how closely the chord approximates the real curve. The geometrical characteristics of a surface element may be obtained by applying the above analysis in the surface principal directions, see e.g. [1], i.e. the directions in which the surface curvature assumes its minimum and maximum values. The two principle directions, di ∈ R3 , can be defined as J(r)v˙ i di = q v˙ Ti Gv˙ i

i = 1, 2 ,

where J(r) is the Jacobian matrix of the mapping r = r(u, v), i.e. ∂r ∂r J(r) = , , ∂u ∂v

(13)

(14)

G = JT J is the first fundamental matrix of the surface Σ, and v i ∈ R2 the principal directions in the parametric space u − v. Notice that, by construction, the direction d1 and d2 are orthogonal. The element sizes hi associated to the directions di are, then, computed as p p ε(2 − ε) , (15) hi = 2ρi ε(2 − ε) = κi where ε is a user prescribed parameter, ρi the principal radii of curvature and κi the principal curvature. Finally, the metric tensor M3 (P ) associated to the surface is computed according to 3 X 1 M3 (P ) = d dT , 2 i i h i=1 i

(16)

where d3 = n and h3 = min(h1 , h2 ). The metric M2 (P ) can be simply expressed as M2 (P ) = J(r)T M3 (P )J(r) ,

(17)

where J is defined in Eq. (14). 4.2

Treatment of “nearly-flat” surfaces

The values of hi have to be occasionally limited, in order to avoid too large or too small elements. If h1 and h2 denote the smallest and the largest spacings, respectively, 5


˜ i are defined as the limited spacings h ˜ 1 = min(max(h1 , h1,m ), h1,M ) h (18) ˜ 2 = min(max(h2 , h2,m ), h2,M , sM h1 ) , h where hi,m and hi,M are the minimum and maximum value that hi can assume. The parameter sM represents the maximum value of element stretching. By choosing sM = 1, an isotropic element will be formed. The quantities hi,m and hi,M are defined as follows. Let denote by lu and by lv , the length of the curves represented parametrically by equations u = [uP , v(s)]T and u = [u(s), vP ]T , respectively. The spacing hu and hv are computed as a fractional part of the lengths lu and lv , i.e. lu lv hu = , hv = , (19) nu nv where the parameters nu and nv (nu , nv ≥ 1) are constants (integer) defined on Σ. The quantities hi,M are computed by using the Euler theorem  1  2 X 2 1 ϕi,u 2  hu  = cos , 1 ϕi,v h2i,M 2 i=1 hv

(20)

where ϕi,u and ϕi,v are the angles formed by the principal direction di with the tangent vectors to the curves u = [uP , v(t)]T and u = [u(t), vP ]T , denoted by tu and tv , respectively. In addition, hi,m are defined as a fraction of hi,M or by taking a constant value on the surface. When the geometric model is too complex, the computation of hi,M is performed by using the concept of background surface spacing characteristics, which are built by the knowledge of the characteristic dimensions of the object to be discretized and of the tangent vectors to the surface. When the surface is a plane (ρi = ∞) or when the principal curvatures are equal, the principal directions of normal curvature are undefined. However the directions t u and tv ˜i still exist and can be chosen as directions di , if tu · tv = 0. In this case the spacings h associated to the directions di are, respectively, hu and hv . if tu · tv 6= 0, d1 is chosen coincident with the direction which specifies the minimum spacing and d2 = d1 × n. In ˜ 1 = min(hu , hv ) and h ˜ 2 is computed by using the equation (20). this case h 4.3

Metric intersection

A prerequisite for the grid generation algorithm is the specification of only one metric at each point. However, the metrics are multi-defined on the boundary curves and at each corner. Thus, in these points, a metric intersection is requested. 6


Let P be a point in the space at which two different metrics exist, say Ma (P ) and Mb (P ), defined as 3 X 1 T Ma (P ) = 2 αi αi δ i=1 i

3 X 1 T Mb (P ) = 2 βi βi . γ i=1 i

(21)

The directions ξ i of the intersection metric are chosen coincident with the directions associated to the metric matrix which specifies the minimum spacing, i.e. ξ i = αi . The new metric Mnew (P ) is defined as 3 X 1 Mnew (P ) = ξ ξT , 2 i i h i=1 i

where

4.4

1 hi = min δi , kTb ξ i k−1

(22)

.

(23)

Metric smoothing

Due to the fact that the metrics on the boundary curves and at each corner are, in general, multi-defined, discontinuity or large gradient may appear in the metric field, which can lead to the generation of bad quality elements. A spring analogy like algorithm is used to propagate and regularize the metric field. It is first assumed that each node i is connected to each adjacent node j by a fictitious spring under the generalized force Fij defined by −1 Fij = Kij (T−1 j − Ti ) ,

(24)

where the scalar Kij is the spring constant, Ti and Tj are the transformation matrices associated to nodes i and j. The spring constant can be expressed as Kij = Kijsm Kijlg ,

(25)

where

ij Aij e1 + Ae2 Ai is the spring constant associated to a Laplacian smoothing algorithm and s 1 (xj − xi )T Mij (xj − xi ) lg Kij = α (xj − xi )T (xj − xi )

Kijsm =

(26)

(27)

ij is related to the locally required grid spacing. Aij e1 and Ae2 are the areas of the two elements adjacent to side ij, Ai is the area of the elements surrounding node i, α a scaling

7


factor, which is taken equal to the maximum value of Kijlg (with α = 1) over the patch of elements surrounding node i, Mij the metric matrix along side ij. The resulting metric field is the solution of the equilibrium system for each vertex i: X Fij = 0 , (28) j∈Ki

where Ki is the set of nodes surrounding i. In this work, a few iterations of the following pointwise scheme is actually used to smooth the metric field: X −1 T−1 = T + Fij,(p−1) i,(p) i,(p−1) j∈Ki

=

T−1 i,(p−1)

+

X

j∈Ki

5

−1 Kij (T−1 j,(p−1) − Ti,(p−1) ) .

(29)

MESH SPACING CONTROL IN THE VOLUME

The mesh spacing in the volume is controlled through a background grid, which is used to interpolate the boundary value of the metric field at any other point of the domain to be discretized. The procedure in order to build the background grid is described in the following. The process starts constructing the Delaunay triangulation [7, 8] of the points which discretize the skin surface and the external boundaries. The circumcenters of the tetrahedra that have at least one node on the skin and that are not inside the body are computed. ¯ il is obtained by solving the minimization problem An “average” direction v min

nj X j=1

v ij ¯ il · 1−v kv ij k

2

,

(30)

where v ij are the vectors that connect the node i on the skin with the set of circumcenters j associated to the elements which have node i as vertex, and nj is the number of circumcenters belonging to the set. Then, the vector v il is constructed connecting the ¯ il with node i with the local node l which is located at the intersection of the direction v the surface formed by the nj circumcenters. The procedure is repeated by taking the new layer of nodes l as external boundary until a maximum of iterations is reached. This results in the construction of a succession of nodes layers which become closer and closer to the skin surface. The background grid is defined by building the Delaunay triangulation of the nodes on the skin, on the external boundary and on additional layers. The transformation matrix Tl at the node l is computed as −1 −1 Tl = ((1 − ω)T−1 i + ωT∞ )

8

ω ∈ [0, 1] ,

(31)


where T∞ is the transformation matrix based on the element characteristics at the external boundary, Ti the transformation matrix associated to the skin of nodes i and ω a parameter defined as 2 kv il k kv il k +b ω = f (kv il k, r∞ , b) = (1 − b) , (32) r∞ r∞ where r∞ is the distance from the skin to the external boundary and the coefficient b the slope of the function at its origin. 6

HYBRID MESH GENERATION

In the case of the simulation of high Reynolds number viscous flows, it is appropriate, knowing the physics of boundary layers, to consider a form of a priori adaptation to reflect the difference in gradients in flowfield variables across a boundary layer as compared with those along the boundary layer in the direction of the flow. If such an approach is followed then elements with high aspect ratios are required. This leads to build quadrilaterals in 2D and prisma in 3D within some regions of the computational domain, in particular close to solid boundaries, by using the method of advancing layers (or advancing normals). The basic approach amounts to growing layers of elements by advancing along lines which are approximatively normal to the boundary. These layers of elements are grown until a specified distance from the boundary is reached. After this, the domain is filled with triangles in 2D and tetrahedra in 3D by using standard approaches. 6.1

Basic algorithm

Let denote by xki the position of node i on the layer k, its new position on layer k + 1 may be computed according to the formula = xki + dki nki xk+1 i

(33)

where nki and dki represent the normal vector and the displacement at node i with respect to layer k, respectively. The vector nki is obtained from the knowledge of the normal vector nke related to the set Eik of elements belonging to the patch Pik of elements surrounding node i by solving the minimisation problem 2 ne X nke k min (34) 1 − ni · k nke k e=1

where ne is the number of elements that belong to the set. In general, the normal vector is subsequently smoothed by using a few iterations of a Laplacian smoothing procedure X k,l k,l k,l k K n − n (35) nk,l+1 = n + 0.5 ω ij i j i i i j∈Kki

9


where Kik is instead the set formed by the nodes on the boundary of Pik and ωik is a parameter that belongs to the interval [0,1]. The parameter ωik may be taken as a constant value or be expressed as a function of the layer indicator k, e.g. k−1 k ,1 , 1 < k l ≤ kt (36) ωi = min kl − 1 where kt and kl are the total number of layers and the indicator of the layer for which ωik reaches unity, respectively. Note that for k equal to one, which corresponds to the initial surface, the normal vector is not altered if formula (36) is used. The smoothing procedure has not to be applied to nodes belonging to curves or to corners for which the geometry presents concavities. The displacement dki may be assumed constant or variable from nodes to nodes depending on the applications. Let denote by d1i the distance of the first layer from the surface, and having establish the distance dmax and the value of the layer indicator kt at which the quasi-structured layer may finish, the distance between layers k + 1 and k is given by dk+1 = κi dki (37) i where the parameter κi represents the geometrical ratio and is solution of dmax = d1i

kt X

(κi )j−1

(38)

j=1

In the case of nodes belonging to curves or to corners for which the geometry is concave, the value of the displacement dki needs to be corrected. A node on a curve belongs by definition to two surfaces and two normal vectors can be associated to it, one for each surface. The normal vector associated to the node is finally computed as the direction that bissects the angle formed by the normal vectors related to the two surfaces. Note that this direction represents also the solution obtained by using an equation similar to equation (34) but based on normal vectors associated to surfaces instead of elements. At the nodes of the curves where the geometry is not smooth the scalar product of those vectors is small. If this does not represent any problem if the geometry is convex, the value of dki needs to be corrected in the case of concavities according to the formula d¯ki = dki /sp

(39)

where sp represents the scalar product between the normal vector at the node and the normal vectors associated to the surfaces. Similar techniques are applied to corners. In particular, the displacement at a corner is computed as the maximum of the displacement values at this node associated to the curves that intersect at the corner. The basic algorithm for building quasi-structured layers may be linked up with the unstructured grid generation approach in two ways: 10


1. After the generation of the quasi-structured grid, the unstructured grid generation process is performed starting from the discretisation of the outer layer. 2. Starting from an already existing unstructured grid, the nodes on the surface are moved in order to recover the position of these nodes on the outer layer of the quasi-structured grid. A node movement algorithm is then used that transfers the new metric of the unstructured grid surface in the 3D grid from the surface to the external boundary. An alternative is to apply this procedure each time a layer is built. 6.2

Weak formulation

Let xi indicates the position vector of node i and Ei the set of elements belonging to the patch Pi of elements surrounding a given node i. Ki is instead the set formed by the nodes on the boundary of Pi . A node movement algorithm, based on the spring analogy, assumed that each node i is connected to each adjacent node j by a fictitious spring under the force Fij defined by Fij = Kij (xj − xi ) (40) where Kij is the spring constant which in general will depend on some local grid features. The resulting mesh is for each node i the solution of the equilibrium system: X Fij = 0 (41) j∈Ki

which may be rewritten as X

j∈Ki

Kij (∆xj − ∆xi ) = −

X

j∈Ki

Kij (x0j − x0i )

(42)

where ∆xi = xi − x0i represents the deviation of the actual position xi from its initial value x0i . In order to modify the grid in such a way to maintain basic characteristics of the original one, additional terms may be introduced according to X X Kij (∆xj − ∆xi ) = Kij (Qij − I) (x0j − x0i ) (43) j∈Ki

j∈Ki

where Q is a transformation matrix which is related to the change between the original and deformed surface grids. Let consider the quasi-structured grid built between two layers whose layer indicators ¯ i the set of are k and k + 1, respectively. For a node i of the layer k + 1, let denote by K ˜ node-pairs ij that belong to the layer k + 1 and by Ki , the set of node-pairs which connect ˜ i being formed by one element. the node i to the corresponding node on layer k; the set K S˜ ¯ Note that Ki ≡ Ki Ki . 11


A revised version of equation (43) may be built according to X X X X Kij (xj − xi ) + Kij (xj − xi ) = Kij Qij (x0j − x0i ) + Kij Sij ¯i j∈K

¯i j∈K

˜i j∈K

(44)

˜i j∈K

which may be rewritten as X X ˜ i] = ˜ i ] (45) Kij (xj − xi ) − Kij (xi − xj ) [j ∈ K Kij Qij (x0j − x0i ) − Kij dj nj [j ∈ K ¯i j∈K

¯i j∈K

where we assumed Sij = − dj nj . If the condition ˜ i ] Kij [j ∈ K ¯ i] Kij [j ∈ K

(46)

is satisfied, equation (45) reduces to xi = x j + d j n j

(47)

This represents the basic idea behind the use of a node movement algorithm for hybrid grid generation. Starting from an initial unstructured grid, it is then proposed to build the hybrid grid layers according to the following procedure 1. Add a layer of zero thickness to the discretised surface grid through the duplication of its nodes and build the associated quasi-structured grid, that is formed by prismatic elements, 2. Compute the metric matrix on the surface grid, taking also into account as the normal distance the prescribed displacement between layers, 3. Compute the stiffness constants Kij on the unstructured and quasi-structured grids, 4. Form the right-hand side of equation (45) by assuming Qij =I; the quantities dj and nj being computed by using the approach described in the previous section, ¯ i represents now the 5. Perform node movement. Referring to equation (45), the set K set Ki related to the unstructured grid, 6. Save the coordinates of the nodes of the new layer and update the interface with the unstructured grid, 7. Repeat steps 1 to 6 starting from the update unstructured grid until the prescribed number of layers is reached.

12


We recommend to normalise the node-pairs stiffness in such a way that they satisfy X ˜ i] = Kij [j ∈ K Kij = 1 (48) ¯i j∈K

that means to give the same weight to the two grids in the mesh deformation process. Due to the interaction between the hybrid and unstructured grids, the weak approach allows to obtain a better control of the overall grid, the unstructured grid offering resistance to the development of the hybrid one, in particular in concave areas. In addition, the use of the node movement technique improves grid smoothness. Condition (48) does not allow in general to get the prescribed spacings within the prismatic layers, a lower value (> 0.5 the prescribed value) being obtained. In order to recover this requirement in a given number of layers close to the body surface, the value of the node-pairs stiffness ¯ i ] may be multiplied by a small number, when constructing those layers. This Kij [j ∈ K indirectly improves orthogonality in the prismatic layers close to the body. 7

EXAMPLE OF APPLICATIONS

Results of the turbulent flow around a wing-body configuration are here illustrated and compared with experimental data. The turbulence Baldwin-Barth one-equation model has been adopted. The grid is formed by around 780000 nodes (13500 on the body) and 2500000 elements, i.e. 1420000 tetrahedra and 1080000 prisms in the boundary layer quasi-structured grid, respectively. This relatively small number of nodes for such computation is due to the anisotropic nature of the grid elements that allows to reduce by an order of magnitude the number of nodes with respect to an isotropic grid based on equivalent mesh spacings. The flow conditions are M∞ = 0.85, α = 2 and Re = 54.2×106 . The “boundary layer” grid is extended from the skin to a distance up to 2% of the value of the aerodynamic mean chord and is formed by 40 layers. Around 150 nodes discretise the wing sections, this number of nodes being approximatively maintained along the overall wing span. Figures 1 to 3 mainly illustrate the grid characteristics while figure 4 presents a comparison between numerical and experimental solutions. The numerical solution is in close agreement with experimental data. This simulation represents at our knowledge the first application of a Navier-Stokes solver on a fully anisotropic unstructured/hybrid grid. 8

CONCLUSIONS

Most of the avaible unstructured grid generators can only construct isotropic meshes. This is a serious limitation since the use of anisotropic meshes can result in significant savings when the flow field or the geometry of the domain are characterized by anisotropic features. This is indeed the case in many applications and is typical in CFD. In this paper an original procedure to anisotropically prescribe the metric field on both the surface and the volume of the computational domain is presented. The spacing on the 13


boundary mesh is computed by the analysis of the principal curvatures and directions of the surface, and the spacing in the domain is obtained by propagating the spacing at the boundaries by means of a background grid. The special treatments required to deal with nearly flat surface and to regularize the spacings are thoroughly discussed in the paper. A reduction of nearly one order of magnitude of nodes with respect to the isotropic case can be obtained by the proposed methodology. In order to build hybrid grids suited to solve Navier-Stokes equations, a spring analogy based node movement method has been here introduced that adds layers of prisms at the skin boundary while maintaining good smoothness characteristics of the overall grid. REFERENCES [1] G. Farin, Curves and Surfaces for Computational Aided Geometric Design, A Practical Guide, Academic Press, 1997. [2] P. Frey, P.L. George, Mesh Generation, Application to finite elements, Hermes, 2000. [3] J. Peraire, K. Morgan, L. Formaggia, J. Peiro, O.C. Zienkiewicz, Finite element Euler computations in three dimensions, International Journal for Numerical Methods in Engineering 1988; 26:765–781. [4] H. Borouchaki, P. Laug, P.L. George, Parametric surface meshing using a combined advancing-front Delaunay approach, International Journal for Numerical Methods in Engineering 2000; 49:233–259. [5] H. Borouchaki, P. Laug, F. Hecht, E. Saltel, P.L. George, Delaunay mesh generation governed by metric specification. Part I: algorithm, Finite Elements in Analysis and Design 1997; 25:61–83. [6] T.J. Baker, Automatic Mesh Generation for Complex Three-Dimensional Regions Using a Constrained Delaunay Triangulation, Engineering with Computers 1989; 5:161–175. [7] A. Bowyer, Computing Dirichlet Tessellation, The Computer Journal 1981; 24:162– 166. [8] D.F. Watson, Computing the n-dimensional Delaunay tessellation with application to voronoi polytopes, The Computer Journal 1981; 24:167–172. [9] A. Ghidoni, E. Pelizzari, S. Rebye and V. Selmin, 3D Anisotropic Unstructured Grid Generation, Int. J. Numer. Methods Fluids, 2006.

14


Surface grid overview

Pressure coefficient contours Figure 1: Skin grid overview and Cp contours

15


Wing root region

Enlargement at leading edge

Wing crank region


Wing tip region


16 Figure 2: Skin grid characteristics


Skin boundary discretisation

Outer viscous layer discretisation

Figure 3: Skin and outer viscous layer discretisations

17


Cp

Cp

numer exper

numer exper

eta=0.216

0

0.2

0.4

eta=0.470

0.6

0.8

1

0

0.2

0.4

x/c

Cp

Cp

numer exper

0.2

0.4

0.6

0.8

1

0

0.2

0.4

Cp

numer exper

0.8

1

numer exper

eta=0.837

0.4

0.6 x/c

Cp

0.2

1

eta=0.750

x/c

0

0.8

numer exper

eta=0.610

0

0.6 x/c

eta=0.950

0.6

0.8

1

0

x/c

0.2

0.4

0.6 x/c

18 Figure 4: Comparison between numerical and experimental data

0.8

1