Discrete Surface Evolution and Mesh Deformation for
Aircraft Icing Applications
David Thompson
∗
, Xiaoling Tong
†
, Qiuhan Arnoldus
‡
, Eric Collins
§
,
David McLaurin
¶
, and Edward Luke
k
Mississippi State University, Mississippi State, MS, 39762, USA
Colin Bidwell
∗∗
NASA Glenn Research Center, Cleveland, OH, 44315, USA
Robust, automated mesh generation for problems with deforming geometries, such as
ice accreting on aerodynamic surfaces, remains a challenging problem. Here we describe a
technique to deform a discrete surface as it evolves due to the accretion of ice. The surface
evolution algorithm is based on a smoothed, face-offseting approach. We also describe a fast
algebraic technique to propagate the computed surface deformations into the surrounding
volume mesh while maintaining geometric mesh quality. Preliminary results presented here
demonstrate the efficacy of the approach for a sphere with a prescribed accretion rate, a
rime ice accretion, and a more complex glaze ice accretion.
I. Introduction
There are several technical challenges associated with mesh generation for simulating ice accretion on a
three-dimensional configuration. First, ice accretion is an evolutionary process; therefore, the mesh must
evolve in response to the growth of the ice shape. Assuming a loosely coupled ice accretion strategy, such
as that used in LEWICE3D,
1
a sequence of quasi-static accretion steps is performed to generate the final
ice shape. Since the ice shape changes, each accretion step requires a new mesh. However, a full mesh
regeneration may be expensive for complex configurations. An alternative strategy is to deform the mesh
in response to the ice growth. The second challenge is that accreted ice shapes can become quite complex.
While the current state of the art in ice accretion prediction does not produce shapes with exceedingly
complicated geometries, the predicted ice shapes can nevertheless present significant challenges for meshing
software. Additionally, as the fidelity of ice accretion prediction increases, the complexity of the numerically-
generated ice shapes will increase. Any mesh generation strategy will necessarily require the ability to handle
such complex geometries if it is to represent a viable, long-term solution.
Currently, there is no automated mesh generation process designed to work with LEWICE3D. The re-
sulting capability gap precludes routine grid-based, multi-time-step simulations of ice accretion on complex
configurations. As part of the NASA Atmospheric Environment Safety Technology Project, an ongoing
effort at Mississippi State University seeks to facilitate routine simulation of ice accretion on realistic, three-
dimensional configurations by developing a suite of meshing tools that will produce unstructured, mixed
element (hybrid) meshes for evolving ice shapes in an automatic, efficient, and robust manner. Such auto-
mated mesh generation is a necessary step in the enhancement of existing ice accretion prediction tools as
well as in the development of the next generation of these tools.
∗
Associate Professor, Department of Aerospace Engineering, PO Box A, Associate Fellow.
†
Assistant Research Professor, Center for Advanced Vehicular Systems, PO Box 5405, Associate Member.
‡
Research Associate II, Center for Advanced Vehicular Systems, PO Box 5405, Nonmember.
§
Postdoctoral Associate, Center for Advanced Vehicular Systems, PO Box 5405, Member.
¶
Assistant Research Professor, Center for Advanced Vehicular Systems, PO Box 5405, Member.
k
Associate Professor, Department of Computer Science and Engineering, PO Box 9637, Senior Member.
∗∗
Aerospace Engineer, Icing Branch, 21000 Brookpark Road, Member.
1 of 20
American Institute of Aeronautics and Astronautics
Any approach that has the potential to make grid-based ice accretion simulations for realistic configura-
tions commonplace occurrences must have the following characteristics:
• Automation: Simulating the evolving ice shape necessarily requires generating a new mesh for each
ice shape. For this approach to be routine, it is necessary that the mesh generation process be as
automated as possible. Once an initial mesh is generated, the user should be removed from the loop,
even in cases where the surface evolution initially produces an invalid mesh or a mesh of poor quality
due to the complexity of the ice shape.
• Efficiency: Although a new mesh must be generated for each ice shape, the process must be efficient.
The simplest approach, completely regenerating the mesh, is potentially a time-consuming task for
complex aircraft configurations and not appropriate for ice accretion simulations.
• Robustness: Any mesh generation tool that is to be employed in an automated analysis environment
must be robust. In the context of mesh generation for ice accretion simulations, robustness implies
that a valid mesh of reasonable quality must be generated for ice shapes of varying complexity. The
challenge here is to ensure that the mesh retains sufficient quality as the ice surface evolves.
In this paper, we briefly describe the algorithms we employ to evolve the discrete surface mesh that
represents the accreting ice and to project these deformations into the volume mesh. We are developing a
solver-neutral interface propagation tool that computes the position of a discrete surface as it evolves under
an accretion rate map specified by LEWICE3D. The surface evolution tool, iceSurf, is designed to produce
a new surface mesh given the current surface mesh, a face-centroid accretion rate map, and the icing time.
The output from iceSurf provides input, in the form of a surface displacement file, to the mesh deformation
tool gridMover to produce a new volume mesh. We include preliminary results for a sphere with a prescribed
accretion rate field, a relatively simple rime ice accretion, and a more complex glaze ice accretion.
II. Background: Surface Evolution
Figure 1. Face offsetting produces unambiguous nodal po-
sitions in two dimensions.
One of the challenges associated with evolving a
faceted, discrete surface is that the normal at a node
is not unique. This is caused by the discontinuous
nature of the discrete representation of the surface.
One possible solution is to define a displacement di-
rection at each node, based on the normals in the
adjacent faces, and displace the surface a prescribed
distance in this direction at each node. However,
there are numerous challenges associated with this
approach not the least of which is conservation of
volume. Alternatively, the surface evolution could be modeled by generating a plane that is parallel to a
given face by extruding a specified distance – the product of the accretion rate and the time step – from the
face centroid in the direction of the face normal as shown for a two-dimensional surface in figure 1. As seen in
the figure, there is no ambiguity in the location of the nodes in the new surface in two dimensions; however,
this is not the case in three dimensions in which any two of these offset planes (not parallel), intersect in a
line while three non-parallel planes intersect at a point. In general, the intersection of four or more planes
is not defined in three dimensions. Except in special cases, the number of faces that share a given node in
a typical triangular surface mesh is usually more than three and, consequently, the node determination is
overspecified. This results in an ambiguity in how the nodal positions are defined in the new surface.
One approach that has shown promise for evolving a surface mesh while conserving volume is the method
developed by Jiao.
2, 3
Jiao employs a singular value decomposition (SVD) to solve a least square problem and
then applies an eigenvalue/eigenvector analysis at each node to resolve its normal motion, which generates
the surface geometry, and its tangential motion, which maintains mesh quality.
The first step is to propagate each evolving face in its normal direction. Since the face velocity is given,
a simple, first-order Euler scheme is chosen to integrate along the face normal. This provides the offset
distance for each face.
The second step is to reconstruct the vertices. After computing the new face positions, a new position
for each node on the surface must be determined. For simplicity, assume that the node under consideration
2 of 20
American Institute of Aeronautics and Astronautics
is located at the origin. Each plane passing through a point p with unit normal n can be expressed by a
linear equation n
T
x = δ, where δ = n
T
p. If there are m faces passing through a node, an m × 3 linear
system will be formed
Nx = a. (1)
Here, each row of the system corresponds to one of the m faces that are incident on the node and elements of
a are the offset distances for each incident face. The linear system given by Eq. 1 can be under-determined or
over-determined depending on the value of m. To address this difficulty, a least square solution is computed.
A point is chosen that minimizes the weighted sum of squared distances to the face planes, which is a solution
of the following 3 × 3 linear system:
Ax = b (2)
where A = N
T
WN, b = N
T
Wa, and W is an m × m diagonal matrix with W
ii
equal to the weight
associated with the i
th
face, which is based on the area of the face incident on node p.
Since the matrix A in Eq. 2 is symmetric and positive semi-definite, it has an eigenvalue decomposition
A = VΛV
T
, where Λ is the diagonal matrix consisting of the eigenvalues of A, which are real and non-
negative, and the corresponding eigenvectors are the columns of V. Since A = N
T
WN, the following
singular value decomposition can be derived
√
WN = U
√
ΛV
T
, (3)
where U is a m × 3 matrix.
The vector space spanned by the eigenvectors corresponding to the larger eigenvalues of A is called
the primary space and the complementary space is the null space. An eigenvalue analysis is performed to
identify the primary space. All of the eigenvectors corresponding to eigenvalues smaller than a threshold will
be filtered to avoid instability due to division by a very small number. The nodal displacement is restricted
to the primary space.
The solution to Equation 2 represents an advective motion in which the resulting surface is the intersection
of the propagated face planes. For wavefront motion, such as that produced by burning, erosion, and
deposition, the displacement in the primary space satisfies an entropy condition. Unfortunately, the exact
displacement for wavefront motion can be difficult to compute. A simple solution is to assume the direction
of the displacement will not change and adjust the displacement to satisfy the required offset.
After the displacement is computed, Jiao improves mesh quality by performing a null space smoothing by
computing a tangential motion t at each vertex v. t is a weighted average of the neighborhood of v projected
onto the null space. This smoothing scheme has been shown to preserve sharp features and to introduce
only very small volume errors. Jiao suggests repeating this step without displacement in the primary space
to incorporate a global smoothing into the algorithm that preserves the accreted volume.
III. Approach
Generating an ice shape for a specified accretion time t
ice
is accomplished by performing a series of
quasi-static, loosely-coupled, ice accretion/flow simulation steps. This approach is necessary because, as the
ice shape evolves, it changes the flow field, which, in turn impacts the local ice accretion rate. We term
each of these quasi-static steps an “ice accretion step” with an associated time interval ∆t. As noted below,
each time interval ∆t may be further subdivided into subintervals ∆t
s
. For each ice accretion step: (1)
a CFD simulation is performed to compute the flow field about the current ice shape, (2) a LEWICE3D
computation is performed to determine the new ice shape, or alternatively, the accretion rate map, and (3)
the surface and volume meshes are evolved based on the ice accretion rate using iceSurf and gridMover,
respectively. Here, we focus on the the process employed to compute the deformation of the computational
mesh in response to the ice accretion, which can be divided into three distinct phases:
• Generate a volume accretion rate map on the wetted surface
• Evolve the surface mesh based on the accretion rate map using iceSurf
• Deform the volume mesh by projecting the surface deformations into the volume mesh using gridMover
3 of 20
American Institute of Aeronautics and Astronautics
Each of these processes is described in the sections that follow.
III.A. Generate Accretion Rate Map using Lofting
Currently, LEWICE3D does not provide an accretion rate map. LEWICE3D generates ice shapes using the
strip-based strategy employed in LEWICE2D,
4
which is based on the Messinger icing model.
5
In future
generations of LEWICE3D, a fully three-dimensional approach will be employed to generate an ice accretion
rate map that will be used by the surface mesh evolution algorithm. Two different strategies were employed
to circumvent this shortcoming.
The first approach, which is applicable only for cold, rime icing accretions, uses the collection efficiency
to estimate a surface icing rate. The underlying assumption employed in this approach is that the droplets
freeze on impact producing a pure deposition problem. The icing rate for each surface element is calculated
assuming that no evaporation or runback occurred and is given by
dv
ice
dt
=
V
∞
× β × LW C × A
ρ
ice
(4)
where V
∞
is the freestream airspeed, β is the local collection efficiency, LW C is the free stream liquid water
context, A is the area of the surface element under consideration, and ρ
ice
is the ice density.
The second approach, which is applicable for warmer, glaze icing conditions, uses a lofting method to
generate the icing rate map. The lofting algorithm assumes that the ice thickness varies linearly along
spanwise lines for wings and circumferentially for bodies of revolution (e.g. inlets, spinners and radomes).
Lofting information is used to generate transformations that facilitate interpolation of ice thickness from
the strip-based ice accretions to the surface. This method is flexible and robust and allows interpolation on
wings with taper, twist, and leading edge curvature and bodies-of-revolution.
The ice patch lofting scheme uses the ice thickness values for the surface nodes, which are interpolated
from the LEWICE3D ice shape values, along with the surface normal at the nodes to generate the new iced
surface. Volume elements are formed from the original surface element and the new displaced iced surface
element. The icing rate is then determined by calculating the volume of these iced volume elements and then
dividing this volume by the icing time. The use of ice thickness interpolated from the LEWICE3D ice shapes
allows a convenient, robust method for generating three-dimensional iced surfaces which have run-back and
evaporation effects.
Two types of lofts are available to describe various surfaces of interest. The first lofting type is a wing-
type lofting. This lofting requires the input of the leading edge and trailing edge of the wing. The second
lofting type is a body of revolution-type lofting. The body of revolution lofting requires the line of rotation
and the leading edge center of rotation. A two-dimensional coordinate system (S,T ) is employed for the
loftings for which S is the axial coordinate and T is either the spanwise coordinate for wing-type lofts or the
circumferential angle for body of revolution-type lofts (figure 2).
Figure 2. LEWICE3D planform types and their associated local coordinate systems – wing(left) and body of revolution
(right).
The position (s,t), and the local coordinate system (u,v,w) at a position (x,y,z) on the lofting are
generated using an iterative process. The up-vector w is the planform normal vector at (S,T ). For wing-
type lofts, the planform normal is generated by taking the cross product of the leading edge direction vector
at T and a vector formed from a line at T between the leading edge and trailing edge lines. For the body of
4 of 20
American Institute of Aeronautics and Astronautics
rotation type geometries the planform normal is the radial vector at T . The spanwise vector, v is generated
by interpolation from the leading edge and trailing edge lines at (S,T ) for wing-type planforms. For body of
rotation type planforms the spanwise vector is the tangential vector at T . The axial vector, u is generated by
taking the cross product of the up-vector and the spanwise vector (figure 2). The local coordinates systems
generated for the surface point and icing cut points are used to transform the local surface points into the
local ice cut point coordinate system for the interpolation of the ice thickness (figure 3).
Figure 3. Local coordinate systems for surface point and ice shape.
Two lines are formed in the local ice cut point system (figure 4). One line segment connects neighboring
ice cut points (l1). The other line is formed using the local ice cut spanwise vector as the slope and the
transformed surface point as the intercept (l2). A set of tests is performed to determine if the minimum
distance between the l1 and l2 occurs within the endpoints of l1 and if this minimum distance is reasonably
small value. If both tests are positive the ice thickness is interpolated linearly from the two ice cut thickness
values at the ice cut endpoints. The local coordinates systems generated for the surface point and icing cut
points are used to transform the local surface points into the local ice cut point coordinate system for the
interpolation of ice thickness (figure 3).
Figure 4. Ice thickness interpolation scheme for surface point.
This procedure is repeated for all of the ice cuts associated with the local surface point. If more than
one intercept is found for the local surface point then a linear interpolation on T is performed from the two
surrounding intercepts (larger T and smaller T than the surface point T ). If the value of T of the surface
point is either greater than the T ’s of all of the ice cut intercepts or is less than the T ’s of all of the ice cut
intercepts the value of ice thickness is set to zero. If an intercept has been found and only one cut has been
associated with the local surface point then an extrapolation is assumed and the ice thickness at the local
5 of 20
American Institute of Aeronautics and Astronautics
