|
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.
|
