Citation |

- Permanent Link:
- https://ufdc.ufl.edu/UF00100807/00001
## Material Information- Title:
- Three-dimensional representation of evolving interfaces in computational fluid dynamics
- Creator:
- Khanna, Abhishek, 1976- (
*Dissertant*) Shyy, Wei (*Thesis advisor*) - Place of Publication:
- Florida
- Publisher:
- State University System of Florida
- Publication Date:
- 2000
- Copyright Date:
- 2000
- Language:
- English
## Subjects- Subjects / Keywords:
- Coordinate systems ( jstor )
Curvature ( jstor ) Curves ( jstor ) Grid refinement ( jstor ) Polynomials ( jstor ) Radius of a sphere ( jstor ) Radius of curvature ( jstor ) Surface contours ( jstor ) Triangles ( jstor ) Triangulation ( jstor ) Aerospace Engineering, Mechanics, and Engineering Science thesis, M.S ( lcsh ) Dissertations, Academic -- Aerospace Engineering, Mechanics, and Engineering Science -- UF ( lcsh ) - Genre:
- bibliography ( marcgt )
theses ( marcgt ) government publication (state, provincial, terriorial, dependent) ( marcgt ) non-fiction ( marcgt )
## Notes- Abstract:
- Moving boundary problems arise in numerous important physical phenomena and they often form complex shapes during their evolution. Numerical methods for the solution of moving boundary problems primarily focus on interface tracking and on the solution of the field equations. The ability to track interfaces in two dimensions is well established. However, modification of the grid, representing the interface, as it evolves in a three-dimensional space introduces additional issues. In the current work, moving interfaces are represented by triangulated surface grids. Continuous spline surfaces are also used as an alternative means of surface representation. The unstructured triangulated grids are restructured and refined based on the shape and size of the triangular elements in the grid. As the interface deforms, points are automatically added and mesh connectivity modified to ensure that the accuracy of the interface representation remains consistent. Normals at the nodes of the surface mesh are required for the movement of the boundary and curvature plays an important role in problems involving capillarity. In the current work, methods for computing geometric features such as normals and curvatures on a triangulated surface are discussed. Curvature calculations are based on the mean of normal curvatures and on approximations to the Dupin indicatrix. Normals and curvatures are also computed on the parametric spline surface. This is done by the calculation of the partial derivatives of the coordinates of the surface nodes with respect to the parameter values. Results are presented to show the complex feature representation capability of the tracker code. A comparison between the different methods, used for computing curvatures, is made here. Finally, the effect of mesh quality and accuracy of the point location on the curvature error is also investigated. ( , )
- Thesis:
- Thesis (M.S.)--University of Florida, 2000.
- Bibliography:
- Includes bibliographical references (p. 89-90).
- System Details:
- System requirements: World Wide Web browser and PDF reader.
- System Details:
- Mode of access: World Wide Web.
- General Note:
- Title from first page of PDF file.
- General Note:
- Document formatted into pages; contains viii, 91 p.; also contains graphics.
- General Note:
- Vita.
- Statement of Responsibility:
- by Abhishek Khanna.
## Record Information- Source Institution:
- University of Florida
- Holding Location:
- University of Florida
- Rights Management:
- All applicable rights reserved by the source institution and holding location.
- Resource Identifier:
- 47681989 ( OCLC )
002678729 ( AlephBibNum ) ANE5956 ( NOTIS )
## UFDC Membership |

Downloads |

## This item has the following downloads: |

Full Text |

THREE-DIMENSIONAL REPEII.l.ETATION OF EVOLVING INTERFACES IN COMPUTATIONAL FLUID DYNAMI( ' By ABIIISIIEK KHANNA A THESIS PRESENTED TO THE GRADUATE SCHOOL OF THE UNIVERSITY OF FLORIDA IN PARTIAL FULFILLMENT OF THE REQUIREMENTS FOR THE DECREE OF MASTER OF SCIENCE UNIVERSITY OF FLORIDA 2000 ACKNOWLED( ;I I I NTS I would sincerely like to thank Dr. Wei Shyy for his useful si'-.f.tii,, belief in my capabilities and for extending my financial support to continue this work. I express my gratitude to Dr. H.S. Udaykumar for his constant support and guidance throughout this work. I would also like to thank him for spending a considerable amount of his valuable time with me during the course of this research work. I wish to thank Dr. Jorg for his help and for the useful discussions we have had together. My thanks are due to Vivek .1, ii ., iiI' for helping me in getting familiar with the tracker code developed by him. I would also like to thank Ramji Kamakoti and Pr.ii:m-.r Vaidyanathan for helping me in making my report. Equally, I wish to express my heartfelt gratitude to my parents for their endless love and encouragement. TABLE OF CONTENTS page ACKNOWLED( ;' I'I'NTS ........................................ ii NOMENCLATURE ............. ................................ v ABSTRACT ........................................ ............ vii CHAPTERS 1 INTRODUCTION .......................................... 1 1,ving-boundary Problems .......................................... 1 Three-dimensional Interface Tracking ..................................... 3 Overview of the Report ........... ............................. 5 2 REVIEW OF CUIRVEI'iS AND SURFACES ............................... 6 Implicit and Parametric Forms .. ........................ ...... ....... 6 Bezier Curves ............................................ 8 Polynomial Bezier Curves .......................... ............. 8 Rational Bezier Curves ........................................ 12 Bezier Surfaces ........................................... 13 Polynomial Tensor Product Bezier Surfaces ............................... 14 Rational Bezier Surfaces ........................... ............ 15 B-Spline Curves ............................................ 16 Polynomial B-Spline Curves ............... .......................... 17 Nonuniform Rational B-Spline Curves ................ ................ 21 B-Spline Surfaces .......................................... 22 Polynomial B-Spline Surfaces ...................... ................. 22 Nonuniform Rational B-Spline Surfaces .................................. 23 Bezier Triangles ............................................ 24 Polynomial Bezier Triangles ................ .......................... 24 Rational Bezier Triangles ........................................ 28 3 INTERFACE TRACKING METHOD ............ ................ 30 C1 Surface Spline Algorithm .. ....................................... 34 1!.i Refinement ........................................... 40 Point Insertion ........................................... 42 Diagonal Swapping ............ ............................. 46 .1,-li Update .............. ................................... .......... 48 l1.. li Update on C1 Spline Surface ............. ...................... 48 l1 i 'i Update on Surface Triangulation ................................. 50 4 RESULTS ............................................. 57 Complex Feature Representation ......................................... 57 Grid Refinement ........................................... 62 No Grid Refinement ......................................... 68 Grid Adaptation ........................................... 70 5 CONCLUSIONS AND FUTURE WORK ................................ 81 APPENDICES A POINTS AND VECTORS ............................................ 83 B BARYCENTRIC COORDINATES ............ ......................... F C DATA STRUCTURES FOR INTERFACE TRACKING ................ 86 Node Data Structure .................................................... 86 Edge Data Structure ......................................... 86 Triangle Data Structure .................................. ............. 87 Surface li \,,Ii Data Structure ............................................ 87 REFERENCES ...................................... ........... 89 BIOGRAPHICAL SKETCH ................................................. 91 NOMENCLATURE The following list gives the meanings of symbols that appear in this report. A, B, C Coefficients of Dupin indicatrix Ai Edge coefficient a, b Real numbers ail, ai2 Blend ratios of cell i Bi,n Bivariate Bernstein polynomial of degree n Bi,, Univariate Bernstein polynomial of degree n C Curve C' Preliminary cell center C Cell center c Cell d Direction vector d, e, f Elements of direction vector d d Distance of projection of new point e Euclidean or point space E, F, G Coefficients in the mean curvature equation el, e2, e3 (1, 0, 0), (0, 1, 0) and (0, 0, 1) respectively F Function to be minimized to get Dupin indicatrix f Real-valued function H -1, I i curvature kl, k2, k3 hi,. ii curvature values at the nodes of a triangle L, M, N Coefficients in the mean curvature equation m, n Integers m Scaling factor of Dupin indicatrix Ni,p B-spline basis function of degree p n Normal vector P Point on the control polygon P Bezier coefficients of a cubic triangular patch p, q Integers Q Intermediate point Q Bezier coefficients of a quadratic triangular patch R' Linear or vector space R Radius of a sphere Ri,, Rational basis function of degree n r, s Integers r Circumradius of a triangle S Surface t Unit tangent vector along an arbitrary direction ti Direction represented by unit tangent vector tl, t2 Tangent directions corresponding to principal curvatures U, V Knot vectors u Point in the parametric plane u, v Vectors U, v Coordinates of a point in the parametric plane V Vertex V Vertex V, Velocity along normal direction V1, V2 Vectors vE, vy, vz Velocity along x, y and z directions W, X, Y Polynomial functions w Weight ., Constant coefficient x Point in Euclidean space x, y, z Coordinates of a point x Xc Circumcenter of a triangle Kn,i Normal curvature in direction ti Kt Normal curvature along direction t r1, K2 Principal curvatures I Angle between arbitrary normal section and tl Angle between an arbitrary plane and tangent plane ai Constant coefficient e Distance between tangent plane and plane of intersection Abstract of Thesis Presented to the Graduate School of the University of Florida in Partial Fulfillment of the Requirements for the Degree of 'i ir r. of Science THREE-DIMLN SIGNAL REPI:I.;S..\TATION OF EVOLVING INTERFACES IN COMPUTATIONAL FLUID DYNAMIC( ' By Abhishek Khanna December 2000 ('lI.ir: Wei Shyy '.. jor Department: Aerospace, Engineering '. li ii i, and Engineering Science livingng boundary problems arise in numerous important pl'ii -i ,i1 phenomena and they often form complex shapes during their evolution. liil iiii, ,1l methods for the solution of moving boundary problems primarily focus on interface tracking and on the solution of the field equations. The ability to track interfaces in two dimensions is well established. However, modification of the grid, representing the interface, as it evolves in a three-dimensional space introduces additional issues. In the current work, moving interfaces are represented by triangulated surface grids. Continuous spline surfaces are also used as an alternative means of surface representation. The unstructured triangulated grids are restructured and refined based on the shape and size of the triangular elements in the grid. As the interface deforms, points are automatically added and mesh connectivity modified to ensure that the accuracy of the interface representation remains consistent. Normals at the nodes of the surface mesh are required for the movement of the boundary and curvature pl1,i" an important role in problems involving capillarity. In the current work, methods for computing geometric features such as normals and curvatures on a triangulated surface are discussed. Curvature calculations are based on the mean of normal curvatures and on approximations to the Dupin indicatrix. Normals and curvatures are also computed on the parametric spline surface. This is done by the calculation of the partial derivatives of the coordinates of the surface nodes with respect to the parameter values. Results are presented to show the complex feature representation capability of the tracker code. A comparison between the different methods, used for computing curvatures, is made here. Filnlly, the effect of mesh quality and accuracy of the point location on the curvature error is also investigated. CHAPTER 1 INTRODUCTION Moving-boundary Problems o1 ving-boundary problems arise in a variety of important 1pl!-i, ,o1 phenomena. Some examples are materials processing, biofluid dynamics, flame propagation and fluid-structure interactions [1]. Such systems are often characterized by internal boundaries or by interfaces demarcating regions with different transport properties and mechanisms. Across these interfaces compositions, phases, material properties and flow features can change rapidly. The movement of interfaces is influenced by the flow field and the interfaces, in turn, affect the behavior of the flow. During their movement, the interfaces in:,, form complex shapes and in',:- change topologically while undergoing severe deformations such as dilatation, compression, fragmentation and collision. The main challenge to numerical methods for the solution of moving boundary problems is that the boundary location and shape must be determined as part of the solution. The evolution of the boundary is closely coupled with transport of mass, continuity, momentum and energy. Within each domain, field equations must be solved with the location of the internal boundary being determined simultaneously. The interface is often a discontinuity and there is need to track it accurately in both time and space within the limits of finite grid resolution. The challenges posed by fluid flow problems involving moving boundaries have resulted in a variety of computational techniques directed toward their solution [1]. These numerical methods can t'lpicOlly be classified as either Eulerian or Lagrangian [2]. Eulerian methods such as the level set method [3, 4] and the volume of fluid (VOF) method [5] do not keep track of the position and shape of the interface explicitly whereas Lagrangian methods such as the marker and cell (MAC) method [6] follow the interface explicitly, keeping track of both position and shape. While this makes the Lagrangian methods more complex, there is a significant advantage in knowing the interface position explicitly. In tracking distorted interfaces with accuracy, a combination of the strengths of both Eulerian and Lagrangian methods has been shown to be useful for problems involving fluid-fluid (Unverdi and Tryggvason [7]) as well as solid-liquid (Udaykumar et al. [8]) interactions. These methods employ fixed structured grids and afford the simplicity and availability of well-tested flow solvers. The only moving component is the interface which is a lower-dimensional surface overlying the fixed grid and is explicitly tracked. Solution of a moving boundary problem involves interaction of the flow solver with the interface tracking components. The effects of internal stationary as well as moving boundaries can be communicated to the flowfield by either incorporating changes in the control volume definitions that take place due to intersection of the interface with the background mesh (cut-cell approach [9]) or appropriate source term contributions (immersed boundary technique [10]). Information from the flowfield is transferred to the boundaries by interpolation. The ELAFINT (Eulerian-Lagrangian Algorithm For INt,-i f,,,- Tracking) method [8] is currently operational for the solution of two-dimensional problems. The current work is aimed at extending the ELAFINT algorithm by giving it three-dimensional interface tracking capabilities. Three-dimensional Interface Tracking Representation of interface in two dimensions is done by connecting a series of points describing the curves. Points on a curve always have two neighbors and thus storing the connectivity between the points does not have to be separate from storing the points themselves if the ordering of points is maintained appropriately. Addition of points to the curve is also straightforward as the new points can be inserted at appropriate locations into an already-existing list of points. In three dimensions, however, the task is much more difficult. To represent arbitrary deformations of a boundary, the surface is tesselated using triangles. In terms of data structures alone, there is no set ordering of points that sli--.t; itself. Thus, connectivity information cannot be stored implicitly if any nl::il iliti in changing the mesh is desired. The need to store connectivity information explicitly calls for the interface to be represented as a surface mesh. If the collection of points is stored as a structured grid, there is no nL::ililit:y in changing the mesh. Arbitrary addition or deletion of points is difficult when a complete regridding procedure has to be carried out at every time instant. The natural choice in such circumstances is to use an unstructured triangulated surface mesh to represent the surface. Using triangular-element unstructured mesh offers several advantages. The grid can be restructured. Individual triangles and sides can be bisected or rearranged to give new grid structures that better represent the evolving interface. Since the number of triangles meeting at a vertex is variable, increased accuracy in one region does not force unnecessary resolution in other areas. Triangles can cover a sphere more satisfactorily than rectangles without cusps or other local representation irregularities. Thus, free surfaces, complicated interfaces and boundaries of immersed objects can be represented accurately and economically. A general-connectivity triangulated mesh involves nontrivial bookkeeping tasks associated with the grid and its connections as all the data structures have to be updated during triangle, edge and vertex additions, reconnections and subtractions. The present approach uses unstructured grids to represent the interface while the background (stationary) grid is structured. The interface is initialized by means of a triangulation provided by standard CAD software. Since the initial connectivity information is available to us, the issues of unstructured grid generation, from discrete point data, involving complex geometries can be bypassed. After an initial interface position and shape is defined, the interface must move based on the velocity field. For the movement of the interface it is often required to know the normal at every node of the surface mesh. For the calculation of the surface forces which are distributed to the background structured mesh (.III1nl, "i by the field equation solver, information about the curvature at all surface points is needed. In this work, normal and curvature calculations are performed on the parametric spline surface as well as on the triangulated mesh and accuracy of these methods is discussed. While the generation of both two-dimensional planar as well as three-dimensional volume unstructured grids [11] is now not uncommon and largely straightforward, it is important to note that unstructured surface grid generation is not. Further, methods to adapt irregular surface grids by adding surface points to an already-existing set of points are still being developed. In order to track arbitrarily moving interfaces successfully, it becomes important to deal with triangulation of sets of points that are not D,1,,-iiil triangulated at every time instant. This capability is demonstrated by the current work. The interface tracking code has been tested on various cases documented here which reveal its versatality and accuracy. Overview of the Report The first chapter of the report presents an introduction to the topic of moving boundary problems and gives a review of the numerical methods used for solving these problems. Interface tracking issues which form a part of the solution of moving boundary problems have also been discussed here. ('l.iprtir 2 deals with the basics of B6zier and B-spline curves and surfaces that form the theoretical basis of the C1 surface spline algorithm discussed later in the report. ('liptri-. 3 presents the interface tracking algorithm. Algorithms for mesh refinement and parametric surface approximation have been discussed in detail. ,-t ri, iii used for computing normals and curvatures at the nodes of the interface have also been highlighted. ('lipr.r1 4 shows several examples which demonstrate the capabilities of the interface-tracking code developed so far. Conclusions drawn from the current work are presented in ( 'lhipirt 5. CHAPTER 2 REVIEW OF CURVES AND SURFACES This chapter reviews some of the modern methods used for representing surfaces in three-dimensions. Our particular interest is in the description and computation of surface features such as normals and curvatures for a general triangulated surface. The most popular methods for this purpose are currently Bezier and B-spline curves and surfaces [12, 13]. These parametric elements of curve and surface description have become the standard in geometric modeling. In this chapter we will detail the mathematical basis of these representations and point out their usefulness or otherwise vis-a-vis the accuracy of surface representation. Implicit and Parametric Forms The two most common methods of representing curves and surfaces in geometric modeling are implicit equations and parametric functions. An implicit equation of a curve in the xy plane has the form f(x, y) = 0. This equation describes an implicit relationship between the x and y coordinates of points lying on the curve. In the parametric form a curve can be represented as C(u) = (x(u), y(u)) where a < u < b. In this form, each of the coordinates of a point on a curve is represented separately as an explicit function of an independent parameter u. The parametric representation of a curve is not unique. It is instructive to think of C(u) = (x(u), y(u)) as the path traced out by a particle as a function of time where u is the time variable and (a, b) is the time interval. The first and second derivatives of C(u) are velocity and acceleration of the particle respectively. A comparison of implicit and parametric methods is given as follows: By adding a z coordinate, parametric method can be extended to represent arbitrary curves in three dimensional space, C(u) = (x(u),y(u),z(u)). The implicit form only specifies curves in xy (or xz or yz) plane. It is cumbersome to represent bounded curve segments (or surface patches) with implicit form. Boundedness is, however, built into the parametric form through the bounds on the parametric interval. The parametric form is more natural for designing and representing a shape in a computer. The coefficients of many parametric functions possess considerable geometric significance and explain many properties of these curves. The complexity of many geometric operations depends greatly on the method of representation. For example, it is difficult to compute a point on a curve or surface in the implicit form whereas in the parametric form, it is not directly possible to determine whether a given point lies on the curve or surface. In the parametric form, sometimes one has to deal with parametric anomalies which are unrelated to the true geometry. For example, poles in a sphere, where all the longitudnal lines meet, are critical points in a parametric form although geometrically they are no different from any other point on the sphere. Parametric curves possess a natural direction of traversal (from C(a) to C(b) if a < u < b) whereas implicit curves do not. From the above comparison, parametric approach appears to be an attractive option for representing geometry in a computer. The primary disadvantage, that is parametric anomalies, can be avoided by using parametric patches. B6zier Curves Polynomial B6zier Curves These are parametric polynomial curves described over a set of given points. An nt degree Bezier curve, with parameter u, is defined as C(u) = B,, (u)Pi 0 < u < 1 (2.1) i=O Here the basis (blending) functions {Bi,,(u)} are classical nth degree Bernstein polynomials given by Bi,(u) -= (1 u)n- (2.2) i!(n i)! The geometric coefficients of this form {Pi} are called control points. The polygon formed by the points {Pi} is called the control Pl.l,.;' Figure 2.1 shows a second- degree (n = 2) Bezier curve along with its control points. In any curve (or surface) representation scheme, the choice of basis functions determines the geometric characteristics of the scheme. Some important properties of basis functions can be listed as follows: Symmetry: For any n, the set of polynomials {Bj,,(u)} is symmetric with respect to the parameter value u = 1/2. PI P = C(O) P2= C(1) Figure 2.1: A second-degree Bezier curve along with its control points. * Bo,n(0) = B,,=(1) = 1 hence, a B6zier curve passes through its first and last control points. For an arbitrary set of points the curve does not pass through interior control points. * Nonnegativity: Bi,,(u) > 0 for all i, n and 0 < u < 1. * Recursive definition: Bi,,(u) = (1-u)Bi,,_l(u)+uBi_1,,_l(u) with condition Bi,,(u) = 0 if i < 0 or i > n. * Bi,,(u) attains exactly one maximum on the interval [0, 1]. * Partition of unity: Z2 Bi,,(u) = 1 for all 0 < u < 1. Partition of unity and nonnegativity property ensure a set of positive basis functions whose sum is equal to 1. From Equation 2.1, it is evident that points on a Bezier curve are convex barycentric combinations (Appendix A) of the specified control points. Therefore, these points lie within the convex hull of the control points of the corresponding Bezier curve. Derivative: B,,(u) = dB- = n(BlI,,n(u) Bi,,_1(u)) with conditions B_,,n_l(u) = 0 and B,,,_l(u) = 0. This property is used to compute the derivative of a Bezier curve. From the properties mentioned above, it is observed that the coefficients of a Bezier curve possess geometric flavor and give insight about the shape of the curve. The properties of Bezier curves are as follows: Transformation invariance: Bezier curves are invariant under the usual affine transformations such as rotations, translations and scalings, that is, one applies the transformations to the curve by applying it to the control polygon. Convex hull property: A Bezier curve is contained in the convex hull of its defining control points. This property ensures that the Bezier curve obtained from its control points is always bounded. A simple consequence of the convex hull property is that a planar control polygon always generates a planar curve. Variation diminishing property: No straight line intersects a given curve more times than it intersects the control polygon of the same curve (in three dimensions, replace 'straight line' with 'plane'). This property expresses the fact that a Bezier curve follows its control polygon rather closely and does not i,- l 1, more than its control net. Using the property of derivative of a basis function and definition of a B6zier curve, the general expression of derivative of a B6zier curve is given by C'(u) = d iB',,(u)Pi = n Bj,,_-(u)(Pij+ Pi) (2.3) du i=o i=0 The above equation is used extensively to compute geometric features such as normals and curvatures at points on a B6zier curve. These points on the curve can be calculated by making use of the recursive definition of basis functions. Denoting a general nth degree B6zier curve by C,(Po,... P,), the following relation holds true: C,(Po,...,P') = (1 u)C,_(Po,---, Pn-1) + uC-,_(PI,... ,P,) (2.4) Fixing u = u0 and denoting Pi by Po,i, a recursive algorithm for computing the point C(uo) = P,o(uo) on an nth degree Bezier curve is obtained as follows: Pk,i(uo) = (1 uo)Pk-l,i(u0) + UOPk-l,i+l(u0) k = ,...,n and i= ,...,n k (2.5) This is the de Casteljau Algorithm for B6zier curves (Figure 2.2). It is a corner cutting process which yields the triangular table of points shown in Table 2.1. The main attraction of the algorithm is the beautiful i ltil]1,-j, between geometry and algebra. According to the de C, .tfli. ii algorithm, points on a B6zier curve can be obtained by performing a sequence of linear interpolations. The geometry underlying the algorithm enables us to infer some important properties of B6zier curves discussed before. Transformation invariance is a direct consequence of the algorithm because linear interpolation is affinely invariant and so is a finite sequence of them. Convex hull property also follows because every intermediate point, Pk,i, is obtained as a convex barycentric combination of the previous set of points, {Pk-,j}, and at no step of the de C,,i t,-ii, algorithm do we produce points outside the convex hull of the specified control points. PP PP P1 1, 1 12 P2, 1 2, 0 P3, 0 P1, 2 P1, o Po 0 uo 1 P3 Figure 2.2: de C,,.],.jii algorithm for a cubic B6zier curve showing intermediate points generated by repeated interpolation at parameter u,. Rational B6zier Curves Although polynomials offer many advantages, there exist a number of important curve and surface types which cannot be represented precisely using polynomials such as circles, ellipses, hyperbolas, cylinders, cones, spheres etc. It is known from classical mathematics that all conic curves, including the circle, can be represented using rational functions which are defined as the ratio of two polynomials. In fact they are represented by rational functions of the form: X(u) Y(u) x(U) = y (u) = (2.6) W (u) W(u) X(u), Y(u) and W(u) are polynomials. For example, a circle of unit radius centered at the origin can be expressed as x(u) = (1- u2)/(1 +t 2), y(u) = (2u)/(1 + u2). On similar lines, we define an nt degree rational Bezier curve by E Bi, (u)., Pi C(u) = = 0 < u < 1 (2.7) i=0 Pi and Bi,,(u) are the control points and Bernstein polynomials as before while {.,' } is a set of scalars called weights. Thus, W(u) = E o Bi,,(u:.., is the common Table 2.1: Points generated by the de C,.tfli,.i algorithm. Po P1,0 P1 P2,0 P1,1 P2 Pn-1,0 .. P,o = C(uo) Pn-1,1 Pn-2 Pl,n-2 Pn-1 P2,n-2 Pl,n-1 Pn denominator function. Normally .. > 0 is assumed for all i which ensures that W(u) > 0 for all u [0, 1]. A rational Bezier curve can be written as C(u) = Ri,,(u)Pi (2.8) i=0 In the above equation {Rj, (u)} are the rational basis functions for this curve form and can be expressed as Ri,n(u) = B, u' (2.9) E Bj,,(u)wj j=0 The Ri,,(u) have properties very similar to those of the polynomial B6zier basis functions, Bi,,(u). In fact, if :, = 1 for all i then Ri,,(u) = Bi,,(u) for all i. Hence, Bi,,(u) are a special case of Ri,,(u) and polynomial B6zier curves are a special case of rational B6zier curves. B6zier Surfaces A curve C(u) is a vector valued function of one parameter. It is a mapping of a straight line segment into Euclidean three-dimensional space. A surface is a vector valued function of two parameters, u and v, and represents a mapping of a region of the uv plane into Euclidean three dimensional space. Thus a surface has a form S(u, v) = (x(u, v), y(u, v), z(u, v)). In this section we shall discuss nonrational as well as rational B6zier surfaces. Polynomial Tensor Product B6zier Surfaces A tensor product scheme is a bidirectional scheme that uses basis functions and geometric coefficients. The basis functions are bivariate functions of u and v which are constructed as products of univariate basis functions. Polynomial B6zier surfaces are obtained by taking a bidirectional net of control points and products of univariate Bernstein polynomials. They can be expressed as n m S(u, v) = Bi,n (u)Bj,.(v)Pi,j 0 < u, v < 1 (2.10) i=0 j=0 Properties of tensor product basis functions follow from corresponding properties of univariate basis functions. Properties of B6zier surfaces follow along the lines of the B6zier curves listed in Subsection 2. It is to be noted that unlike B6zier curves, there is no variation diminishing property for B6zier surfaces. The de Co, tfli., algorithm can be easily extended to compute points on a B6zier surface (Figure 2.3 [13]). Let (u,, Vo) be fixed. Fixing jo, Qjo(uo) = I=o B,,(Uo)Pi,, is the point obtained by applying the de C, .t,-li, i algorithm to the j, row of control points, that is, to {Pi,jo}, i = 0,..., n. Therefore applying the de Co, .tliii algorithm m + 1 times (jo = 0,..., m) yields Cuo(v) and applying it once more to Cuo(v) at v = vo yields Co(vo) = S(uo, vo). This process requires n(n+l)(m+1)/2+m(m+1)/2 linear interpolations. By symmetry, Co (u) can be computed first (n + 1 applications of de Co.ftli.,ii) and then Co(Uo) = S(uo, o). This, on the other hand, requires m(m- + 1)(n + 1)/2 +n(n+ 1)/2 linear interpolations. Thus if n > m compute Cv (u) first, then C o(u,). Otherwise, compute Cuo(v) first, then Cuo(vo). z S"w P2,0 Figure 2.3: de Co,-t.li]ii algorithm for a Bezier surface showing intermediate points generated by repeated linear interpolations at parameters (uo, Vo). Rational B1zier Surfaces Rational surfaces, similar to rational curves, are useful in representing specific shapes such as spheres, ellipsoids etc. which cannot be represented as polynomial surfaces. A rational Bezier surface with degree n in u direction and m in v direction can be defined as n m SBi,.(u)Bj,.(v'..,, P,j i=0j=0 S(u v) = 0 iE EBj )Bj - i=0j=o0 0 < u,v < 1 (2.11) {.: i} are the weights as before. The above given equation can be rewritten as shown below: n m S(u, v)= EE Rij(u, v)P, (2.12) i=o j=o {Rij(u, v)} are the rational functions but they are not products of other basis functions. They can be expressed in terms of Bernstein coefficients as follows: R,j (u, v) = ,(u)B,(v'. (2.13) E Br,,n()Bs,m(v)Wr,s r=0 s=0 Assuming .., > 0 for all i and j, the properties listed previously for nonrational Bezier surfaces (and product functions Bi,,(u)Bj,,(v)) extend naturally to rational Bezier surfaces. Furthermore, if .', = 1 for all i and j then Rij(u, v) = Bi,n(u)Bj,m(v) and the corresponding surface is nonrational. B-Spline Curves Curves consisting of just one polynomial or rational segment are often inadequate. Some of their shortcomings are listed below: A high degree is required to satisfy a large number of constraints. For example, an n 1 degree polynomial is required that passes through n given points. This task involves inversion of an n x n matrix and can be computationll'1: expensive if n is large. Hence, high degree curves are inefficient to process. They are also numerically unstable. * A high degree is required to accurately fit some complex shapes. Single segment curves and surfaces can be shaped by means of their control points and weights but the control is not local. For example, movement of one control point of a Bezier curve alters the location of all the points on the curve changing the shape of the entire object being represented by the Bdzier curve. A possible solution is to use piecewise polynomial or rational curves and surfaces which join with some level of continuity. B-splines are attractive in this regard. Their piecewise polynomial or rational basis functions have nice analytic properties which follow along the lines of Bdzier basis functions. This ensures that B-Splines have good geometric properties similar to their Bdzier counterparts. Polynomial B-Spline Curves A Bdzier curve is a one segment curve and therefore the degree of a Bdzier curve is one less than the number of control points. However, since the B-spline curve is a piecewise curve, the degree of the B-spline segments has to be specified. A pth degree B-spline curve is defined by C(u) = N ,pPi a < u < b (2.14) i=O Here {Pi} are the control points and {Ni,p(u)} are the pt^ degree B-spline basis functions defined on a knot vector U = {a,...,a, p, +,..., Um-p-_, b,..., b} as P+1 p-+ 1 if ui < u < Ui+ 0 otherwise N -- Ni ( i)+ p+1 -- I Ni,p -= N,pl(u) + u l-uN1,p_1(u) (2.15) ti+p -- Ui Ui+p+l -- i+1 Hence, computation of a set of basis functions requires specification of a knot vector, U, and degree, p. The basis function equation can yield a quotient 0/0 which is defined to be zero. The half open interval, [ui, ui+l), is called the ith knot span and it can have zero length since knots need not be distinct. Three steps are required to compute a point on a B-spline curve at a fixed given parameter value. They are enumerated as follows: 1. Find the knot span in which parameter u lies. 2. Compute the nonzero basis functions. 3. hiill iply the values of the nonzero basis functions with the corresponding control points. {N,p (u)} are piecewise polynomials defined on the entire real line. Some properties of these basis functions are listed as follows: Local support property: {Ni,p(u) = 0} if u is outside the interval [ui, ui+-p+). This is illustrated by the following triangular scheme: N1,0 N0,2 \ N1,1 No,3 N2,0 N1,2 N2,1 N1,3 N3,0 N2,2 N3,1 N2,3 N4,0 N3,2 Notice that N1,3 is a combination of N1,0, N2,0, -V,,, and N4,0. Thus N1,3 is nonzero only for u E [Ul, u5). * In any given knot span, [uj, Uj+l), almost p + 1 of the N,p are nonzero, namely functions Nj_,p,,..., Nj,p. This is shown as follows: N2,0 N3,0 N4,0 Ni,i No,3 N1,2 N2,1 N1,3 N2,2 N3,1 N2,3 N3,2 N4,1 N3,3 On [3, U4) the only nonzero zeroth degree function is .V,,. Hence, the only cubic functions not zero on IU3, U4) are N0,3, ..., V, . * Partition of unity: For an arbitrary span, [ui, ui+1), i-p Njp(u) all u E [u, Ui+l). 1 for * Except for the case p = 0, Ni,,(u) attains exactly one maximum value. * Nonnegativity: N,p > 0 for all i, p and u. * Derivative: The derivative of a basis function has a recursive definition where N',p = N u) N p u) P i+p--ji i+p+l--i+l * A knot vector of the form U = {0,..., 0, 1,..., 1} yields Bernstein polynomials p+1 p+1 of degree p. If there are m + 1 knots, there are n + 1 basis functions where m = n + p + 1. Also, .,,(a) = 1 and N,,,(b) = 1. S.1,1 iy properties of B-spline curves follow from those of their basis functions. Some of them are listed as follows: Transformation invariance: B-spline curves, like Bezier curves, are invariant under affine transformations such as rotations, translations and scalings. If n = p and U = {0,...,0,1,...,1} then C(u) is a Bezier curve. +1- p-+1 Endpoint interpolation: C(0)= Po and C(1) = P. Convex hull property: A B-spline curve is contained in the convex hull of its control polygon. Local modification scheme: '.1,oving Pi changes C(u) only in the interval [u, Ui+p+l). This follows from the fact that Ni,p(u) = 0 for u u [ui, Ui+p+,). '.loving along the curve from u = 0 to u = 1, {Np,,(u)} functions act like switches. As u moves past a knot, one Ni,,(u) (and hence the corresponding Pi) switches off and the next one switches on. Variation diminishing property: No plane has more intersections with the curve than with the control polygon (for two dimensional curves, replace the word plane with line). Derivative of a B-spline curve C(u) = E- o' Ni,p(u)Pi defined on knot vector U = {0,..., +,..., Um-p-, 1,...,1} is given by p+1 p+1 C'(u) = PpNip- (u) (2.16) i=0 Ui+p+1 Ui+1 The basis functions in the above equation have been computed on knot vector U' obtained by dropping the first and last knots from the original knot vector U, that is, U'= {0,...,0,Up,...,Um-p-2,1, **,1} P P Nonuniform Rational B-Spline Curves NonUniform Rational B-Splines (NURBS) have become the industry standard for the representation, design and data exchange of geometric information processed by computers. The enormous success of NURBS can be attributed to the fact that they provide a unified mathematical basis for representing both analytic shapes such as conic sections as well as free-form entities such as car bodies and ship hulls and this feature makes NURBS an attractive tool for representing geometry in moving boundary problems. These elements are the general form of polynomial B-Splines and therefore polynomial Bezier elements. NURBS curves are defined on a nonuniform knot vector. A pth degree NURBS curve is defined by E N,p(u'.- Pi C(u) = a < < b (2.17) i=0 {Pi} are the control points (forming a control polygon), {.. } are the weights and {N,p(u)} the pth degree B-spline basis functions defined on a nonperiodic knot vector U = {a,..., a, Up+1,..., Umn-p-1, b,..., b}. Unless otherwise stated, a = 0, b = 1 and p+1 p+- .', > 0 for all i. A NURBS curve can also be rewritten as C(u) = E Ri,p(u)Pi (2.18) i=O {Rp(u)} are rational basis functions. They are piecewise rational functions on u G [0, 1] and can be expressed in terms of polynomial B-spline basis functions as Ri,,(u) = ,' "' (2.19) F' Np,,(u) wj j=0 {Ri,p(u)} have properties similar to those of {Ni,p(u)} hence properties of NURBS are no different from those of polynomial B-spline curves. B-Spline Surfaces Like polynomial Bezier surfaces, polynomial B-spline surfaces are tensor product patches of their basis functions. Rational B-Spline surfaces, on the other hand, are obtained by using piecewise rational functions. Polynomial B-Spline Surfaces A polynomial B-spline surface, just like a polynomial Bezier surface, is a tensor product structured patch. It is obtained by taking a bidirectional net of control points, two knot vectors and products of the univariate B-spline basis functions as shown in the following equation: n m S(U, V) = NE E N,p(u)N,q(v)Pi, (2.20) i=0 j=0 This surface is defined on knot vectors U = {0,... 0, +, ..., Ur-p-l, 1,..., 1} p+1 p+1 and V = {O,...,O, qv+,..., vq_-,1,...,1}. U has r+ 1 knots and V has s+ 1. q+l q+1 Hencer = n + p + hands = m +q +1. Five steps are required to compute a point on a B-spline surface at fixed (u, v) parameter values: 1. Find the knot span in which u lies, say u E [ui, Ui-). 2. Compute the nonzero basis functions Ni_,,,(u),..., Ni,(u). 3. Find the knot span in which v lies, say v E [vj, Vj+-). 4. Compute the nonzero basis functions Nj_,,,q(),..., Nj,,(v). 5. ilt Iply the values of the nonzero basis functions with the corresponding control points. Properties of tensor product basis functions follow from corresponding properties of univariate basis functions and properties of B-spline surfaces are similar to those of B-spline curves. The only exception to this is the variation diminishing property which is not known to exist for a B-spline surface. Nonuniform Rational B-Spline Surfaces A NURBS surface of degree p in the u direction and degree q in the v direction is a bivariate vector valued piecewise rational function of the form: n m E N,,p(u)Nj,q(V'.' iPi,j S(u,v) = i=j= m 0< u,v < 1 (2.21) E E Np(U)Njq(V). i=0j=0 {Pi,j} form a bidirectional control net, {.'. i} are weights, {NA,p(u)}, {Nj,q(v)} are nonrational basis functions on knot vectors U = {0,... 0, U, +, ..., Ur-p-,, 1,..., 1} p+1 p+1 and V = {0,...,0, V ..., _q_, 1,...,1} where r = n +p+ 1 and s = m + q+ 1. q+l q+l The surface can be rewritten as n m S(u, v) = E R,(u, )P, (2.22) i=o j=o Here {Ri,j(u, v)} are the rational basis functions. They are not a tensor product of the polynomial B-spline basis functions and are related to these functions as shown below: Ni,p(U)Nj,q(v':,: Ri,j(u, v) =n p(u (2.23) E E Nk,,p(u)NiN,(, ,(,v k=01=0 Assuming :, > 0 for all i and j, the properties listed previously for nonrational B-spline surfaces (and product functions N,p(u)Nj,q(v)) extend to rational B-spline surfaces. Furthermore, if : i = 1 for all i and j then Rij(u, v) = N,p(u)Nj,q(v) and the corresponding surface is nonrational. B6zier Triangles Bezier triangles are a direct generalization of Bezier curves and the de C',.1ii',.l algorithm for triangular patches is an extension of the corresponding algorithm for curves. The control net now is of a triangular structure. The unstructured nature of the control net opens the possibility of using Bezier triangular patches along with an existing surface triangulation for representing an interface in a moving boundary problem. Tensor product surface patches cannot be used in this case because the surface information is stored in the form of an unstructured mesh. Polynomial B6zier Triangles A polynomial Bezier triangle, as mentioned previously, is defined on a control net that has triangular structure. Figure 2.4 gives an example of a cubic patch (shaded patch) with its control net. The numbering scheme used for control points is also shown. P030 P. p S. P.,,,, P003 Figure 2.4: A cubic Bezier triangle with its control net. With n as the degree of the patch, the control net consists of (n + 1)(n + 2)/2 vertices. The numbers (n + 1)(n + 2)/2 are called triangle numbers. An n" degree Bdzier triangle is defined as 0 < u,v, 1- u -< 1 Here, the point Pijk has been denoted by Pi. i means i + j + k = n, always assuming i,j, k > 0. u in the parametric plane having barycentric coordinates the bivariate Bernstein polynomials defined by (2.24) Si + + +k, therefore, i = n = (u,v, 1- -v) is a point (Appendix B). {Bi, (u)} are Bi,(u) = niv(1 u v)k 1il = n (2.25) Bi, (u) = 0 if i, j, k < 0 or i, j, k > n for some subscripts. This follows standard convention for trinomial coefficients. Some properties of bivariate polynomials are stated as follows: S(u) = 7 Bi,(u)Pi If=n Endpoint interpolation: A Bezier triangular patch passes through endpoints of the control net as Boon0, = Bono,n = B oo0, = 1 for all 0 < u, v, 1 u v < 1. Nonnegativity: Bi,(u) > 0 for all i, n and 0 < u, v, 1 u v < 1. Partition of unity: = Bi,n(u) = 1 for all 0 < u, v, 1 u v < 1. Recursive definition: Taking el = (1, 0, 0), e2 = (0, 1, 0) and e3 = (0, 0, 1), Bi;,(u) = uBi-el,n-1(u) + vB,-e2,n-(u) + (1 u v)Bi-e3,n-l(u) relation holds for all Bernstein polynomials when i = n. .I ,iy properties of triangular Bezier triangles follow from those of the bivariate Bernstein polynomials. Some of them are listed below: Transformation invariance: Bezier triangles are invariant under usual affine transformations such as rotations, translations and scalings. Boundary curves: For a triangular patch, boundary curves are determined by control vertices at the boundaries (having at least one zero as a subscript). For example, a point on the boundary curve P,(u, 0, 1 u) can be obtained by P,,i(u, 0, 1 u) = uPr-1,i+el + (1 u)P,-1,i+e3, which is the de Coi-t, li"l algorithm for Bezier curves. Convex hull property: This property is guaranteed since for 0 < u, v, 1 u- v < 1, each P,.; is a convex combination of the previous Pr-1,i. The de Co, -tli., 1 algorithm for triangular patches is a direct generalization of the corresponding algorithm for curves (Figure 2.5). Figure 2.5: de Co.t(-liiil algorithm for a B6zier triangle showing the intermediate points generated. Like the curve algorithm, the triangle algorithm uses repeated linear interpolation. Given a triangular array of points, {Pi}, and i = n, the recursive algorithm follows: Pr,i(u) = uPr-l,i+el(u) + vPr-l,i+e2(u) + (1 U v)Pr-l,i+e3(u) r = ,...,n and i =n- r (2.26) Here Po,i = Pi and P,,o is the point with parameter value u on the triangular Bezier patch P,(u) = S(u). As stated before, the triangular de C,,f-tli,,ii algorithm is completely analogous to the univariate one. While discussing derivatives for tensor product patches, partial were considered because they are easily computed for those surfaces. The situation is different for triangular patches as the appropriate derivative here is the directional derivative. Let ul and u2 be two points in the domain. Their difference d = u2 ul defines a vector. The directional derivative of a surface at x(u) with respect to the direction vector, d = (d, e, f), is given by Ddx(u) = dxu(u) + exv(u) + fxlu- (u) (2.27) A geometric interpretation: in the domain, draw the straight line through u that is parallel to d. This straight line is mapped to a curve on the patch. Its tangent vector at x(u) is the desired directional derivative. The partial of a Bezier patch are easy to compute and the directional derivative is given by DdP(u) = n [dPi+el + Pi+e2 + fPi+e3]Bi,n-(u) (2.28) i|=n-l Thus a directional derivative is obtained by performing one step of the de C', -t -l i, ,l algorithm with respect to the direction vector d and n 1 more with respect to the position u. In general, for an rt directional derivative, one has to perform r steps of the de C', .ftlil, algorithm with respect to d and the rest n r with respect to u. l i::i- derivatives can also be calculated this way by applying the de C, .t-.li ji algorithm along different directions. The order in which the operations are performed is irrelevant. Rational B1zier Triangles Following the familiar theme of generating rational curve and surface schemes, a rational Bezier triangle is defined as E Bi,n(u)wiPi s(u) = (2.29) E Bi,.(u)wi Here, as usual, {wi} are the weights associated with the control vertices Pi. For positive weights, the properties of polynomial Bezier triangles extend to rational Bezier triangles and if all weights are equal to one, the two forms are exactly the same. This chapter presented methods for representing surfaces, in a parametric form, in three dimensions. Bezier and B-Spline tensor product surface patches result in smooth representations of geometric shapes. However, these are structured patches and cannot be used for representing an unstructured triangulated surface available to us. B6zier triangular patches require a control net that has a triangular structure. This requirement poses a restriction as it is not be possible to break up an an arbitrary surface triangulation, where arbitrary number of triangles meet at every node, into triangular control nets. Therefore, the methods discussed so far are not directly applicable to us. We need an algorithm that uses the available surface triangulation to generate control nets for the B6zier triangles, such that these triangular patches join with some level of continuity. CHAPTER 3 INTERFACE TRACKING METHOD The challenges posed by moving boundary problems have resulted in a number of computational techniques directed towards their solution. 1,r lil I, based on the combined Eulerian and Lagrangian approach require the interaction of the flow solver with the interface tracking components. The flow solver solves the governing flow equations over a stationary background mesh. The velocities computed over the cells of the background mesh are used to determine the velocity field of the moving surface. The tracker keeps track of the shape and position of the moving interface. It also computes geometric information such as normals and curvatures for the flow solver. This chapter discusses a three-dimensional surface tracker, originally developed by .1,i.',liII, et al. [14], for tracking interfaces in moving boundary problems. The moving interface is represented by a triangular-element unstructured mesh which offers flexibility in changing the mesh. Individual triangles, which do not satisfy the mesh quality criteria, can be bisected or rearranged to give new grid structures which better represent the interface. The purpose of the tracker is to refine and modify an existing mesh as it evolves in time in a moving boundary problem, hence, it does not include the unstructured grid generation methods used for triangulating a discrete set of points in space. The interface tracking algorithm is presented in Figure 3.1 and a description of its steps is given in the rest of the rest of the chapter. Read control file Read grid data and establish surface mesh connectivity Set surface mesh quality criteria ^ No t < tmax YYese Generate control mesh from surface mesh data Refine surface and control mesh Update surface mesh Calculate velocities for surface mesh Move surface mesh Print surface mesh information Destroy control mesh Increment time Destroy surface mesh Figure 3.1: Interface tracking algorithm. The interface tracking algorithm starts by reading the number of time steps and the time step size value from the control file. In the next step, the initial unstructured surface triangulation data is applied to the tracker. This unstructured mesh, initializing the moving interface, is obtained from a standard CAD software. The number of nodes and their coordinates, number of triangles and the node numbers of every triangle are required to establish the surface mesh connectivity. With the information given above, we develop the node, edge, triangle and surface mesh data structures described in Appendix C. Having provided an initial grid information, a mesh quality criterion is set. The tracker refines I-, tliiiIhl by breaking them into smaller triangles. A triangle is considered large enough to be bisected if its area exceeds the ,,.vi,... triangle area" times an -.,. '. ti, ll" (which is set in this step). In order to get accurate results from the governing fluid flow equations, it is desirable to have triangles having a uniform size and a good aspect ratio. The length of edges of the surface mesh should t'l-pilly be of the order of the grid spacing of the background field equation solver mesh. Thus the value of the average triangle area is computed on the basis of the size of the background Cartesian cell. For now, in order to test the interface tracker alone, the average triangle area is the mean of areas of all the triangles in the mesh. The first step in the time loop is generation of a control mesh from surface mesh data. This is done using the C' Surface Spline algorithm [15] presented later in this chapter. The motive behind using a surface modeling scheme is to use the control mesh to generate a smooth surface which gives accurate normal and curvature values at every node. M-..1- refinement according to the size and shape of the triangles of the surface mesh follows next. A restructuring of the mesh does not generally involve movement of any of the vertices but does include many triangle, vertex and edge additions, reconnections and subtractions. As mentioned previously, the surface tracker bisects big triangles by breaking them up into smaller triangles. A triangle has a distorted shape if its circumcircle contains an opposite node of an ,dji;,p-.~t, triangle. This is called the circumcircle criterion and the edge common to the two triangles is a l. I i- l,i" that has to be swapped. The details of these mesh refinement operations are included in this chapter. Once mesh refinement is complete, the surface mesh is updated. This involves computations of normals and curvatures. The interface points move in the direction of their normals with a velocity which is obtained from the field equation solver. For the calculation of surface forces to be later distributed to the background structured mesh employed by the solver, information about the curvature at all the points is required. Details of normal computations and various methods used for computing curvature on the triangulated surface mesh as well as on the parametric surface obtained by the application of the C' Surface Spline algorithm are given in this chapter. In order to evolve the interface in time, velocity is required at all the nodes of the mesh. For testing purposes, the tracker is currently working with user-input velocity functions. With the velocity field available and normals known at all points, the interface is advected by a distance equal to velocity times the time increment value, obtained from the control file, in the normal direction. The new interface position can be printed by storing the grid data in the output file. If required, normal and curvature values can also be plotted. This completes the time cycle. The control mesh, generated by the C1 Surface Spline algorithm, is destroyed to be created again during the next cycle according to the new positions of the surface mesh points. The iteration counter is incremented before the control goes to the next cycle. C' Surface Spline Algorithm The generation of control mesh is done using the C1 Surface Spline algorithm proposed by Peters [15]. The method is applicable to irregular meshes with non- quadrilateral cells and with arbitrary number of cells meeting at each mesh node. Arbitrary freer-form surfaces are allowed to be modeled in the same conceptual frame work as tensor product B-Splines. In this report, however, the scheme is discussed for a triangulated surface. The algorithm requires mesh points, mesh connectivity and blend ratios (knot spacings) as input. The output is a set of three-sided, cubic patches expressed in the Bernstein-Bezier form [16]. The mesh points serve as the control points of a smooth piecewise polynomial surface representation that is local, evaluates by averaging and obeys the convex hull property. The algorithm defines a spline space over irregular meshes by local averaging as do the other algorithms by Peters [17, 18]. The main improvement of the present approach over the earlier algorithms is in the number of triangular patches generated per mesh triangle. While earlier schemes generate 48 three-sided Bernstein-Bezier patches for each mesh triangle, the current algorithm reduces that number to 12. The C1 smooth surface can be generated by applying the triangular de C,. tli. i algorithm on each of these patches. The representation obtained has the following properties: Free form modeling capability: There are no restrictions on the number of triangles meeting at a mesh node. Smoothness: Surface splines form smooth surface parameterizations and in order to manipulate surface splines, it suffices to add, subtract or move mesh points locally. Evaluation by averaging: The coefficients in Bernstein-B6zier form can be obtained by performing averaging operations on the input mesh. The algorithm is local and can be interpreted as a rule for cutting an input mesh such that the control mesh obtained gives a spline surface. Convex hull property: The surface lies locally and globally in the convex hull of the input mesh. Intuitive shape parameters: Averaging process is geometrically intuitive. Blend ratios, analogous to knot distances, govern the depth of cuts into the input mesh to generate the control mesh. Smaller cuts result in a surface that follows the input mesh more closely and changes the normal direction more rapidly across the boundary. Low degree parametrization: The surface is parametrized by low degree polynomial patches and the representation can be extended to rational patches. A six step algorithm for generating a surface spline control mesh is given below. The first three steps generate a refined mesh consisting of quadrilateral cells such that each original node of the unstructured surface mesh is surrounded by vertices of degree four. The last three steps generate the surface parametrization. In the following steps, all the coefficients A, C, P, Q, V are points in space. All indices are interpreted modulo n where n is the number of quadrilateral cells, edges or Bernstein- Bezier patches meeting at a node. Step 1: Mesh Refinement. Referring to Figure 3.2, for each edge with node vertices Vi, Vj, add an edge vertex Vij = (Vi + Vj)/2 and for each triangle having node vertices Vi, Vj, Vk, add a triangle vertex Vijk = (i + Vj + Vk)/3. Connect triangle vertex to all edge vertices of a triangle. This leads to the creation of a refined mesh of quadrilateral cells. Each triangle now contains three cells. All the edge vertices have a degree four whereas triangle vertices have degree three. Vk Vij k Vi Figure 3.2: ',,.ii refinement step where the addition of edge vertices and triangle vertices leads to creation of refined mesh of quadrilateral cells. Step 2: Edge Cutting. For every cell ci, starting with the node vertex that corresponds to a point on the surface mesh and taking vertices VI, V2, V3, V4 (Figure 3.3), define a preliminary cell center Ci subject to blend ratios 0 < ail,ai2 < 1 as follows: Ci = (1 ail)(1 ai2)Vl + (1 -- ail)ai2 + ailai2V3 + ail(l -- ai2)V4 (3.1) \::,L derive cell centers C1 from Ci such that all (Ci+Ci+1)/2 surrounding a vertex lie in the same plane. For each vertex V surrounded by less than five cells, Ci = Ci. When the number of cells, in anticlockwise order, is greater than four (cl,..., c.), the cell center corresponding to ci is given by n C = V + 2.: /n cos(27rj/n)C j=1 (3.2) Here V = l/n Ei' Ci and 0 < .:, < 1. Also, .., and ,- = 1/(2cos(ir/n)) if n is odd. 1/(1 + cos(27r/n)) if n is even Figure 3.3: Edge cutting step where edges of the surface mesh are cut according to the blend ratios. Step 3: Quadratic Meshing. For each vertex V surrounded in order by cell centers C1,..., Cn, compute the final vertex location V = 1/n Y-, Ci. Note that Fi_' Ci = F_', C'. For each edge separating two cells with centers Ci_1, Ci, create an edge coefficient Ai = (Ci_- + Ci)/2. This step leads to a refined mesh such that three points are associated with each cell edge (Figure 3.4). These can be interpreted as the Bernstein-Bezier coefficients of a quadratic curve segment. Note that each vertex V lies in the same plane as the surrounding Ai's. Ci-1 Ai Figure 3.4: Quadratic meshing step which yields Bernstein-Bezier coefficients of a quadratic curve segment. Step 4: Quadratic Patching. Each quadrilateral cell is made of four second- degree Bernstein-Bezier triangular patches as shown in Figure 3.5. For each vertex V set Qo2o,i = V, QlOi = Ai where i = 1,..., n and compute the interior coefficients by averaging these values. Set coefficients Qo1,i = Qilo,i+i = (Qllo,i + Qllo,i+1)/2 and Q002 = (Q01,i + o011,i+1)/2. 8 neighbors 6 neighbors Q 020, i-1 Q 200, i+l o101, i+1 . S1:002 Q 110, i+l 011, i 101, i . 011, i+1 = Q020, i Ai Ql10, i Q200, i 2n neighbors 8 neighbors Figure 3.5: Quadratic patching step where internal coefficients are computed to give quadratic Bezier triangle patches. Step 5: Degree Raising. The quadratic patches are reinterpreted as cubic patches using the following formulas: Po30 Q020 p21 Qo20+2Qoll P021 3 P120 p012 Qo20+2Qllo Qoo2+2Qoll 120 012 3 3 p111 p _Qiio+Qioi+Qoii Q) P111 P003 ---+-3 -o+Qo Q002 p p Q200o+2Qlo Qoo2+2Qiol 210 102 3 3 p,21 Q200+2Qo10 P201 3 P300 Q200 The coefficients P111, P102, P012 and POO3 are recomputed in the Step 6. Step 6: Twist Adjustment. At vertex V = Po30 = Qo2o surrounded by n quadrilateral cells, set c = cos(27r/n). Using Ci, the center of ith cell, compute coefficient P11,i = [(2 c)Q10,i + (1 + c)Q200,i + 2Ci + Qo20,i]/6. Finlly, counting with j around the four-neighbor vertex Qoo2, P102,j = P, -l i+1 = (P111,j + Pll,j+1)/2, Poo3,j = (P102,o + P102,2)/2. Steps 5 and 6 are not necessary if n = 4. This step is required for getting C1 continuous Bezier triangular patches. The surface obtained by application of the surface spline algorithm is parametrized by Bernstein-B6zier patches that can be, depending on users choice, either rational or polynomial. There are 12 Bernstein-B6zier triangular patches generated for evry triangle of the unstructured surface mesh. Three sided patches ensure tangent plane continuity of the surface. Mesh Refinement As the interface evolves in time, the positions of the mesh nodes change which might lead to large or skewed triangles having a large aspect ratio. For accuracy of the numerical computations, the spacing of the interface points should be of the order of the background Cartesian mesh spacing. Hence, it becomes important to refine the mesh periodically and retain an accurate representation of the interface. This is done according to the mesh refinement algorithm shown in Figure 3.6. The algorithm starts by first creating a list of -IW, ti,,.~ilr"., that is, a list of triangles whose area is greater than a specified value. As long as the list is not .llln , the algorithm breaks up the big triangles into smaller parts by adding new points and it works on the '-i2-'2-;t triangle first. By experience, this results in a better quality mesh that gives accurate curvature results as compared to results obtained by the refinement of any arbitrary "II.5 liiiiiil Once all the 'II; tliiiiii.." are broken into smaller parts resulting in a surface made up of triangles having a uniform size, the tracker works on the skewed triangles. Diagonal swap method, which is essentially a mesh reconnection method, is .lii'1ii,' for this purpose. The details of the point insertion and diagonal swapping procedures are discussed in detail in this section. Start Create a list of big triangles on the basis of the maximum area criteria Is list empty? Yes No Pick the biggest triangle from the list Do point insertion Check new triangles formed for quality and update the list of big triangles Create a list of bad edges on the basis of the circumcircle criteria Yes -Is list empty? No Pop an edge from the list Do diagonal swapping Check new edges formed for quality and update the list of bad edges Figure 3.6: 1, !i refinement algorithm. \-.-- points are inserted into "I'i ti I. i,,ll', whose area exceeds a limit specified previously (based on the average triangle size and an area factor), according to the Bowyer-Watson point replenishment algorithm [19, 20] which was concieved for two- dimensional D,-1,Lililil triangulation [21]. Once all the big triangles are bisected into smaller ones, the refinement algorithm works on the distorted triangular elements that do not satisfy the circumcircle criterion. "Bad I,. l" are identified and diagonal swap operation is performed on the quadrilateral formed by two triangles sharing a bad edge. The details of point insertion and diagonal swap technique are presented below. Point Insertion Insertion of new point is done to break a big triangle into smaller triangles. The algorithm followed is shown in Figure 3.7. Start Calculate coordinates of new point to be inserted Delete triangles using the circumcircle criteria and form insertion cavity Reform cavity and form new triangles by connecting the new point with the nodes making the insertion cavity Figure 3.7: Point insertopn algorithm. Figure 3.7: Point insertion algorithm. To begin with, the original grid data of the interface includes coordinates of the nodes and node numbers of every triangle in the mesh. The connectivity information provided is of topological nature and the exact geometry of the mesh (how nodes in the mesh are connected to their respective neighbors) is not available to the tracker. Connection of nodes by straight lines is actually an assumption of the nil,.i" surface. In such circumstances, it is not possible to find a unique location of the new point to be inserted. The tracker computes normals and curvatures at all nodes of the mesh for the field equation solver and the results of these calculations, especially curvature calculations, are very sensitive to the coordinates of the new point. Placement of the new point in the plane of the big triangle leads to zero curvature at that location whereas slight displacement of the point from the plane causes rapid variation of curvature. Keeping this in mind, the new point is inserted based on the local curvature information obtained from its neighborhood. The new point is made to lie on a sphere which passes through the nodes of the big triangle and whose curvature is equal to the average of the curvature value at the same three nodes. Let k1, k2 and k3 be the curvatures at the node locations of the triangle under consideration. Then k = (kl+1! ".+ :"')/3 is the curvature at the new point location. If R is the radius of the sphere on which the new point is located then R = 1/ I I I Let r be the circumradius of the triangle being refined. Since the sphere passes through the nodes of this triangle, the normal to the plane containing the triangle is in the same direction as the line that passes through center of the sphere and circumcenter of the triangle. The new point has to be placed so that it lies on the sphere, hence, it should be projected by a distance d = R- R2 -r2 along the normal from the circumcenter. The point can be projected above or below the triangle because theoretically there are two spheres of radius R passing through three given points. The appropriate choice is made on the basis of sign of local surface curvature. The advantage of using a circumcenter is that, after projection, the new point is equidistant from all the nodes of the triangle. With this method of calculating the coordinates of the new point, three cases arise which are enumerated as follows: 1. k = 0 and as a result the radius, R, tends to infinity. This implies that the local surface around the new point is planar. Hence the new point is inserted at the circumcenter of the triangle. 2. R > r and in this case the point placement is done by projecting a point along the normal from the circumcenter by a distance d. 3. R < r which means that the surface triangle lies in a region of very high curvature (since R = 1/111 ||). In such a situation, the new point is placed at the circumcenter of the triangle. Once the coordinates of the new point are available, the subsequent step in the point insertion algorithm is the deletion of triangles according to the Bowyer-Watson criterion and formation of an "illi'i tion (,i it ". The details of these operations are given below. The Bowyer-Watson algorithm is based on the circumcircle property which ensures that no point of a D.1,]ii i,;' triangulation can lie within the circle circumscribed to any triangle. It is essentially a reconnection method since it computes how an existing triangulation is to be modified because of insertion of a new point. The Bowyer-Watson algorithm is a right choice for the current application as it is required to maintain an accurate representation of a deforming interface by adding points where necessary. Even though the algorithm is for two dimensional triangulations, its modification is being used for three-dimensional surfaces. The circumcircle criterion is illustrated in Figure 3.8. AABC is a big triangle that needs to be refined because of its large size. Initially, its edges AB, BC and CA are "iI',i" because each of them is associated with two triangles of the mesh. Since AABC is to be bisected, it is marked for deletion and all its edges are declared "dti l in, lt because now each of them is a part of only one triangle of the surface mesh. Point F is the new point inserted to break this triangle into smaller triangles. Distance between the circumcenter of AACE (point H) and F is less than the circumradius of AACE hence that triangle is marked for deletion. Edges CE and EA change state from alive to dormant whereas edge AC, which was dormant previously, now becomes "(t., Il". AADB satisfies the circumcircle criterion discussed above and is not marked for deletion. E D H G* C Figure 3.8: Bowyer-Watson criterion for marking triangles to be deleted in order to form an insertion cavity. Each time a new point is inserted into the mesh, a search is carried out starting from the 'illi-,'tion tiill i"h which is marked for deletion. The search then travels across the three edges of this triangle to its neighbors. If the neighboring triangles fail to satisfy the circumcircle criterion discussed above, they are marked for deletion. Any edges common to two deleted triangles are also identified (declared dead) to be removed later. The recursive process continues up to a specified number of levels which, by experience, is chosen such that it covers all the triangles whose circumradius is less than the distance between its circumcenter and the new point. Once triangles, edges and nodes to be deleted are identified, an iill' tion (,I':, it is constructed. The process is depicted in Figure 3.9 where AABC, AADB, ACBE and edges AB, BC are deleted according to the Bowyer-Watson criterion. Edges AD, DB, BE, EC and CA are in a dormant state because each of them is a part of only one triangle of the surface mesh. These edges form the insertion cavity. F is the new point to be inserted and new triangles are formed by connecting this new point with points A, D, B, E, C that form the insertion cavity. This completes the point insertion algorithm which results in triangles having relatively uniform sizes. Diagonal Swapping This mesh restructuring technique is used after point insertion is complete and none of the triangles have area larger than the set criterion. Diagonal swapping is basically attempted to mend distorted triangles which do not satisfy the circumcircle criterion. The procedure can be clearly understood with the help of Figure 3.10. The distance between point C and the circumcenter of AABD is less that the circumradius of AABD. Hence, edge BD, common to the two triangles under consideration, (a) (b) (c) Figure 3.9: Bowyer-Watson algorithm: (a) Existing triangulation; (b) Formation of insertion cavity; (c) Reconnection of an existing grid around a newly inserted point. becomes a 'ldiIiI ,.li." that is to be swapped. A diagonal swap connects A to C and the newly formed triangles, AABC and AACD, have a better aspect ratio than AABD and ABCD in the original triangulation. \D B (a) (b) Figure 3.10: Diagonal swapping operation: (a) Triangles before diagonal ,.-,I,]ili'-: (b) Triangles after diagonal swapping. Every time a point insertion or diagonal swap operation is performed on the surface triangulation, corresponding changes are made to the control mesh. Deletion of a triangle causes deletion of three cells in the control mesh and creation of new triangles leads to formation of new cells by local application of the surface spline algorithm discussed in the previous section. The control mesh is thus coupled dynamically with the surface triangulation. Mesh Update After an initial interface position and shape is defined, the interface moves on the basis of the velocity field. The normal at all the nodes is required for the movement of the interface. Curvature pl1,j, an important role in the estimation of interface characteristics as well as interfacial fluid dynamics. At the end of the mesh refinement step, we have a surface approximated by an unstructured mesh. We also have a C1 smooth surface composed of parametric Bernstein-Bezier patches obtained by the application of the triangular de Co-t.li.ii algorithm on the control mesh. The control mesh here is the outcome of the C' Surface Spline algorithm applied on the unstructured surface mesh. ,lr i IId of computing normals and curvatures on these surface representations are quite different because one representation is implicit while the other is parametric. The details of the methods used are given below. Mesh Update on C' Spline Surface The application of C' Spline Algorithm results in a surface expressed in the Bernstein-Bezier form. Each parametric patch of the surface is of the form: x = x(u, v) = (x(u, v), y(u, v), z(u, v)) (3.3) Here cartesian coordinates x, y, z of a surface are differentiable functions of the parameter u and v. For a Bezier triangle the parametric domain is a triangle and (u, v, 1 u v) are the barycentric coordinates of a point in the reference triangle. Numerical calculation of normal. The partial x, and x, at a point x span the tangent plane to the surface at x. The normal of the tangent plane, n, coincides with the normal to the surface at x and is given by X, X X, n = X (3.4) X, X X, The normal is of unit length and is perpendicular to x, and x,. The normal together with unnormalized vectors x, and x, form a local coordinate system. If n triangles abut at a particular node in the surface triangulation, 2n cubic Bezier patches meet at the corresponding node in the control mesh. A normal at a point P is computed by taking the normalized mean of normals obtained from all the patches meeting at P. Numerical calculation of curvature. Let n denote a unit vector normal to the surface S at a point P and consider the intersections, with S, of planes containing n. Each of these planes meets the surface in a curve of intersection that is referred to as a normal section at P. A direction can be associated to each normal section, the direction of the unit tangent vector, t, to the curve at P. Every normal section also has a normal surface curvature, the curvature, Kt, of the curve at P. These normal surface curvatures vary as the normal cutting plane rotates around n. The maximum and minimum values of it, which occur in orthogonal directions, are referred to as the principal curvatures at P. If i1 and K2 denote the principal curvatures, the mean curvature, H, is given by L 1 +K2 NE-2MF+LG 2 2(EG F2) Here E = xxu, F = xx,, G = x,x,, L = nx u, M = nx,,, N = nx,. Derivatives xu, x,, x U, xU,, x, can be obtained by application of the de C',.ti 1ili algorithm for Bezier triangles with respect to the directions of derivative and the position u = (u, v, 1 U v). Both t1 and r2 change sign if the normal, n, is reversed and H is affected by such a reversal. Mesh Update on Surface Triangulation Triangulated mesh yields data which is of implicit nature. The coordinates of the nodes of the mesh along with the connectivity information are available to us. The methods for computing normals and curvatures on such a surface are different from those applicable for a smooth parametric patch. This section is devoted to the ways of computing normals and curvatures on triangulated surfaces. Numerical calculation of normal. Normal at a point P is basically a vector which is perpendicular to the tangent plane at P. On a triangulated surface, planes containing triangles meeting at P can be considered as tangent planes at P. Due to the absence of a unique tangent plane, there is no unique way of computing a surface normal in this case. The method used is based on geometrical averaging of normals obtained from triangles meeting at every node of the mesh. Consider a triangle having nodes x_-, x0 and x+i. Define vectors vl = x_- xo and v2 = X+1 x0. The unit normal to the triangle is now simply V1 X V2 n =V (3.6) V1l X v.|| A normal at a point P is computed by normalizing the mean of normals obtained from all the triangles meeting at P. Surface curvature using the Dupin indicatrix. For an arbitrary interface, a method for calculating curvatures for a given distribution of discrete grid points on the interface has been proposed by Todd and 'I, I.,,,d [22]. This method uses the Dupin indicatrix from given surface point data and is significantly more accurate than those using local quadric interpolants for the structured grid cases considered by Todd and '.1 1, I.-d. The approach along with some modifications made for computing surface curvatures is discussed here. Let I1 and K2 denote the principal curvatures and tx, t2 denote the corresponding directions. Further, let T denote the angle between the direction t of an arbitrary normal section and the direction tx. Euler showed in 1760 [23, 24] that the normal curvature Kt in an arbitrary direction t is related to the principal curvatures r1 and K2 by Kt K= 1 cos2 t + K2 sin2 2 (3.7) If X = x cos /(l t )1/2 and y = sin/( K(t )1/2 then the Dupin indicatrix of the surface S at the point P is given by KlX2 a+ .2 ./- = 1 (3.8) If the principal curvatures are of same sign, the normal curvature in any direction is also of that sign and the surface, locally, lies entirely on one side of the tangent plane at P. In this case, Equation (3.8) provides a nontrivial curve which is an ellipse. When the principal curvatures are of differing sign, that is, when P is hyperbolic or saddle point, Dupin indicatrix is a pair of hyperbolas. If axes other than those in the direction of the principal curvatures are used, Dupin indicatrix takes the form: Ax2 + 2Cxy + By2 = 1 (3.9) Hence if Dupin indicatrix at a point is known, so are principal curvatures and their directions. The mean curvature can be computed by taking the average of the principal curvature values. In turn, knowing estimates of normal curvatures in particular directions, some approximate fit to that data by a central conic of the form of Equation (3.9) can be made and from there, via a rotation, the Dupin indicatrix can be estimated. From the coordinates of points on the surface, however, the normal curvatures cannot be estimated directly. To circumvent this difficulty, a result due to \li',-iili n, ,- be used. Let K be the curvature, at P, of a curve obtained by the intersection of an arbitrary plane with the surface S and let Kt be the normal curvature in the same direction. If is the angle between the arbitrary plane and the normal plane then, according to llii -i'i l's theorem [23, 24], we get the following result: Kt = K cos O (3.10) The surface curvature at any point P can be computed by estimating the planar curvatures in plane sections through P. Three distinct points in a three-dimensional space determine a plane. Within this plane these three points determine a circle and hence curvature and tangent vector along the circle. Let the three distinct points be x_-, xo and x+l. Performing a translation so that xo is now at the origin, it is straightforward to obtain the circumcenter, x, = (a, b, c), of the triangle formed by the three points. Knowing the circumcenter, the curvature vector is given by (a, b, c) (a2 + b2 + C2) Defining us define vectors u and v as u = x_1 xo and v = x+1 xo. The unit tangent vector to the circle passing through x_-, xo and x+l, at the location xo, is given by (u X v) X (3.12) t = (3.12) (u x v) x ,. 1 Given sets of three points, we can find the tangent and curvature vectors of the circles passing through each of these sets. Hence for any point on the surface, two neighboring points can be selected and from these three points an estimation of the curvature of the planar section through the initial point can be obtained. As far as implementation is concerned, xo and pairs of points lying on opposite sides of xo are chosen. For each of these sets of three points, let ti and K, be the corresponding tangent and planar curvature vectors as indicated before. Each tangent vector ti is tangent at the point xo to the circle defined by the point set under consideration. Hence each pair of tangent vectors, ti and tj, defines a plane that is tangent to two circles, each of which passes through the point xo. This plane is an approximation of the tangent plane to the surface at xo. Since the unit normal of the surface S at point P is known, \Ii'-inilrs theorem can be applied to give (3.13) Kn,i -= "i n Here Kn,i is an estimate of the normal curvature in the direction ti. ('liii i.IIg a coordinate frame such that the x y plane is perpendicular to the vector n, define t- x = (i, y) = )1/2(3.14) In the above equation, i = 1,..., n where n is the number of triangles meeting at point P. If Ax2 + 2Cxy + By2 = 1 is the Dupin indicatrix then in the absence of numerical error we get Ax2 + 2Cxzy, + By2 = sgn(K,,i) (3.15) Here again i = 1,...,n. The above given equation has three unknowns (A, B and C) and if n > 3, normal curvatures with their associated errors cannot all be expected to lie on a conic of this form. In this case, a least-squares approximation to the Dupin indicatrix can be made by minimizing the following function: F = (Ai + 2Czxiy + By sgn(,n))2 (3.16) i=1 The minimum is given by values of A, B and C which satisfy the linear system of equations: Zx3i xYi zxy A xisgn(Kn,i) EXi E>' :l E:.' 2C = Ei i sgn(,i) (3.17) E': E xy3 E yf B E y sgn(Kr,i) All summations are over i where i varies from 1 to n. Now, the principal curvatures KI, K2 and the principal directions of curvature are eigenvalues and eigenvectors of the following matrix: A C C B Kr and K2 can be obtained by performing a rotation in the plane by an angle where tan = 2C/(A B). The Dupin indicatrix has a geometric interpretation as illustrated in Figure 3.11. A surface at x can be approximated by a paraboloid, that is, a Taylor expansion with terms up to second order. This leads to a simple interpretation of the Dupin indicatrix: the indicatrix, scaled down by a factor of m, can be viewed as the intersection of the surface with planes parallel to the tangent plane at x in the distance e = 1/(2m2). The Dupin indicatrix can be assigned a sign depending on whether it is ',,ilve" or "b, -.i. the tangent plane. Figure 3.11: Dupin indicatrix (scale 1 : m) obtained by the intersection of the surface with the lower plane. With all this information, different methods tried for numerically computing the mean curvature value are presented. All these methods, which are largely based on the Todd and '.1 I .,L-d approach, are discussed below: 1. Mean of Normal Curvatures: According to this method, surface curvature is estimated by computing the mean of the normal curvatures, {Kn,i}, obtained by choosing sets of three points. This is a quick way of computing curvatures and the accuracy of values obtained will depend on the number of normal curvature values used and directions along which normal curvatures are calculated. 2. Todd and McLeod Method: In this method, the mean curvature value is obtained by taking the average of the maximum and minimum normal curvature of a conic whose curvature, K, vs angle, T, curve "( I, -l,: fits the set of normal curvature values, {K,,i}, obtained from the local geometry around point P. The details of obtaining a conic by method of least squares and principal curvature values by matrix rotation have already been discussed. 3. Tangent Plane Intersection Method: This method is based on fitting a conic at P through points obtained by the intersection of two planes parallel to the tangent plane at P and the local surface around that point. The two intersecting planes are at equal distance on both sides of tangent plane. This distance is chosen such that one of the planes passes through the midpoint of the edge having shortest z component in the coordinate system where the x y plane is perpendicular to the normal vector, n, at P. Once the intersection points are available, an estimate of Dupin indicatrix is made by following the least squares method discussed before. The sign of curvature now depends on the intersection plane used to get the intersection point. Rest of the procedure for obtaining mean curvature from the Dupin indicatrix is followed as outlined earlier. CHAPTER 4 I:L>1I LTS Several examples showing the potentialities of the code developed are given below. Complex interface representation capability is demonstrated first. Grid refinement ability, where points are inserted on a stationary sphere based on the area criteria, is shown next. Various methods of computing curvatures are compared later by expanding a sphere, uniformly in its normal direction, without any refinement. The effects of grid adaption on normal and curvature calculations of a sphere and an ellipsoid are also shown. These tests demonstrate the capabilities of the method as an important tool in interface representation. Complex Feature Representation The ability of the tracker to represent complex shapes is demonstrated in this section. In Figure 4.1, velocity functions are prescribed to a sphere that result in the interface forming complex features. Figure 4.1a shows the initial sphere used to generate rest of the geometries. Figure 4.1b is a uniformly expanded sphere obtained by the velocity function: vn = 1 where v, is the velocity in the direction of the surface normal. An ellipsoid (Figure 4.1c) is obtained by the following velocity function: v., = x (velocity is positive when the x coordinate is positive and negative when the x coordinate is negative), vy = 0, vz = 0. Finally, Figure 4.1(d) shows the final shape obtained from the initial sphere by the velocity function: vn = sin( xl)sin(yl )sin( z|). In all these cases, the features of interface are clearly captured by unstructured surface grids used. It is also clear that the features could not be represented as well if structured grids were used for the initial surface. It would be difficult to refine and insert new points to a structured surface as compared to an unstructured one. .1! refinement is based on the triangle quality and is done according to the mesh refinement algorithm discussed earlier in the report. We are associating quality values with every node of the surface mesh in order to investigate the effect of mesh quality on the geometric features such as normals and curvatures. The mesh quality value is set in two ways: first, by computing the ratio of the areas of the largest and the smallest triangles meeting at a node (area ratio) and second, by computing the mean of the aspect ratios of all the triangles meeting at a node (average aspect ratio). Here, aspect ratio of a triangle is the ratio of its circumradius and inradius which is equal to 2 for an equilateral triangle. Figure 4.2 shows the area ratio contours of the shapes discussed before. The original sphere has triangles of unequal areas which leads to nonhomogenous point insertion. Considering this fact, the results are fairly good. The average aspect ratio contours are depicted in Figure 4.3. It is observed that, even after considerable deformations, the average aspect ratio values do not change significantly. This clearly demonstrates the mesh refinement capability of the interface tracker. The use of Bowyer-Watson algorithm for point insertion and circumcircle criterion for diagonal swapping results in a good quality mesh. The mesh quality values, especially the average aspect ratio values, have an effect on the numerical curvature calculations performed. This correlation is highlighted in the rest of this chapter. S(a) b) S (C) d) Figure 4.1: Grid refinement as a sphere expands into different shapes: (a) Initial sphere (402 points); (b) Fili,,] sphere (4,926 points); (c) Ellipsoid (3,794 points); (d) Fili,,] shape (2,454 points). ; ::r: ;' i, ::z . ,a) Area Ratio 1.244 1.212 S1.179 S1.146 S1.114 1.081 1.049 1.016 I -v (C) Area Ratio 2.811 2.570 2.328 2.087 1.845 1.604 1.362 1.121 (b) Area Ratio 2.637 2.419 2.201 1.983 1.765 1.546 1.328 1.110 (d) Area Ratio 2.538 2.333 2.128 1.923 1.718 1.513 1.308 1.103 Figure 4.2: ,l;:]] : i to minimum area ratio contours as a sphere expands into different shapes: (a) Initial sphere; (b) Fin,,] sphere; (c) Ellipsoid; (d) Fin,,] shape. *^ a) (b) Avg AR ,- Avg AR 2.379 v2.908 2.329 1. 2.789 2.279 2.670 2.229 2.552 S2.179 2.433 2.079 2.196 2.030 2.077 S (C)d) Avg AR Avg AR 3.099 2.920 S2.953 2.799 2.807 2.679 2.661 2.558 2.515 2.438 2.370 2.317 2.224 2.197 2.078 2.076 Figure 4.3: Average aspect ratio contours as a sphere expands into different shapes: (a) Initial sphere; (b) Fiii,1] sphere; (c) Illip,-,i Ii (d) Fiii,1] shape. T,i1.1, 4.1: E iT- r of area factor on the I. N s radius and curvature errors as mesh refinement takes place on a stationary sphere. Area factor \ i:ii Curvature Radius of nodes ('.. I.'.s) error (. I.'.s) error 0.2 3,262 1. -. 6.940 x 10-3 0.4 1,502 0.996 7.067 x 10-3 0.6 1,030 1.002 7.306 x 10-3 0.8 770 2.561 ,. x 10-3 1.0 658 2.259 2.259 x 10-3 Grid Refinement In order to test the accuracy of the point-placement strategy and the curvature values obtained, consider a stationary sphere that is refined on the basis of the area criteria. Therefore, points are inserted on a sphere that remains fixed in space. Table 4.1 shows the effect of the area factor (a triangle is refined if its area is greater than the average triangle area times an area factor) on the root mean square (I', I.N) radius and curvature errors of a sphere. The initial radius and the theoretical value of curvature are used to normalize the errors (expressed in percentage). Initially, that is, before any mesh refinement takes place, we have the following: Curvature ('. I s) error = 0..~, and Radius (.. I:?)S) error = 3.184 x 10-5. Curvature calculations are done using the Todd and .1, 1.,Ld method and the initial error in curvature, even though points accurately lie on a sphere, is because of the numerical computation of curvature from available discrete node data. Area factor does not have appreciable effect on the radius (.. I :'.:S) error as this error increases marginally with the increase in the number of nodes. This information is an indication of the fact that point placement is highly accurate. As far as curvature ('. I:.'s;) error is concerned, it shows erratic behavior. Since the curvature error not only depends on the location of points but also on the way curvature is computed numerically, we have attempted to correlate the curvature errors with mesh quality measures discussed before. We have also compared the configuration of the sphere that gives accurate curvatures (area factor 0.4) with the one that results in relatively large curvature errors (area factor 0.8). Figure 4.4 shows the initial configuration of a sphere about to be refined. Figure 4.4a shows the surface mesh composed of eight spherical patches. Here, curvature is computed using the Todd and '. 1 .ll d approach. Normalized curvature (I: -' ;) error contour is shown in Figure 4.4b. It is seen that there is some error in the zones where the surface patches meet. Figure 4.4c shows the normalized radius (I:'. ;) error contour. Radius error is very low in this case. The area ratio contour is shown in Figure 4.4d and the average aspect ratio contour is illustrated in Figure 4.4e. It is strikingly clear that there is a direct correlation between the curvature error and the aspect ratio. Even though points of the surface mesh lie on a sphere, as is evident from the negligible radius error, we still get appreciable curvature error. Hence, besides the new surface-mesh nodes being placed accurately, surface-mesh triangles should have good aspect ratios as well. This is the conclusion we can draw for getting accurate curvature results from the discrete node and their connectivity information available to us. (b) K Error 0.057 0.047 0.036 0.026 0.015 0.005 -0.005 -0.016 (d) Area Ratio 1.244 1.212 1 .179 / 1.146 I 1.114 1.081 1.049 1.016 ee) Avg AR 2.979 2.852 2.726 2.599 2.472 2.345 2.218 2.091 Figure 4.4: Initial configuration before mesh refinement takes place on a stationary sphere: (a) Initial surface mesh; (b) Contour of normalized curvature (It ')l) error; (c) Contour of normalized radius (I' I N ) error; (d) Contour of area ratio; (e) Contour of average aspect ratio. (a) (c) 9- R Error 5.7E-07 4.0E-07 2.4E-07 6.9E-08 -9.7E-08 -2.6E-07 -4.3E-07 -6.0E-07 Figure 4.5 shows the final configuration of a stationary sphere that is refined with an area factor = 0.8. The initial uniform unstructured mesh loses its symmetry due to the insertion of a few points. The number of new points added is less than the number of points in the mesh originally. This results in a nonuniform triangulation as shown in Figure 4.5a. This fact is also evident from the area ratio (Figure 4.5d) and aspect ratio (Figure 4.5e) contours drawn on a stationary sphere that is refined by point insertion and diagonal swap operations. Here again, the normalized curvature error (Figure 4.5b) is high in a zone of triangles having high aspect ratios. The case of a stationary sphere that is refined with an area factor = 0.4 (Figure 4.6) is quite different from the one discussed above. Since the maximum allowable area of a triangle is less in this case, there are more number of points inserted into the original surface mesh. The final configuration has more than double the number of points in the original triangulation, hence, every triangle in the initial unstructured surface mesh gets bisected resulting in a more uniform-looking triangulation (Figure 4.6a). This is also deducible from the area ratio (Figure 4.6d) and aspect ratio (Figure 4.6e) contours. Since the resulting triangulation is uniform, curvature errors (Figure 4.6b) are low. The conclusion drawn from the cases of stationary grid refinement is that the mean curvature of the surface not only depends on the accuracy of the placement of the points on the surface but also on the quality of the surface mesh. The quality of the unstructured mesh affects the mean curvature value because we do not have a local analytical representation of the surface and curvature is being computed numerically from discrete point and their connectivity information available to us. _Cb) K Error 0.057 0.047 0.037 0.028 0.018 0.009 -0.001 -0.011 (d) / Area Ratio 1.921 1.798 1.675 O/ 1.553 1.430 1.307 1.184 1.061 ee) Avg AR 2.979 2.852 2.726 2.599 2.472 2.345 2.218 2.091 Figure 4.5: Fi,,] configuration of a stationary sphere that is refined with an area factor = 0.8: (a) Finl,] surface mesh; (b) Contour of normalized curvature (It'.l%) error; (c) Contour of normalized radius (Iti\l) error; (d) Contour of area ratio; (e) Contour of average aspect ratio. (a) (c) R Error 8.9E-05 2.9E-06 -8.4E-05 -1.7E-04 -2.6E-04 -3.4E-04 -4.3E-04 -5.2E-04 \ I ~g\s (a) A1 A)24r (c) R Error 8.9E-05 2.9E-06 -8.4E-05 -1.7E-04 -2.6E-04 -3.4E-04 -4.3E-04 -5.2E-04 (e) Avg AR 2.979 2.852 2.726 2.599 2.472 2.345 2.218 2.091 Figure 4.6: Fin,,] configuration of a stationary sphere that is refined with an area factor = 0.4: (a) Finl,] surface mesh; (b) Contour of normalized curvature (It'.l%) error; (c) Contour of normalized radius (Iti\l) error; (d) Contour of area ratio; (e) Contour of average aspect ratio. _Cb) K Error 0.057 0.047 0.036 0.026 0.015 0.005 -0.005 -0.016 d) Area Ratio 2.735 2.510 2.286 2.061 1.836 1.612 1.387 1.163 *Ue Various methods adopted for evaluating the mean curvature are compared by uniformly expanding a unit sphere. Nodes on the surface are simply advected and there is no mesh refinement taking place in the cases discussed next. No Grid Refinement Figure 4.7 shows the curvature ('.. I:MlS) errors plotted against the number of time steps as a sphere uniformly expands in time. The curvature calculated over the parametric surface, obtained after application of the C' surface spline algorithm, is fairly accurate. The surface in this case is composed of Bernstein-BeMier triangular patches. Control points of these patches are used to compute derivatives that are required in curvature calculations. Blend ratios of values 0.83 each lead to this result. The error, however, increases with a different choice of ratios. Actually, the blend ratios determine the depth of cuts into the polytope outlined by the control mesh. Smaller cuts result in a surface that follows the input mesh more closely and changes the normal direction more rapidly across the boundary. The choice of the best blend ratios, hence, becomes an issue and the inability to determine them automatically makes the C' surface spline algorithm unuseable. l,.iii, of normal curvatures and Todd and '. 1, l,-d methods yield more accurate results. The curvature-error plots obtained from these two methods are almost identical in the absence of any mesh refinement. The curvature error grows in time with the application of the tangent plane intersection method and this method turns out to be highly inaccurate. The results obtained show that either the mean of normal curvatures or the Todd and '.1, 1.,-)d approach should be used for computing mean curvatures on arbitrary interfaces. 30 Plane Cutting 25 201- 15 10- C1 Surface Spline Todd and McLeod, Mean of Normal Curvatures I I I I I )10 20 30 40 5( Time Step Figure 4.7: Comparison of different methods used for computing curvatures on a sphere when it expands without any mesh refinement. Table 4.2: EfiT- i of time step size on the radius and curvature (.I I. 's) errors of a sphere expanding in the absence of any mesh refinement. Time \ii iil.-i of CFL Curvature Radius increment time steps (. I ,S1) error (' IMS,1) error 0.025 40 0.1 0.264 2 ,s x 10-3 0.125 8 0.5 0.264 2.760 x 10-3 0.250 4 1.0 0.265 2.638 x 10-3 1.000 1 4.0 0.263 2.131 x 10-3 Table 4.2 shows the effect of time step size (CFL) on curvature and radius errors. Curvature, in this case, is computed using the Todd and .1, I. L-,d approach which is based on the numerical approximation of the Dupin indicatrix. The final configuration of the sphere is same as it expands for unit time in all the cases. From the results obtained, it is clear that time step size does not have any significant effect in the absence of mesh refinement. The radius error is negligible and the curvature error is also insignificantly same in all these cases. These results also indicate that the normals, used to advect the nodes of the mesh, are computed accurately. Grid Adaptation Figure 4.8 shows the contour plots of components of normals of various shapes. The contours of z components of normals of the initial sphere (Figure 4.8a) and the final sphere (Figure 4.8b) are circles as expected. Figure 4.8c shows the contours of y components of normals of an ellipsoid and Figure 4.8d depicts the contours of z components of normals computed numerically on the final shape. I(a) (b) NZ NZ 0.875 0.875 S0.625 0.625 0.375 0.375 0.125 0.125 -0.125 -0.125 -0.375 -0.375 -0.625 -0.625 -0.875 -0.875 S(C)(d) N.Y N.Z 0.875 0.875 0.625 0.625 0.375 0.375 0.125 0.125 -0.125 -0.125 -0.375 -0.375 -0.625 -0.625 -0.875 -0.875 Figure 4.8: Contours of components of normals obtained after mesh refinement: (a) z component of normal of initial sphere; (b) z component of normal of final sphere; (c) y component of normal of ellipsoid; (d) z component of normal of final shape. Table 4.3: EI~Tf i of time step size on the radius and curvature sphere that gets refined during uniform expansion. (.I I:'.1s) errors of a CFL N\ illii,,.i of nodes 0.1 1,5s- 0.5 1,609 1.0 1,613 4.0 1,592 Curvature (, It',IS) error 0.999 0.929 0 I 12 1.416 Radius (.. I.l;) error 7.658 x 10-3 7.365 x 10-3 7.311 x 10-3 7.440 x 10-3 Table 4.3 shows the effect of time step size (CFL) on the curvature and radius errors of a uniformly expanding sphere. The final configuration of the sphere is the same in all these cases. Area factor = 1.5 here. Curvature is computed using the Todd and '.I, I.-,d method. The curvature ItP'sl errors show an erratic trend but they approximately have the same value. The randomness in the result is attributed to the quality of the mesh obtained after all these cases. Radius I t'; errors, on the other hand, are hardly affected by the time step size. Figure 4.9 shows the final configuration of a sphere (initial radius = 1 and area factor = 1.5), that gets refined during uniform expansion, after 50 time steps (Time step size = 0.05). The resulting surface mesh is shown in Figure 4.9a. The new point placement is accurate as is evident from the radius error (Figure 4.9c) contour. Area ratio contours of the surface mesh are depicted in Figure 4.9d. A careful observation shows that high curvature error (Figure 4.9b) occurs in the regions where triangles have high aspect ratios (Figure 4.9e). Time increment 0.025 0.125 0.250 1.000 1 1 ili i-' of time steps 40 8 4 1 73 (a) (b) K Error 0.049 S0.038 0.026 0.014 0.002 -0.009 -0.021 -0.033 c) d) Srror Area Ratio / A t RError ., , I -1.6E-05 2.419 1.3E--5 2.201 -4.5E-05 2.201 S' -7.4E-05 1.983 -1.0E-04 1.765 S- -1.3E-04 '1.546 -1.6E-04 1.328 -1.9E-04 1.110 X 1.110 ,,-. Avg AR S ," 2908 1Z 2.789 2.670 .z 2.552 2.314 2.196 2.077 Figure 4.9: Finl] configuration of a sphere (initial radius = 1), that gets refined during uniform expansion, after 50 time steps (dt = 0.05): (a) Finil surface mesh; (b) Contour of normalized curvature (It: ) error; (c) Contour of normalized radius (It'.ls) error; (d) Contour of area ratio; (e) Contour of average aspect ratio. Table 4.4: E IT-, i of area factor on the radius and curvature ('.. I.'.si) errors as mesh refinement takes place on a sphere expanding uniformly in the outward direction. Area factor \nili.i CFL Curvature Radius of nodes (. I ,1'\1) error ('.. I,1S) error 0.5 4,733 0.37 1.235 4.179 x 10-3 1.0 2,501 0.26 1.364 6.191 x 10-3 1.5 1,586 0.22 0.866 7.440 x 10-3 2.0 1,242 0.19 0.971 8.459 x 10-3 2.5 962 0.17 1.247 8.970 x 10-3 The effect of the area factor, on the curvature and radius I '. errors, is shown in Table 4.4. The case is that of a uniformly expanding sphere that is refined as it grows in size. The time step size = 0.05 and the number of time steps = 50. Curvature calculations are done using the Todd and I'. I.-,d approach. Curvature I '.1%- error remains almost the same. The variations in these values are because of the quality of the mesh obtained at the end. The radius error increases slightly with an increase in the area factor but this error is insignificant because of the accurate placement of the new points. A comparison between the Todd and '.1 I .,ld method and the mean of normal curvature method is made below. For a uniformly expanding sphere (Figure 4.10), the mean of normal curvatures computed along a few directions at every node results in more accurate results than those obtained by the Todd and i'.1 I.,-,d approach. The initial radius of the sphere = 1, area factor = 1.5, time step size = 0.05 and the number of time steps = 50. The results in both the cases depend on the choice of points used for computing normal curvatures. It is observed that a pair of points lying opposite to each other with the given point in the middle form a set of three points which results in good normal curvature values. The accuracy decreases if the pair of points lie on the same side of the concerned node. For a sphere, the graphs of Todd and '.1 I.,- )d curvature error and mean of normal curvature error have the same trend. In the latter case, the curvature error peaks are, however, lower and this is because of the averaging involved which smoothens out the curvature error. Figure 4.11 shows the variation of normalized curvature ('.. I's.) errors with time steps for an ellipsoid. The ellipsoid is obtained from a sphere (initial radius = 1) expanding with a time step size of 0.05. The area factor = 1.5 here. Theoretically, the x coordinate of a point on an ellipsoid can be expressed in terms of its y and z coordinates as: x = V/1 -y2 (1 + At)n. Here At is the time step size and n is the number of time steps. This relation is used to compute the exact curvature and therefore curvature error of an ellipsoid. The graph shows that the Todd and S1, I.,, d method works better for an ellipsoid where normal curvatures, at arbitrary nodes on the surface, change with the tangential directions. The mean of normal curvature is not as accurate because in principle, the mean of the normal curvatures obtained at a few normal directions ino," not be equal to the average of the maximum and minimum normal curvatures at the concerned point. The error depends on the directions chosen for computing normal curvatures. As a conclusion, we can say that the Todd and '. 1, dd method should be used for numerically computing curvature over arbitrary moving interfaces. 3 2.5 2 W 2 S1.5 T Todd SM ean 0.5 0 10 20 30 40 50 Time Step Figure 4.10: A comparison between the Todd and 1I, l-,,d approach and the mean of normal curvatures method used for numerically computing curvature of a uniformly expanding sphere. 16 15 14 13 S- Mean 12 w 11- w 10 - 9 S 8- S. Todd ) 7 S6- 0 4- 3 - 2 1- 0 10 20 30 40 50 Time Step Figure 4.11: A comparison between the Todd and I. -nd approach and the mean of normal curvatures method used for numerically computing curvature of an ellipsoid. Figure 4.12 shows the final configuration of an ellipsoid, obtained from a sphere of unit radius, after 20 time steps (time step size = 0.05). As the ellipsoid grows in time, it gets refined by the point insertion and diagonal swap operations. Point insertion, performed according to the local sphere approximation (Subsection 3), introduces an error in the radial distance of an arbitrary point on an ellipsoid from the origin (we call this a i',, li. error"). This is shown in Figure 4.12c. The resulting triangulation has good area ratios (Figure 4.12d) and aspect ratios (Figure 4.12e) though. The curvature error of an ellipsoid (Figure 4.12b) appears to be related to the radius errors. A close observation shows that regions of high radius error correspond to zones of high curvature error. This behavior is different from that observed for a sphere where the main contribution to the curvature error came from bad aspect ratios. We have tried to circumvent the problem arising due to radius errors by projecting points exactly on an ellipsoid. The mesh is refined at every time step by insertion of points according to the local sphere approximation and the new points are exactly placed on an ellipsoid by listingg their x coordinate according to their y and z coordinates. Figure 4.13 shows the final configuration of such an ellipsoid, obtained by the expansion of a sphere (radius = 1), after 20 time steps (time step size = 0.05). From a visual inspection of the final surface mesh (Figure 4.13a), we observe that the ellipsoid has high area ratios (Figure 4.13d) and aspect ratios (Figure 4.13e) towards its center. The curvature errors (Figure 4.13b) here seem to be directly related to the aspect ratios. K Error 0.082 I t 0.040 S -0.003 -0.045 -0.087 -0.130 -0.172 -0.214 R Error 0.015 0.013 0.011 I 0.008 0.006 0.004 0.002 I-0.001 - . Area Ratio 2.313 2.138 1.963 1.788 1.613 1.438 1.263 1.088 Avg AR 2.861 2.747 2.633 2.518 2.404 2.175 2.061 Figure 4.12: Fill configuration of an ellipsoid, obtained from a sphere (radius = 1), after 20 time steps (dt = 0.05): (a) Fil, surface mesh; (b) Contour of normalized curvature (I: \,) error; (c) Contour of normalized radius (I: \,1) error; (d) Contour of area ratio; (e) Contour of average aspect ratio. S. 9- SI R Error 0.015 0.013 0.011 0.008 0.006 0.004 0.002 -0.001 K Error 0.082 0.040 / -0.003 -0.045 -0.087 -0.130 -0.172 -0.214 4 , Area Ratio 2.313 2.138 1.963 1.788 1.613 1.438 1.263 1.088 Avg AR 2.861 2.747 2.633 2.518 2.404 2.175 2.061 Figure 4.13: Fill configuration of an ellipsoid, obtained from a sphere (radius = 1) such that there is no radius error, after 20 time steps (dt = 0.05): (a) Fiil, surface mesh; (b) Contour of normalized curvature (Itl% ) error; (c) Contour of normalized radius (It: IH) error; (d) Contour of area ratio; (e) Contour of average aspect ratio. .. -_ ..---_- .-- _- 1,f,- i ,-r 4A CHAPTER 5 CONCLUSIONS AND FUTURE WORK The algorithm for tracking arbitrary interfaces in three dimensions is developed in object oriented C++. Some of the difficulties associated with unstructured grids are related to the numerous bookkeeping operations that frequently need to be performed during node, edge and triangle additions, reconnections and subtractions. It becomes important to design appropriate data structures that optimize both search and storage procedures. Our code takes full advantage of useful data structures such as stacks and linked lists and makes extensive use of the Standard Template Library (STL) for performing operations on them. The decision of using unstructured grids for representing a moving interface turns out to be a judicious one. Grid refinement, which is possible in this case, is very useful in preserving complex interface features as they develop. Point insertion by the application of the Bowyer-Watson algorithm and diagonal swapping by using the circumcircle criterion results in D.1,,iiii,,' triangulation having fairly good aspect ratios and area ratios. Calculation of normals and surface curvatures is an important part of the interface- tracking code. Computation of the surface normal by taking the mean of normals of triangles meeting at the concerned point results in accurate values. As far as curvature calculations are concerned, the Todd and ., I1 ld and the mean of normal curvatures methods give highly accurate results when compared with the values obtained by the plane cutting and C' surface spline methods. For an arbitrary interface, where the normal curvature value varies with orientation of the normal plane, the numerical approximation to the Dupin indicatrix by the Todd and '.. 1, lnd approach is found to be more accurate than the mean of normal curvature values obtained along a few directions. Curvature values are related to the quality of the mesh and significant errors are obtained at locations where triangles having bad aspect ratios meet. Incorrect placement of the new point during the point insertion operation also contributes to this error. The local sphere approximation, used to calculate the coordinates of these new points, leads to placement of the point on a symmetric spherical patch having uniform curvature values along all directions. This technique, therefore, is crude because, in reality, for an arbitrary surface the curvature at the new node location varies with the orientation of the normal plane. The efforts in the future, thus, should be directed towards insertion of the point at an appropriate location. We need to place the new point on a smooth parametric surface patch whose control points, unlike the C' spline control points, coincide with the points of the unstructured triangulated surface mesh. APPENDIX A POINTS AND VECTORS Points are elements of three-dimensional euclidean (or point) space E3. They identify a location, often relative to other objects. Examples are the midpoint of a straight line segment or the center of gravity of a pl-'ii ,'1 object. Vectors are elements of three-dimensional linear (or vector) space R3. They have a magnitude and a direction associated with them. Although both points and vectors are described by triples of real numbers, there is a clear distinction between them. For any two points a and b, there is a unique vector v that points from a to b. This is produced by componentwise subtraction v = b a. On the other hand, given a vector v, there are infinitely many pairs of points a, b such that v = b a. For if a, b is one such pair and if w is an arbitrary vector, then, a + w, b + w is another such pair since v = (b + w) (a + w). Assigning the point a+w to every point a E 3 is called a translation and vectors are invariant under translations while points are not. Elements of point space E3 can only be subtracted from each other and this operation yields a vector. They cannot be added as this operation is not defined for points (it is defined for vectors). However, addition like operations are defined for points and they are called barycentric combinations 1. These are weighted sum of points where the weights sum to one: 1Also called F. ~.. combinations. b = Yj-o ajbj where ,,, + .- + a, = 1. Barycentric combinations generate points from points. An important special case of barycentric combinations are convex combinations. These are barycentric combinations where the coefficients aj, in addition to summing to one, are also nonnegative. A convex combination of points is always I II I" those points. This observation leads to the definition of the convex hull of a point set which is the set formed by all convex combinations of a given point set. Figure A.1 gives an example of a convex hull. Figure A.1: Convex hull of a point set (a polygon) shown shaded. The convex hull of a point set is a convex set. Such a set is characterized by the following: for any two points in the set, the straight line connecting them is also contained in the set. APPENDIX B BARYCENTRIC COORDINATES Consider a triangle with vertices a, b, c and a fourth point p, all in E2. It is always possible to write p as a barycentic combination of a, b and c as p = ua + vb + we where u+v+w = 1. The coefficients, u = (u, v, w), are called barycentric coordinates of p with respect to a, b and c. If four points a, b, c and p are given, we can find the barycentric coordinates of p by Cramer's rule: area(p, b, c) area(a, p, c) area(a, b, p).1) area(a, b, c)' area(a, b, c)' area(a, b, c) The Cramer's rule makes use of determinants and they are related to the areas by the identity: ax bx cx area(a, b, c) = a by c (B.2) 1 1 1 APPENDIX C DATA STRUCTURES FOR INTERFACE TRACKING Node Data Structure * Index a unique number for every node. * Coordinates x, y and z coordinates of the node. * Curvature value of mean curvature at the node location. * Normal a unit vector normal to the surface at the node location. * Velocity a vector for the velocity of the node. * EdgeList list of edges meeting at the node. * TriangleList list of triangles meeting at the node. Edge Data Structure * Index a unique number for every edge. * Nodl,-'2] nodes of this edge. * Triangle[2] triangles sharing this edge. * El -.i r [i] lists of edges storing the other edges (in order) of the two triangles of this edge. 87 * Direction 2] values to indicate whether the direction of motion in the two triangles is correct or reverse as we move from the first to the second node of this edge. * State a number that indicates whether the edge is ";iI]'i., "dt.iiiiiit" or "di-,dl depending on the number of triangles (two, one or zero respectively) this edge is part of. Triangle Data Structure * Index a unique number for every triangle. * Area area of this triangle. * Fi t EdI an edge that is a handle to the nodes and other edges of this triangle. Surface Mesh Data Structure * niii i iOfNodes total number of nodes in the surface mesh. * nil iiOfEdges total number of edges in the surface mesh. * \iili i.i OfTriangles total number of triangles in the surface mesh. * NodeList list of all the nodes in the surface mesh. * EdgeList list of all the edges in the surface mesh. * TriangleList list of all the triangles in the surface mesh. * \, l ::i]] iiiiArea area criteria value beyond which a triangle is bisected into smaller triangles. * BadEdgeList list of edges to be swapped through the diagonal swap operation. * BigTriangleList list of triangles undergoing the Bowyer-Watson point insertion algorithm. REFERENCES [1] W. S. Shyy, H. S. Udaykumar, M. M. Rao and R. W. Smith, Computational fluid dynamics with moving boundaries, Taylor and Francis, Washington, DC (1996). [2] J. M. Floryan and H. Rasmussen, \n iliil ,, i methods for viscous flows with moving boundaries, Applied Mechanics Review 42, 323-341 (1989). [3] S. Osher and J. Sethian, Fronts propagating with curvature dependent speed: Algorithms based on Hamilton-Jacobi formulations, Journal of Computational Physics 79, 12-49 (i'-s). [4] J. A. Sethian, Level set methods: Evolving interfaces in geometry, fluid me- chanics, computer vision and materials science, Cambridge University Press, Cambridge (1996). [5] C. W. Hirt and B. D. \?i i, i, Volume of fluid (VOF) method for the dynamics of free boundaries, Journal of Computational Physics 39, 201-225 (1981). [6] F. H. Harlow and J. E. Welch, \Niln,. II ,iI calculation of time-dependent viscous incompressible flow of fluid with free surface, Physics of Fluids 8(12), 2182-2189 (1 ,7,). [7] S. O. Unverdi and G. Tryggvason, A front tracking method for viscous, incom- pressible, multi-fluid flows, Journal of Computational Physics 100, 25-37 (1992). [8] H. S. Udaykumar, W. S. Shyy and M. M. Rao, ELAFINT: A mixed Eulerian- Lagrangian method for fluid flows with complex and moving boundaries, Inter- national Journal of Nii i;.-, ,,1 Methods in Fluids 22(8), 691-712 (1996). [9] H. S. Udaykumar and W. S. Shyy, Simulation of morphological instabilities dur- ing solidification; Part I: Conduction and capillary effects, International Journal of Heat and Mass Transfer 38, 2057-2073 (1995). [10] C. S. Peskin, N\il:i.iI ,li analysis of blood flow in the heart, Journal of Compu- tational Physics 25, 220-252 (1977). [11] P. L. George, Automatic mesh generation: Application to finite element methods, John Wiley & Sons (1991). [12] G. Farin, Curves and surfaces for computer aided geometric design, Academic Press, N\.-- York, 4th edition (1997). [13] L. Piegl and W. Tiller, TI,. NURBS book, Monographs in Visual Communica- tion, Springer Verlag, N\.-- York, 2nd edition (1997). [14] V. .1,i: .,i'iiiiii H. S. Udaykumar and W. S. Shyy, Adaptive unstructured grid for three-dimensional interface representation, Nu iirn i o, I Heat Transfer, Part B 32, 247-265 (1997). [15] J. Peters, C' surface splines, SIAM Journal of Nii, L- ;1, Analysis 32(2), 645-666 (1993). [16] W. Boehm, G. Farin and J. Kahmann, A survey of curve and surface methods in CAGD, Computer Aided Geometric Design 1, 1-60 (1984). [17] J. Peters, Constructing C' surfaces of arbitrary topology using biquadratic and bicubic splines, in Designing fair curves and surfaces, N. Sapidis (ed.), SIAM, Philadelphia (1994). [18] J. Peters, Smooth free-form surfaces over irregular meshes generalizing quadratic splines, Computer Aided Geometric Design 10, 347-361 (1993). [19] A. Bowyer, Computing Dirichlet tesselations, Th,- Computer Journal 24(2), 162- 166 (1981). [20] D. F. Watson, Computing the n-dimensional D,-1.]iljii, tesselation with applica- tion to Voronoi polytopes, Th,- Computer Journal 24(2), 167-172 (1981). [21] N. P. Weatherill, D,.1.]i,',: triangulation in computational fluid dynamics, Com- puters and Mathematics with Applications 24(5/6), 129-150 (1992). [22] P. H. Todd and R. J. Y. `i 1.l, d, \Niliiiilii ,,i estimation of the curvature of surfaces, Computer Aided Design Journal 18(1), 33-37 (1986). [23] D. J. Struik, Lectures on classical differential geometry, Addison-Wesley, Cam- bridge, MA (1950). [24] T. J. \\Vlliiii,1, An introduction to differential geometry, Oxford University Press, N.-. York (1959). BIOGRAPHICAL SKETCH Abhishek Khanna was born in 1976 in \-,.- Delhi, India, where he spent his entire childhood. He did his schooling from Don Bosco School, Delhi, and Delhi Public School, R. K. Puram, \I-.- Delhi. He received his bachelor's degree in mechanical engineering from the Indian Institute of Technology, B iilJ,, -. After graduation, he came to the United States of America to pursue a master of science degree in the field of engineering mechanics at the University of Florida. As part of his thesis work, he has been developing a three-dimensional tracker code used for tracking interfaces in moving boundary problems. |