Visual Simulation of Glazed Frost Considering the Supercooled StateTomokazu Ishikawa* Yonghao Yue† Taichi Watanabe* Kei Iwasaki‡
Yoshinori Dobashi§ Masanori Kakimoto* Kunio Kondo* Tomoyuki Nishita¶Tokyo University of Technology †Columbia University (JSPS) ‡Wakayama University/UEI Research§
Hokkaido University/UEI Research ¶UEI Research/Hiroshima Shudo University*ABSTRACTWe propose a method for simulating glazed frost bycomputing heat transfers between water droplets and thesurrounding air. Motions of air and water droplets arecomputed based on a FLIP (Fluid-Implicit-Particle)based fluid solver. Glazed frost is a crystal clear ice andformed from supercooled raindrops that freeze when
they hit object surfaces. We propose a method to createan animation of glazed frost formation by taking into
account the heat transfer between particles and the outside grids. The time to freeze is calculated by considering the heat flux on the surface of the raindrop until
freezing from the adhesion on the obstacles. We propose
a technique which is highly compatible for a FLIP
method to solve the heat conduction equation and wereproduce the formation of glazed frost.
1. INTRODUCTIONGlazed frost is a cluster of clear ice crystals and formed
from supercooled raindrops that freeze when they hit
object surfaces such as the ground and branches of trees(see Fig. 1). Simulation methods for formation of ice
crystal, such as frost, on the surface of objects have been
proposed by Kim et al. [1]. However, a supercooled
state has to be considered for simulating freezing rain,
and fluid simulation is required for reproducing the effect of raindrops running down on the ice surfaces. To
our best knowledge, there has been no research presenting glazed frost by using a fluid simulation. We use thefluid simulation based on a FLIP method [2]. Hence,
raindrops and obstacles are represented by particles
which are used to calculate the advection term, and theupdate of velocity field is calculated by using grids except for advection term. We propose a method to createan animation of glazed frost formation by taking into
account the heat transfer between particles and the grid.
2. RELATED WORKIn the computer graphics (CG) field, research to represent ice and snow began long ago. Kim et al. proposed amethod for generating the crystal geometry of ice by(a) Branches attached by glazedfrost.

(b) Glazed frost attached on one
side of branches.

Fig. 1: Glazed frost (photograph).

solving a heat conduction equation in a two-dimensional
plane [1]. The result of Kim et al. makes it clear thatheat conduction plays an important part in creating ice
crystals.
Takahashi and Fujishiro proposed a method for expressing a state of compacted snow [3]. In this method,they calculated an expression behavior of snow attaching to surface of object by introducing an effect of sintering into SPH (Smoothed Particle Method). In thispaper, we propose to reproduce an icing phenomenon by
calculating the ice buildup as an aggregate of particles
and introducing a sticking power term into the equation.
Iwasaki et al. proposed a method to calculate ice-toice interaction at high speeds using SPH [4]. They
showed that ice can be simulated as a rigid body at highspeed using a GPU. However, although this method canexpress a process in which an icicle forms or water
freezes in a blink, the supercooling phenomenon is nottaken into consideration.
In our research, we adopt a FLIP method for representing a fluid phenomenon. The effectiveness and feasibility of the FLIP method in the field of CG has beensuggested by Zhu et al. [2]. Recently, it has been actively improved for use in animation, and it is broadly utilized in various industries for reproducing a fluid statesuch as water or flame.
Numerical calculation of supercooled water has also
been executed in the field of computational physics. Inaddition, a report has verified the thermal ambience of
air against the actually observed glazed frost phenomenon [6]. In this paper, we confirm that temperature, humidity, wind velocity, and precipitation amount create
the conditions for developing glazed frost. InFig. 2: Process flow in each timestep.

accordance with these facts, we propose a feasible model to realize a visual simulation of glazed frost.
3. PROPOSED SIMULATION METHOD
In this section, we explain the modeling of glazed frost
and calculation of heat transfer using a FLIP method.
3.1. Fluid Calculation of Rainwater
The fluid behavior of rainwater is calculated by the Navier–Stokes equations. We solve the equations using theFLIP method, a hybrid method to solve the following
Navier–Stokes equations using a grid and particles:(1)∇ ⋅ u = 0,∂u
1(2)
+ ( u ⋅∇ ) u = − ∇p + µ∇ 2u + F,∂tρ
where u is velocity, t is time, ρ is density, p is pressure,µ is dynamic coefficient of viscosity, and F is external
force such as wind and adhesion by surface tension (details are described in Sec 3.3). In the FLIP method, updating the velocity field by the external force term andpressure term is calculated using a grid method, and theadvection term is calculated using a particle method.
In our research, rainwater is represented as particles.
Velocity vi, temperature Tp,i, and state Si are stored in thecenter of the i-th particle with a radius rd. Particles entering the simulation space from the top represent
raindrops. The size of the raindrop is determined by thefollowing Marshall–Palmer distribution [5].
(3)=N ( D ) N 0 exp ( −λ D ) ,where D is the size of the raindrops, and N0 and λ are
the initial and gradient parameters, respectively. A
raindrop with particle size D is approximated by an appropriate number of particles that can be equal to thevolume of a raindrop. That is, we determine a number ofparticles for raindrops of particle size D to a maximumof natural number n that satisfies the following nonstrict
inequation:
D3(4)
n≤ 3.
8rd
The raindrops of small size D are not considered when
the natural number n satisfying the above nonstrict inequation does not exist.
Water particles in a supercooled state are determined to be solidified or not by considering the amountof heat. Some water particles get solidified at the time ofimpact with an obstacle. Supercooled water that is not
solidified at the time of impact can flow down. We also
process both changing of ice to water and water to ice
by ambient air temperature. Sublimation (changing fromwater vapor to ice) can result in ice crystal formation,
but we do not consider it since the amount of ice crystal
formation by sublimation is quite small. We represent
objects covered by glazed frost (e.g., plants and branches) with particles that are referred to as object particles.
The only particles updated by the FLIP method are thewater particles; object particles and accreted ice particles are set as fixed boundaries in fluid calculation.
Temperature should be maintained even for object particles, and then we calculate a heat transfer with the ambient air. Furthermore, we do not include any computation
regarding a process for creating supercooled water, but
we start a simulation in which all rain falling in the simulation space is considered to be in a supercooled state.
3.2. Heat Transfer in a FLIP Method
We propose a method to calculate heat conduction thatis suitable for the FLIP method. The process flow during
a single timestep is shown in Fig. 2.
We distribute velocity from particle to lattice in thefirst step of FLIP method. Heat is distributed from particle to grid as the velocity field is distributed. Heat transfer is not conducted between particles, but instead iscalculated between lattices for computational efficiency.
Heat is distributed from grid to particles after the heattransfer calculation. In addition, adhesion is calculated
after supercooled water collides with object particles.
The time from adhesion to a collided object’s surface tobeing frozen is calculated with consideration of heat flux
for glazed frost surface.
The heat distribution between particles and grid iscalculated according to Fourier’s law:
∆Qi =− ki A

Tair − Tp.i
dTi,=− ki Adxrd(5)where ∆Qi is heat quantity transferring between the i-thparticle and grid including the particle, ki is thermalconductivity, A is surface area of the particle, dx is theinfinitesimal distance, dTi is the temperature difference
between the particle and ambient air, and Tair is temperature at the grid point near the particle. The temperatureof each particle is updated by the heat quantity,∆Q(6)∆Ti = i ,mi ci
where mi is mass of particle i and ci is specific heat ofparticle. Mass and specific heat of particle are set todifferent values depending on whether the particle represents a water droplet or an obstacle.
Heat transfer is calculated on the grid by solving
heat transfer equations.
∂T∂ 2T(7)=a 2 ,∂t∂x
where a is thermal diffusivity coefficient and x is a position vector. After calculating heat transfer, the heat is
distributed to the particles. At update of the velocity
field, a timestep is set according to the Courant–Friedrichs–Lewy conditions. To use this timestep in thenumerical solution of heat transfer equations, we employ
the Crank–Nicolson method, an implicit scheme that canstably solve at arbitrary time intervals, to handle the heattransfer equation.

3.3. Calculation of adhesion force
For describing an ice accretion when rain in the supercooled state collides with an object such as plant or withalready-formed glazed frost, we calculate a surface tension and interfacial tension for object surface. We take
into account the surface tension and interfacial tension
as external forces in Equation (2). These forces are calculated by the following equation.
p j − pi
(8),
Fisticking = ∑ κ j2j∈Nip j − pi
where pi is position of particle i and κi is coefficient ofadhesion force. Ni is the set of particles which exists in alattice with particle i and 6 lattices adjacent to it, and theadhesion should be considered with an influence fromparticle j with the accumulation. We improve the calculation efficiency due to collision detection limited in Ni.
3.4. Freezing Condition
The freezing time is calculated by considering the heat
flux on the surface of the raindrop until it freezes fromadhesion to the obstacles. The freezing speed is calculated by the following sum of heat flux:(9)
Q = Qs + Ql + Q f ,where Qs is a sensible heat flux, Ql is a latent heat flux,
and Qf is the heat flux necessary to freeze a water film.
These quantities are calculated using the following equations:(10)
Qs =−π ha ∆T ,Ql =−π Le hv ∆ρv ,(11)Q f = L f w,

(12)

where ha is a heat exchange coefficient of air in Equation (10). In Equation (11), hv is a water vapor exchange
coefficient, Le is evaporative latent heat, and ∆ρv is difference of water vapor concentration. In Equation (12),Fig. 3: Simulation space for verification experimentFig. 4: Evaluation by RMSE between exact solution and calculation result.
Table 1: Spatial discretization and computation time.

Lf is a latent heat for freezing and precipitation flux, andw is calculated by the following equation using the relationship of rain rate R and wind speed V.

{w=( R / 3.6 ) + ( 0.067 R 0.846V )2} ×10 .

12 2−3(13)It is known that the smaller the sum of flux Q in
negative value, the faster the freezing speed [6]. So weuse the time integrals Qs and Ql and the required amountof heat Hf as the condition to determine whether the water freezes.
(14)
H f ≤ Qs + Ql dt.

∫Integration calculation is limited to the time the water isattached to an obstacle or other glazed frost. The water
does not freeze if it drifts by other raindrops and is affected by gravity or is not attached to the surface ofglaze frost.
4. RESULTS AND CONCLUSION
For the simulation, we used C++ and the nVIDIA
CUDA on a standard PC (CPU: Intel Core i7 3.2 GHz,
RAM: 16.0 GB, GPU: nVIDIA GeForce GTX 580). For
rendering the results, we used Autodesk 3ds Max.
First, we conducted experiments to verify the calculationof heat transfer by the proposed method. We calculated
the error of numerical calculation in cases that we canexactly solve a heat conduction in two dimensions, suchas the following: A circular object is located at thecenter of the simulation space, the transferred heat
between the outside air and the object is calculated (seeFig. 3), and then the temporal change of the temperature
is calculated. The exact solution at time t in polar
coordinates (r, θ) on the circle in this case is calculatedwith the following formula:T ( r ,θ , t )∞∞∑∑ en 0=m 1=− µ m ,n t

J n ( µm ,n r ) ( Amn cos nθ + Bmn sin nθ ),(15)where Jn is n-th Bessel function and μm,n is m-th value
from the smallest positive zero point of Jn. The proposed
method discretizes the simulation space with particles
and a grid. Table 1 compares space discretization andcalculation time for both methods in the case of thetemperature value with the exact solutionIn addition, thecomparison of the RMSE (root-mean-square error)
between the exact solution and our calculation result isshown in Fig. 4. The results confirm that the proposed
method is comparably as stable as the numericcalculation with a grid method, and it can compute witha small error. Since the error in the exact solution isespecially significant near the boundary, the boundary
condition of the proposed method needs moreconsideration.
In Fig. 5, the brown particles represent tree
branches and the deep blue particles represent ice(glaze). The number of particles used for the simulationis almost 20,000, including the objects such as fluid andbranches. The grid resolution of the simulation space is128 × 128 × 128, and the calculation time for 1 timestep
is 20 msec. We show parameters used in our simulationin Table 2.
Fig. 6 shows rendering results produced by themental ray rendering engine used with 3ds Max. The
surfaces of the glazed frost are calculated by assigning
particles representing glazes to anisotropic kernel
density distributions proposed by Yu et al. [7] andextracting isosurfaces. The surface extraction took 30
seconds on average, and the rendering took almost 3
minutes for 1 frame.
Supercooled water in some cases freezes instantly
after collision, and in some cases gradually condenses in
the process of running down objects. This effect can be
expressed by calculating the thermal environment,sensible heat, and latent heat. Our method can reproduce
the ice covering a branch.
By changing the parameter of the thermalenvironment, the situation was simulated that when theenvironmental temperature is higher than 0°C, water did
not adhere at the top part of a branch and instead froze
while it was flowing down (Fig. 7). Because thesupercooled state is considered, the icing when colliding
with the object surface and the fluid phenomenon of thechange from ice to water can be calculated to the formation of icicles. Moreover, external forces such as

wind can be considered using the FLIP method. The
velocity field of the simulation space was defined where
the wind blew from the right side of the image. The statein which more glaze is formed on the windward side of
the tree branch can also be reproduced (Fig. 8).
We have proposed a technique that is highly compatible with a FLIP method to solve the heat conduction
equation, and we used it to reproduce formation ofglazed frost. By considering the supercooled conditions,
our method can simulate formation of icicles by calculating ice buildup after collision with obstacles and thefluid phenomenon of changing of ice to water.
Table 2: Parameters in our simulation.
parameter

meaning

value

unit

dt

timestep

1.0sρ

density

1.0×103kg/m3µ

dynamic coefficient ofviscosity
initial value of Marshall
Palmer distribution
gradient of Marshall
Palmer distribution

1.792×10－6

m2/s

8.0×103mm－1m－3

0.1192×103

m－1rd

radius of particle

1.0×10－3

mk

thermal conductivity

0.569

W/(m・K)N0
λc

specific heat

4217

J/(kg・K)

thermal diffusivity coefficient
coefficient of adhesion
force (liquid)
coefficient of adhesion
force (solid)
heat exchange
coefficient of air
water vapor exchange
coefficient

0.1349×10－6

m2/s

0.5

N

2.0N

5.69

W/(m2・K)

0.05m/sLeevaporative latent heat

2.498

J/kg

∆ρv

difference of water vapor
concentration

1.61×10－3

kg/m3

Lf

latent heat for freezing

333.4×103

J/kg

R

precipitation intensity

25×10－3m/h

Hfrequired heat quantity
for freezing

333.4

J/m3aκκhahv(a) 10 timesteps(b) 7,000 timesteps

Fig. 5: Simulation results of generating glazed frost.

(a) 1,000 timesteps(b) 4,000 timesteps

Fig. 7: Simulation result of formation of icicle shapes.

Fig. 8: Simulation result of glazed frost adhering on the windward sideby the wind from right to left.

REFERENCES[1][2][3][4](c) 7,000 timesteps
[5]
[6]

[7]

T. Kim, M. Henson and M. C. Lin, “A Hybrid Algorithm for
Modeling Ice Formation”, Proceedings of the 2004 ACM
SIGGRAPH/Eurographics Symposium on Computer Animation,
pp.305–314, 2004.
Y. Zhu and R. Bridson, "Animating Sand as a Fluid", Proceedings of the ACM SIGGRAPH 2005, pp.965-072, 2005.
T. Takahashi and I. Fujishiro, “Particle-based simulation ofsnow trampling taking sintering effect into account”, ACM
SIGGRAPH 2012 Posters, No. 3, 2012.
K. Iwasaki, H. Uchida, Y. Dobashi, T. Nishita, "Fast Particlebased Visual Simulation of Ice Melting", Computer Graphics
Forum, Vol. 29, pp.2215-2223, 2010.
J.S. Marshall and W.M. Palmer, "The distribution of raindrops
with size," Journal of Meteorology, Vol. 5, pp. 165-166, 1948.
K. F. Jones, "Ice accretion in freezing rain", Tech. rep., DTIC
Document, 1996.
J. Yu, G. Turk, "Reconstructing Surfaces of Particle-BasedFluids Using Anisotropic Kernels", Proceedings of the 2010
ACM SIGGRAPH/Eurographics Symposium on Computer Animation, pp.217-225, 2010.

APPENDIX(d) 10,000 timesteps
Fig. 6: Simulation and rendering results of generation of glazed frost.

A. COEFFICIENTS OF THE EXACT SOLUTION
IN THE THERMAL DIFFUSION EQUATION OF
THE DISK AREAIn Equation (15), Jn is n-th Bessel function and definedby the following equation,

( −1)  z n+ 2 k .
 k = 0 k !( n + k ) !  2 k∞Jn ( z ) = ∑(A1)

Parameters in Equation (15) are the constants determined by the following equations.
Amn =2π J n +1 ( µm ,n )2∫ ( ∫ π T ( r ,θ , 0 ) J ( µ r ) cos nθ dθ ) rdr ,1π0−nm,n(A2)

Bmn =

2

π J n +1 ( µm ,n )2

∫ ( ∫ π T ( r ,θ , 0 ) J ( µ r ) sin nθ dθ ) rdr.
1π0−nm,n(A3)