Title: Three-dimensional representation of evolving interfaces in computational fluid dynamics
CITATION PDF VIEWER THUMBNAILS PAGE IMAGE ZOOMABLE
Full Citation
STANDARD VIEW MARC VIEW
Permanent Link: http://ufdc.ufl.edu/UF00100807/00001
 Material Information
Title: Three-dimensional representation of evolving interfaces in computational fluid dynamics
Physical Description: Book
Language: English
Creator: Khanna, Abhishek, 1976-
Publisher: State University System of Florida
Place of Publication: Florida
Florida
Publication Date: 2000
Copyright Date: 2000
 Subjects
Subject: Aerospace Engineering, Mechanics, and Engineering Science thesis, M.S   ( lcsh )
Dissertations, Academic -- Aerospace Engineering, Mechanics, and Engineering Science -- UF   ( lcsh )
Genre: government publication (state, provincial, terriorial, dependent)   ( marcgt )
bibliography   ( marcgt )
theses   ( marcgt )
non-fiction   ( marcgt )
 Notes
Summary: 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.
Summary: ABSTRACT (cont.): 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.
Statement of Responsibility: by Abhishek Khanna.
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.
 Record Information
Bibliographic ID: UF00100807
Volume ID: VID00001
Source Institution: University of Florida
Holding Location: University of Florida
Rights Management: All rights reserved by the source institution and holding location.
Resource Identifier: oclc - 47681989
alephbibnum - 002678729
notis - ANE5956

Downloads

This item has the following downloads:

MSReport ( PDF )


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.




University of Florida Home Page
© 2004 - 2010 University of Florida George A. Smathers Libraries.
All rights reserved.

Acceptable Use, Copyright, and Disclaimer Statement
Last updated October 10, 2010 - - mvs