以粒子為基礎的物體融化模擬與流動控制

國

立

交

通

大

學

多媒體工程研究所

碩

士

論

文

以粒子為基礎的物體融化模擬與流動控制

A Particle-based Simulation for

Object Melting and Flowing Control

研 究 生：羅建凱

指導教授：莊榮宏 教授

指導教授：

林文杰 教授

以粒子為基礎的物體融化模擬與流動控制

A Particle-based Simulation for Object Melting and Flowing Control

研 究 生：羅建凱 Student：Chien-Kai Lo

指導教授：莊榮宏 Advisor：Jung-Hong Chuang

指導教授：

林文杰

Advisor：

Wen-Chieh Lin

國 立 交 通 大 學

多 媒 體 工 程 研 究 所

碩 士 論 文

A Thesis

Submitted to Institute of MultimediaEngineering College of Computer Science

National Chiao Tung University in partial Fulfillment of the Requirements

for the Degree of Master

in

Computer Science

September 2010

Hsinchu, Taiwan, Republic of China

i

以粒子為基礎的物體融化模擬與流動控制

研究生 : 羅建凱 指導教授 : 莊榮宏 博士

研究生 : 羅建凱 指導教授 :

林文杰 博士

國立交通大學

資訊學院多媒體工程研究所

摘 要

A Particle-based Simulation for

Object Melting and Flowing Control

Student: Chien-Kai Lo

Advisor: Dr. Jung-Hong Chuang

Dr. Wen-Chieh Lin

Institute of Multimedia Engineering

College of Computer Science

National Chiao Tung University

ABSTRACT

We purpose an SPH method for object melting simulation and a feedback control algorithm to manipulate the melting flow. By adding an elastic stress term to the Navier-Stokes equation, we can simulate the melting behavior which has both fluid and solid features. We develop a thermals transfer system to simulate the heat flow between particles, and based on the temper-ature of particles to determine the viscosity and elastic stress coefficient automatically. Also, we introduce a feedback control algorithm which can be used to design the flowing behavior intuitively. Users can draw flowing path directly on the surface to design the shape of melting fluid. Our experiments shows that the proposed approach can preserve more details of melting than previous methods, and be useful for authority a desired a melting scene.

Acknowledgments

I would like to thank my advisors, Professor Jung-Hong Chuang and Wen-Chieh Lin, for their guidance, inspirations, and encouragement. I am also grateful to my senior colleague, Tan-Chi Ho, for his advices and suggestions on my thesis research. And then I want to express my gratitude to all my colleagues in CGGM lab: Jau-An Yang and Ta-En Chen discusses with me on research issues and helps me figure out the questions, Ying-Yi Chou and Chien Yun give me several thesis advices and some necessary information, Chia-Ming Liu and Yu-Chen Wu often share some useful programming skills and knowledge with me, Tsung-Shian Huang help me understand much rendering knowledge, and all my junior colleagues’ kind assistances. Also, I want to thank my dear girlfriend Yi-Yu Lin for her support through all these years. Lastly, I would like to thank my parents and my sister for their love, and support. Without their love, I could not pass all the frustrations and pain in my life.

Contents

1 Introduction 1

1.1 Contribution . . . 2

1.2 Organization of the thesis . . . 3

2 Related Work 4 2.1 SPH-based Fluids . . . 4

2.2 Melting and Flowing . . . 7

2.3 Fluid Control System . . . 8

3 Simulation and Control Framework 11 3.1 Overview . . . 11

3.2 Fluid Simulation Framework . . . 13

3.2.1 SPH formulations . . . 13

3.2.2 Viscoelastic Fluid Model . . . 16

3.3 Heat Transfer System . . . 18

3.3.1 Previous Thermal Diffusion Method . . . 23

3.4 Adaptive Particle Sampling . . . 24

3.5 Flowing Control System . . . 29

3.5.1 Path Control . . . 31

3.5.2 Shape Control . . . 34

3.6 Implementation Details . . . 37

4 Results 40 4.1 Simulation Statistics . . . 40

4.2 Simulation Results . . . 41

4.2.1 Object Melting Scenes . . . 42

4.2.2 Thermal Transfer Simulation Results . . . 43

4.2.3 Fluid Control Results . . . 45

4.2.4 Comparison with Previous Methods . . . 48

5 Conclusions 50 5.1 Summary . . . 50

5.2 Limitations and Future Work . . . 51

Bibliography 52

List of Figures

1.1 Real photo of candle melting vs. Our candle melting result . . . 2

2.1 Three smoothing kernels Wpoly6, Wspikyand Wviscosity(from left to right). The thick lines are the kernels, the thin lines are the gradients in the direction to-wards the center, and the dash lines are the Laplacian [MCG03] . . . 5

2.2 Extended local feature size [APKG07]. . . 6

2.3 Thermal transfer in [MGG+10] . . . 8

2.4 Particle-based control system [TKPR06] . . . 10

3.1 System flow of our algorithm . . . 12

3.2 Fluid solver flowchart . . . 14

3.3 Thermal transfer system flowchart. . . 18

3.4 Thermal simulation in our approach. . . 19

3.5 Surface particle and its normal. . . 20

3.6 Thermal radiation diagram . . . 21

3.7 Candle melting without thermal radiation and dissipation . . . 24

3.8 Adaptive particle sampling (on the right) compared to uniform sampling (on the left). . . 27

3.9 Particle resampling scheme in [APKG07]. . . 28

3.10 Flowing control flowchart. . . 30

3.11 Two path control forces. . . 32

3.12 Modify the adoption coefficient around the path starting point . . . 34

3.13 Shape control force for each node. . . 35

3.14 Shape node two level influence range . . . 36

3.15 Signed distance value in four different object . . . 38

4.1 Ice cream melting . . . 42

4.2 Chocolate bunny melting . . . 43

4.3 Thermal simulation in butter scene . . . 44

4.4 Final image in butter scene . . . 45

4.5 Candle melting in path control . . . 46

4.6 The shape control line drawn by user . . . 47

4.7 Butter melting in shape control . . . 47

4.8 Wax melting in [CMVHT02] . . . 48

4.9 Candle melting in [PPLT09] . . . 49

4.10 Candle melting in our approach . . . 49

List of Tables

4.1 Simulation statistics of different scenes . . . 41 4.2 The percent of time consuming in different stages . . . 41

C H A P T E R

1

Introduction

Physics-based simulation has been developing in Computer Graphics field for many years. In particular, fluid simulation has been playing an important role in those works. With the help of these photorealistic results, more special effects can be designed by artists, even it is not exist in the reality. Although many stunning results had been presented in recent years, there are still many fields need to be improved. Such as the fluid detail preserving, intuitive fluid control method, and real-time rendering applications. These topics still remain a great challenge for researchers. This thesis focus on the object melting behavior and flowing control. The melting process is a complex fluid behavior in real world, this natural phenomenon including the flowing of the melting fluid, phase change between solid and fluid, and the heat transfer in the simulation elements.

In the past few years, many different methods had been purposed which may deal with the melting manner. But there are several problems remaining in previous works. First, the previous simulation didn’t capture the detail of fluid, such as the droplet flow along the candle surface as in Figure 1.1. Second, rare of them considered the thermal radiation effect, which is a dominant part in thermology. And third, those methods were lack of control mechanics, thus it was hard

1.1 Contribution 2

for animators to design the way for object melting and the shape they want to achieve.

Figure 1.1: Real photo of candle melting vs. Our candle melting result

In our method, We develop a particle-based fluid solver and a physics-inspired heat transfer system to solve the object melting problems as mention above. Also, we provide a straight-forward fluid control system, which allows animators to design the flowing pattern and the melting paths. For the purpose of preserving more melting detail, we hire an adaptive particle sampling scheme, enabling us to allocate the available computational resources on the regions of interest dynamically.

1.1

Contribution

The contributions of our melting simulation can be summarized as follows:

• Use an adaptive sampling scheme to enrich the melting detail without using too much particle numbers.

1.2 Organization of the thesis 3

• Provide an intuitive controlling framework, which allows user to design the melting flow and the final shape easily.

• May perform a wide range of material types from high viscoelastic fluid to normal New-tonian fluid.

1.2

Organization of the thesis

The rest of the chapters are organized as follows. Chapter 2 gives a literature review for SPH fluids, melting and flowing ,and related fluid control approaches. Chapter 3 illustrates the main algorithm in our method, includes the particle-based simulation framework, thermal transfer simulation in fluid elements, adaptively particle sampling scheme in object melting, and a feed-back control of fluid behavior. Chapter 4 shows the experimental results of our approach. Lastly, we summarize a conclusion and discuss the future work in chapter 5.

C H A P T E R

2

Related Work

In this chapter, we briefly describe the related works by several parts: SPH-based fluids, melting and flowing, fluid control system.

2.1

SPH-based Fluids

The Smoothed Particle Hydrodynamics(SPH) method, originally devised for the simulation of stars [Mon05], was first used by Desbrun et al. [DC96] in computer graphics field to solve the interaction forces of large deformation. This method can handle complex large topological deformation and satisfies the mass conservation equation easily.

Later, M¨uller et al. [MCG03] introduced a particle-based fluid simulation scheme for inter-active applications. By deriving the equation of pressure, viscosity and surface tension term, and using three different smoothing kernels (Figure 2.1), they show that SPH can simulate the fluid behavior appropriately. Adams et al. [APKG07] provides an adaptively sampled particle scheme, that performs the particle merge and split operation based on the extended local feature size(ELFS). The classical local feature size of a point x,LFS(x), on a (closed) manifold M is

2.1 SPH-based Fluids 5

Figure 2.1: Three smoothing kernels Wpoly6, Wspikyand Wviscosity(from left to right). The thick

lines are the kernels, the thin lines are the gradients in the direction towards the center, and the dash lines are the Laplacian [MCG03]

.

defined as the distance of x to the closest point on the manifold’s medial axis. As show in Figure 2.2, they use particle approximation of the medial axis to compute the ELFS value. Red regions in Figure 2.2 (c) have a lower feature size and should be sampled densely while blue regions have a larger feature size and can be sampled coarsely. The extended local feature size of x, ELFS(x), is defined as the minimal sum of the distance from x to the point y on the manifold and the classical local feature size at point y (Equation 2.1). By utilizing the ELFS value, the region with different local feature level can be clearly classified. They used larger particles in the region has larger ELFS value, which means the region is inside the volume or near a thick flat surface, to reduce the number of particles. Thus, the simulation speed may get a significant improvement.

2.1 SPH-based Fluids 6

Figure 2.2: Extended local feature size [APKG07].

ELF S(x) = min

y∈M{kx − yk + LF S(y)}, x ∈ ν. (2.1)

Figure 2.2 (a) shows the fluid volume with the medial axis. The color shows the distance to the fluid surface at any given position. Figure 2.2 (b) approximates the medial axis with particles. And Figure 2.2 (c) shows the color-coded ELFS value.

Moreover, there are several studies on more complicated fluid phenomena. Lenaerts et al. [LAD08] introduced an SPH algorithm to simulate the flow through the porous deformable material and showed some interesting demos, such as the simulation of sponge-like elastic bodies and water-absorbing sticky cloth. Lenaerts et al. [LD09] also simulated the interactions between fluids and granular materials by purposing a unified SPH framework. They treat the granular material as a continuous volume sampled by particles, which can be animated like SPH fluids. Solenthaler et al. [SP09] enforced incompressibility by purposing a prediction-correction scheme to solve particle pressure. They avoided the computational cost of solving a pressure Poisson equation, and achieved an order of magnitude speed-up compared to the commonly used weakly compressible SPH.

SPH has been more and more widely used in fluid simulation. By representing the fluid with a set of particles, we may access individual material properties directly and move the particles according to the Navier-Stoke equation easily. A detailed and extensive introduction to SPH can be found in [Mon05].

2.2 Melting and Flowing 7

2.2

Melting and Flowing

Object melting and flowing has many interesting effects, and it is still remaining a challeng-ing field for researchers. Carlson et al. [CMVHT02] modified the Marker-and-Cell (MAC) algorithm to compute the fluid variable and arbitrarily high viscosity. They also introduced a heat diffusion equation to handle the change in temperature, and varies the material’s viscosity according to the temperature. Keiser et al. [KAG+05] proposed a unified model to solve the de-formation of solids and fluids, by combining the solid mechanics with Navier-Stokes equation in a Lagrangian form. As shown in Equation 2.2, they applied the same concept as heat transfer equation in [CMVHT02] but in a SPH manner.

DTi Dt = η X j (Tj − Ti)∇2W (r − rj, h) mj ρj (2.2)

Paiva et al. [PPLT06] [PPLT09] simulated the object melting effect by adding a viscoplas-tic stress tensor term to the traditional Lagrangian fluid formulation. They used the same SPH concept in [MCG03] but with different smoothing kernels. Instead of three different smoothing kernels, one unified smoothing kernekl and a XSPH correction scheme is used to reduce the particles instability due to the arbitrarily small distance between particles. They modeled the stress tensor with the Generalized Newtonian Fluid model, and simplified the viscosity control by introducing a new rheological parameter called jump number. Chang et al. [CBL+09] ani-mated the viscoelastic fluid by adding an elastic stress to the momentum conservation equation (a modified version of Navier-Stokes equation).Their method also showed the ability to simu-late the object melting and flowing behavior by using the concept similar in [CMVHT02] where the viscosity coefficient µ and elastic stress coefficient µsare influenced by the temperature. All

the methods mentioned above solve the heat diffusion equation and varied the viscosity depend-ing on the temperature, but none of them considered the heat radiation effect, which plays an important role in heat transfer system. Wei et al. [WLK03] employed a linear 3D cellular au-tomata approach on the solid volumetric models to animate the melting process. By reducing

2.3 Fluid Control System 8

the physical simulation rule, they may animate the smooth fluid behavior in an interactive rate. Mar´echal et al. [MGG+10] developed a physically based heat transfer simulation system, to synthesize realistic winter sceneries. They defined the environmental model and three major types of thermal transfers in different environmental elements. Figure 2.3 shows the thermal transfer they had taken into account, including heat conduction, heat convection, and heat ra-diation in different materials. Our thermal transfer system is similar to their’s, but without the heat convection term that is less important in normal object melting.

Figure 2.3: Thermal transfer in [MGG+_10]

2.3

Fluid Control System

Simulation of realistic fluid behavior has been always a hot field in computer graphics. Re-cently, more and more researchers have focused on the fluid control mechanism for artistic reasons in artist. Foster et al. [FF01] was the first to apply control scheme on fluids. They produced 3D parametric space curves to locally control the velocity and pressure of the fluid flow. Rasmussen et al. [REN+04] defined four types of control particles according to the liquid variables, including velocity, viscosity, level set, and the velocity divergence particles, where velocity particles controlled the movement of the liquid interface; viscosity particles for viscous melting behavior; level set particles is used to modified the level set representation in liquid

in-2.3 Fluid Control System 9

terface; divergence particles modeled the liquid expansion and contraction effects. Treuille et al. [TMPS03] enforced the smoke to match the given velocity and density keyframes by building an optimization technique to solve the control parameters. McNamara et al. [MTPS04] used the adjoint method to improve the efficiency of solving the nonlinear optimization problem. The key to their approach is to utilize the linear duality for variables substitution. The adjoint method exploited this powerful aspect of duality to accelerate the gradient calculation in opti-mization process. Fattal et al. [FL04] controlled the smoke shape with a sequence of predefined target smoke states. They added two terms to the standard flow equations : driving force term and smoke gathering term. While the driving force compelled the fluid to a particular smoke shape, the gathering term gathered the fluid to maintain the density value. Instead of solving the control forces optimization problem, Hong et al. [HK04] created a geometric potential field to exert the external force, and coerce the fluid to form the target shape. Shi et al. [SY05] controlled the liquids to match with the rapidly changing targets by building two different ex-ternal force fields. The first one is a feedback force field applied to both shape and velocity to compensate the discrepancies with the target field. The second one is the gradient field of an adaptive geometric potential field defined by the target object shape and skeleton informa-tion. Compared to the original geometric potential field in [HK04], their adaptive potential field could handle different complex target shape with appropriate force magnitude. Later, Th¨urey et al. [TKPR06] introduced a control particles system to control the liquid without losing the detail part. As show in Figure 2.4, they first generated control particles set using some prede-fined criteria, then computed the attraction force and velocity force acting on the general fluid elements separately. In the end, they combined the control forces with weighting sum to evalu-ate the final forces for each element. In order to allevievalu-ate the virtual viscosity effect, they also applied high and low pass filter to preserve the detail motion. More recently, Dobashi et al. [DKNY08] developed a cloud formation control system. They provided a geometric potential field similar to the previous method [HK04] to solve the horizontal shape, but used a different feedback control scheme to control the cloud vertical height. Instead of feedback to the force strength, their feedback system controlled a water vapor supplier and a latent heat controller,

2.3 Fluid Control System 10

Figure 2.4: Particle-based control system [TKPR06]

allowing them to deal with the phase change effects and the cloud generation.

Though there are many fluid control methods, none of them had dealt with object melting and flowing problem, which is more complex because of the heat transfer involved in fluid elements and the melting behavior.

C H A P T E R

3

Simulation and Control

Framework

3.1

Overview

The aim of our algorithm is to simulate the object melting naturally, and also provide a control system to design the melting scene. First, the melting behavior should be triggered when the temperature reaches the melting point, and the material could transform from solid into fluid smoothly. Second, the simulation ought to preserve the detail drops flow on the surface. Third, although the simulation follows the physics law, animators still need certain control ability to design the flow as they wish.

Figure 3.1 is the system flow of our algorithm. In the initial stage, user can author a melting scene as a input to our control system. In the simulation part, we first solve the Navier-Stokes equation using SPH method. Also, by adding an elastic force term to the momentum conser-vation equation, we are able to simulate the high viscous fluids and viscoelastic fluids. After solving SPH, we simulate the thermal transfer between particle elements. We have modeled

3.1 Overview 12

3.2 Fluid Simulation Framework 13

three types of thermal transfer methods: thermal radiation, thermal conduction, and thermal dissipation. According to the temperature of the particle, the phase change effect can be trig-gered by changing the viscosity and elastic coefficients. Besides, we employ an adaptive parti-cle sampling scheme to improve the efficiency, allowing us to use more computational resource for the melting region. Based on the output of the simulation part, we feedback to the control mechanism to adjust the force field automatically. Then, we can move the particles by an Euler scheme. In the end, we reconstruct the surface by the marching cube algorithm and render the polygon mesh for output images.

3.2

Fluid Simulation Framework

In this section, we describe the physics framework in the fluid solver. Figure 3.2 shows the system flow of fluid solver. We first solve the Navier-Stokes equation with viscoelastic force in the SPH manner and an adaptive sampling scheme.

3.2.1

SPH formulations

SPH is an interpolation method defined on the particle system. In SPH, fluids are represented by discrete particles and all the fluid properties are carried on those particles. As shown in Equation 3.1, one can evaluate a particular quantity A at any position r according to the weighted sum of contributions from the neighboring particles j.

A(r) =X j Aj mj ρj W (r − rj, h) (3.1)

Where r and rjrepresent a any given position and the neighboring particle j position, mj means

the mass define on the neighboring particle j, ρj is the density define on the particle j, Aj is a

scalar value define on the particle j. The W is the smoothing kernel with smoothing length h. The gradient and the Laplacian of scalar A can be calculated by using the following gradient and the Laplacian of the smoothing kernel W (Equation 3.2, Equation 3.3):

3.2 Fluid Simulation Framework 14

Figure 3.2: Fluid solver flowchart

∇A(r) =X j Aj mj ρj ∇W (r − rj, h) (3.2) ∇2_{A(r) =}X j Aj mj ρj ∇2_{W (r − r} j, h) (3.3)

Before calculating the force fields, we need to have the density and the pressure information for every particle. The density ρ (in Equation 3.4) and the pressure (in Equation 3.5) of the

3.2 Fluid Simulation Framework 15

particle i are calculated as in [MCG03]:

ρi =

X

j

mjW (rij, hi) (3.4)

Pi = k(ρi− ρ0) (3.5)

where rij = ri − rj, k corresponds to the gas constant, ρi is the density of particle i and ρ0 is

the rest density. The pressure of each particle is given according to the ideal gas state equation [DG96].

The fluids in SPH method need to satisfy the Navier-Stokes equation, including the mass conservation equation (Equation 3.6) and the momentum conservation equation. (Equation 3.7). We add an elastic stress tensor force to the traditional momentum conservation equation in order to model the high viscosity and viscoelastic fluids properties.

∂ρ

∂t + ∇ · (ρv) = 0 (3.6)

ρ(Dv

Dt) = −∇p + µ∇

2_{v + f}

surf ace+ f + felastic (3.7)

There are five force field in the right hand side of Equation 3.7. −∇p is the pressure term, µ∇2_{v is the viscosity, f}

surf aceis the surface tension force field, f is often thought as the external

force (ex: gravity force, f = ρg), and felasticis the elastic strain tensor force field which will be

detailed in subsection 3.2.2. Following the adaptively sampled particles scheme in [APKG07], we use a symmetric fluid force to avoid asymmetrical force distribution. We will explain our adaptive sampling system in the next section. The fluid force acting form particle pj on pi in

Equation 3.7 can be obtained as follows:

f_ijpressure = −∇p = −mi ρi mj ρj (Pi+ Pj) 2 (∇W (rij, hi) + ∇W (rij, hj)) (3.8) f_ijviscosity = µ∇2v = −mi ρi mj ρj (vj − vi) 2 (µi∇ 2_{W (r} ij, hi) + µj∇2W (rij, hj)) (3.9)

3.2 Fluid Simulation Framework 16 f_ijsurf ace=    −σ 2( ni knik + nj knjk)∇ 2_{W (r} ij, hi), if knijk ≥ δ 0, if knijk < δ (3.10) where nij = m_ρii mj

ρj∇W (rij, hi), and P is the pressure, rij = ri − rj, v is the velocity

vector, µ is the viscosity coefficient, σ is the surface tension coefficient, and ni is the gradient

of the color field as introduced in [MCG03]. If the magnitude of n is larger than a threshold δ, particle i is closed to the fluid surface and we need to compute the surface tension force. Moreover, we also store these particles as surface particles in order to conveniently compute the thermal radiation gathered from the outside heat sources and the thermal dissipation to the contact medium from surface.

3.2.2

Viscoelastic Fluid Model

In order to compute the elastic term in Equation 3.7, our method follows the concept in [GBO04] and [CBL+09], which is to decompose the total strain into elastic and plastic component as in Equation 3.11),

T otal = P lastic+ Elastic (3.11) and compute each component by integrating their time derivatives. The total strain at a given time t can be formulated as in Equation 3.12.

T otal_t+∆t = T otal_t +d

T otal

dt ∆t (3.12) where ∆t is the time interval for each iteration. We can get the P lasticand Elasticby a similar form. Assuming the total strain and the plastic strain are both zero in the initial condition, we can evaluate the elastic strain rate by computing the difference of the total strain rate and plastic strain rate as in Equation 3.13:

dElastic dt = dT otal dt − dP lastic dt (3.13)

3.2 Fluid Simulation Framework 17

The total strain rate is given as follows:

dT otal dt ∆t =

∇v + (∇v)T

2 (3.14) where ∇v is the velocity tensor field defined on each particle and can be computed by Equation 3.15: ∇vi = X j mj ρj (vj − vi) ⊗ ∇W (rij, hi) (3.15)

We can calculate the plastic strain according to the Von Mises’s criterion: there is no plastic flow as long as the Frobenius norm of the elastic strain deviation 0 is still below to a yield threshold γ. Once the magnitude is greater than the yield point, the rate of the plastic flow is proportional to the difference of the Frobenius norm and γ.

0 = Elastic−T race( Elastic₎ 3 I (3.16) dP lastic dt = ξ 0 k0_kmax(0, k 0_{k − γ)} (3.17)

where 0 is the elastic strain deviation, T race(A) = a11+ a22+ ... + ann =Pn_i=1aii, I is the

identity matrix, γ is the yield threshold, and ξ is the material’s elastic decay rate. Finally, we can have the elastic strain rate from the Equation 3.13, and we can iterate the elastic strain tensor for a given time step by integrating the time derivatives. Since we get the elastic strain for each particle, the elastic stress force acting from particle pj on pi can be computed as follows:

f_ijElastic = mi ρi mj ρj κi+ κj 2  Elastic i · (∇W (rij, hi) + ∇W (rij, hj)) (3.18)

3.3 Heat Transfer System 18

3.3

Heat Transfer System

The melting process is triggered when the temperature reaches the melting point. Thus, we require a thermal transfer simulation for the heat energy flow in the fluid elements, which are SPH particles in our approach. We have modeled several types of thermal exchange types as listed here:

• Thermal radiation from heat source to surface particles,

• Thermal conduction from particle to its neighbor particles,

• Thermal dissipation from surface particles to cold air.

3.3 Heat Transfer System 19

Figure 3.4: Thermal simulation in our approach.

The thermal information is stored on each particle and heat energy can be exchanged with the neighborhood in the particle supporting range. Every particle has the heat properties as follows: total heat energy, temperature, absorb energy, latent heat threshold, conductive rate, emissive rate, and adopt heat energy coefficient. The total heat energy shows the current heat energy amount in this particle and can be converted into temperature using Equation 3.19:

Ti =

Qi

miCi

(3.19)

where Tiis the temperature, Qiis the heat energy stored in particle i, and Ciis the specific heat

3.3 Heat Transfer System 20

the elastic stress coefficients for the melting process from solid to liquid. When the temperature reaches the melting temperature, which invokes the phase change condition, particles start ab-sorbing energy to destruct the solid structure to become liquid. The viscosity and elastic stress coefficient are inversely proportional to the ratio of absorb energy and latent heat. Conductive rate defines the conduction between particles and emissive rate represents the emissivity from the surface particles to the air. Finally, adopt heat energy coefficient is used to adjust the amount of heat energy obtained from heat sources per iteration.

Besides, to simulate the thermal radiation and the thermal dissipation from the object sur-face, we have to identify the surface particles first. We make use of the gradient of the color field in Equation ??, which also represents the surface normal direction. By detecting whether the magnitude exceeds a given threshold, we may extract the surface particles and the surface normal as shown in Figure 3.5.

Figure 3.5: Surface particle and its normal.

Thermal Radiation transmits energy from heat sources to the surface particles in a radiative way. In this simulation framework, we only consider the direct radiative ray from the sources.

3.3 Heat Transfer System 21

Indirect radiation (multi-bounce of radiative ray) is ignored in our approach since its effect is insignificant. Therefore, we first detect the visibility by ray casting between surface particles and the heat sources. If the particle pass the visibility test, we compute the angle between the surface normal vector and the vector from particle to heat source position (see Figure 3.6). We than calculate the energy value depending the Equation 3.20.

Figure 3.6: Thermal radiation diagram

∆Q =X

s

δm

3.3 Heat Transfer System 22

where ∆Q is the energy variation in each iteration, δ is the adopt heat energy coefficient, d is the distance between the source and the particle position, cos θ is the angle we mentioned above, Qsis a function representing the heat energy value generated by the heat source. Qs(Ts, ∆t) ∝

Ts∆t in our implementation.

We have modeled two kinds of heat sources in our simulation. The first one is the point heat source as in the Figure 3.4. Another one is the area heat source, we modeled this as a plane which provides the heat energy for one side. The area heat source is used to model the heat source like a hot pan or hot ground.

Thermal Conduction diffuses the energy between particle and particle. According the physics of thermology, the heat energy flows from region of high temperature to region of lower tem-perature until it reaches the equilibrium state. Therefore, we model the heat conduction only if the particle temperature is higher than its neighborhood and the diffusion range is defined as the particle’s supporting range. The conduction is proportional not only to the ∆T , but also to the mass and the inverse of distance between particles, as shown in the following Equation 3.21 and 3.22: Qc= λiQi mj(Ti−Tj) dij Wsum (3.21) Wsum = X j mj(Ti− Tj) dij (3.22)

where Qcis the heat energy through conduction, dij is the distance between particle i and

particle j, and λ is the conductivity we have mentioned above. We distribute the energy flow fairly to the neighborhood particles. Equation 3.21 computes the heat energy conducts from par-ticle i to its neighborhood parpar-ticle j, and Equation 3.22 evaluates the weighting of conduction by considering the distance, particle mass, and the difference of temperature of two particles.

Thermal Dissipation is another thermal radiation in which the energy is radiated from the surface to the outside medium. The general equation of thermal radiation can be expressed as

3.3 Heat Transfer System 23

Equation 3.23 because most of the materials in nature are considered as grey bodies, we need to consider the emissivity . This factor has to be multiplied with the radiation spectrum formula before integration. For simplicity, we can simplify the factor as a constant similar to [MGG+10] did.

Q =  · σ · A · T4 (3.23) Equation 3.23 is the equation of the general thermal radiation in physics, where factor σ is the Stefan-Boltzmann constant,  is the emissivity, A is the radiating surface area, and T is the object temperature. Based on Equation 3.23, we formulate the thermal dissipation in our simulation as the following form:

Qemissive = σVp(Tp4− T 4

air) (3.24)

Equation 3.24 computes the heat energy dissipate from the surface to the contact medium. We use the particle volume Vp to approximate the area term since it is much easier to compute

in our approach.

After all the computation of energy transfer, we may update the temperature of particles through Equation 3.19, also update the viscosity according to the ratio of energy absorbtion and latent heat value.

3.3.1

Previous Thermal Diffusion Method

In the previous methods for object melting problem ([CMVHT02], [PPLT06], [CBL+09], [WLK03]), neither of them had implemented the thermal radiation effect, and some of them use the SPH interpolation scheme to simulate the thermal diffusion between particles. But in the real world, thermal radiation is actually a domination part in object melting. And the thermal conduction usually propagate the heat energy much slower than the radiation. In their demos, we found it is hard to capture the detail drops flowing and the melting behavior caused by the exposure to heat sources. Also, the thermal radiation allows us to determine the melting region and the melting

3.4 Adaptive Particle Sampling 24

speed in a natural way, but previous methods require users to define particle temperature in other rules.

Figure 3.7: Candle melting without thermal radiation and dissipation

Figure 3.7 shows the melting simulation involving only thermal conduction, which is done by using the SPH interpolation as in [CBL+_{09]. The temperature is increased by adding a}

constant value per iteration around the heat sources.

3.4

Adaptive Particle Sampling

The goal in adaptive particle sampling is to reduce the particle number in less important regions. Therefore, the computational resources could focus on the regions of interest; for example, the

3.4 Adaptive Particle Sampling 25

melting particles. Our adaptive scheme is basically based on the works in [APKG07], but we additionally take temperature influence into account. We adapt the particle sampling density during the simulation, and classify them into different levels l, where l = 0, 1, 2, 3.... The lowest level 0 particles have a mass m and support radius r. The particles in other levels have mass and support range computed according to their level as follows:

mi = 2lim (3.25)

ri =

3

√

2li_{r (3.26)}

The criterion of splitting or merging not only consider the ELFS value but also the tem-perature and the distance to the heat sources. After updating the ELFS and temtem-perature value, particles are added or removed depending on the conditions below:

A particle pi is split or merged if the following condition are all true:

• Split 1. Li > 1 2. ELF S(xi) < αri 3. Ti < α_LT_i or kpi− Hsk < α_Ld_i • Merge 1. ELF S(xi) > βri or (Ti < β_LT_i and Li < 3) 2. kpi− Hsk > β_Ld_i

where Li is the level of particle i, Ti is the temperature of particle i, kpi − Hsk represents

the distance between the particle and the heat source. The rest of six α, αT, αd, β, βT, βd

are different coefficients respectively. The split operation happened when the ELFS value is below to a certain threshold, and the temperature is increased over the value defined by the

3.4 Adaptive Particle Sampling 26

inverse level. To avoid the split and merge operations take place too frequently due to the dramatically change in temperature, we also consider the distance between particle and heat source since the region close to heat source may has higher possibility to vary the temperature. On the other hand, the particles merge if the ELFS value or the temperature threshold satisfy our condition. We limit the merge level to prevent the excessive merge, and this level upper bound is an experimental result in our simulation.

Figure 3.8 shows the particle distribution and the reconstruct surface in the adaptive and the general version. The first column is the general sampling condition, which uses 17685 particles during the simulation iterations. The second column is our adaptive sampling scheme, the par-ticle number can be lower at most to minimum 5994 parpar-ticles, achieving a great advantage of computation efficiency. Figure 3.8 (A) are the reconstructed polygon meshes using marching cube algorithm, although the particle numbers differ greatly, the surface is still similar to each other. Figure 3.8 (B) is the particle view for our simulation, and the color from red to blue rep-resent the ELFS value from low to high. Figure 3.8 (C) is the color-coded particle temperature, we aggressively merge the particles where the temperature is still low.

3.4 Adaptive Particle Sampling 27

Figure 3.8: Adaptive particle sampling (on the right) compared to uniform sampling (on the left).

3.4 Adaptive Particle Sampling 28

Our implementation of split and merge operation is a modified version of [APKG07]. Figure 3.9 (a) and (b) indicate the resampling operations. Figure 3.9 (a) shows how the split operation works. While a particle need to perform the split operation, it will split two particles place on the symmetric position of its original position. After the splitting, the original particle will be deleted. Figure 3.9 (b) indicates the merge operation. While two particles are all tagged to be merged, the new particle will be generated at the position of the centroid of two particles and the two merged particles will be deleted.

Figure 3.9: Particle resampling scheme in [APKG07].

Split When a particle p is split, the new generated particles p1, p2are positioned symmetrically

around p. The position of new particles should be chosen carefully to avoid the large pressure forces caused by the neighboring particles. Figure 3.9 (c) shows the valid region to sample the split particles. The possible valid region is defined as a sphere with radius d to the p’s center position. We check all the neighboring particles ,which the distance between p and the neighboring particles are less than 2d, in order to narrow down the valid regions which can place the new particles. The regions are eliminated around particle p if they are too close to the neighboring particles, which means the distance is lower than d(the red arcs in Figure 3.9 (c)). For the rest of regions become valid region if the two split particles can be placed symmetrically around the particle p position(the green arcs in Figure 3.9 (c)). If there is no valid region can

3.5 Flowing Control System 29

be found, the distance d will be lowered down iteratively until we have a valid region. After finding the region we can sample on, we may have many possible symmetric pair candidates, and the algorithm just randomly pick one of them.

As for the parameters of particles p1 and p2, their levels l1 and l2 are decided as li − 1.

The parameter about the particle’s properties (e.g. conductivity, emissivity, temperature, elastic stress coefficient, viscosity coefficient, and velocity ...etc) is the same as the original particle p, and the energy store in p (e.g. absorb energy, heat energy) can be distributed evenly to p1 and

p2 since their mass are equal.

Merge Two particles pi and pj can be merged if they are neighbor to each other and both of

them pass the conditions we have mentioned above. A new particle pkwill be generated at the

position xk = xi+xj

2 , and the level is lk = li + lj. The position should be chosen carefully to

inside the fluid volume and also avoid placing too close to other particles, which may cause large pressure forces. If there is no merge pair pass the above criteria, then the merge operation will be canceled in this iteration.

The parameters for the merged particle pk can be a weighted sum of that for pi and pj

according to the mass. By a weighted sum computation, we can correctly handle the merge operation for particles in different level or with different material.

3.5

Flowing Control System

So far, our method provides an reasonable fluid simulation framework and a thermal transfer scheme for simulating the heat flow, which can trigger the melting event. Although we have the ability to simulate the object melting, the animators may expect to design the flowing pattern for object melting, including the shape of the fluid. To this end, we develop a flowing control system, aiming to provide two types of control mechanism for different purpose. The first one is called flowing path control, which is used for guiding the direction of melting flow. The second one is named as melting shape control, which generate an approximate geometric potential field

3.5 Flowing Control System 30

to constrain the flow shape. We utilize the control particles concept as in [TKPR06] to guide the fluid volume.

Figure 3.10 shows the control system flowchart in our approach.

3.5 Flowing Control System 31

At the initial stage, we sample a set of control particles through the SPH particles. We define an influence range for each control particle to control the nearby particles. Each control particle has two forces to affect the movement of SPH particles. Attraction force drags the nearby particles to the control particles and the velocity force gives the nearby particles the forces based on the control particle’s moving direction. The strength of force is inversely proportional to the distance between the control particle and other SPH particles. We make use of these control particles to guide all the SPH particles move as the user designed. Furthermore, we also utilize the attraction force to enhance small droplet formation during the melting process.

3.5.1

Path Control

For path control, user can draw several flowing paths directly on the object surface or other constrainted surface. These path are used for guiding the nearby control particles to flow along the path. Instead of apply forces on all the normal SPH particles, we consider the control particles, which is a subset of SPH particles, are enough to guide the fluid direction. Moreover, by lowering the strength of force with respect to the distance to the center of control particle, more fluid detail can be preserved. Each path is defined by a set of path nodes generated along a user-draw trajectory, the path nodes have the information about the influence range R, and two forces direction: attraction force and velocity force, which can be used to drag the control particles close to the path and give the control particles a direction along the path. In the end, each path may give the forces to control particles that will move along the path trajectory.

3.5 Flowing Control System 32

(a) Attraction force (b) Velocity force

Figure 3.11: Two path control forces.

Attraction Force Our method compute the attraction force based on the vector Va and Vb,

which is shown in Figure 3.11(a). The control particle is attracted by the path along the direction −Vasin θ, the angle θ is the angle between Vaand Vb. Equation 3.27 is the final attraction force

for each control particle:

Fa=

X

i

(− ~Vasin θ)i (3.27)

where i indicates the index of each path, Fais the attraction force.

Velocity Force Velocity force provides the necessary force for fluid to flow along the pre-defined path. Figure 3.11 (b) represents how we generate the velocity force for each control particle. The path direction is decided by the vector from the current node to the next node position, and the magnitude of velocity force is proportional to the distance between the control particle and the path node. Equation 3.28 shows the formulation we use in our method:

3.5 Flowing Control System 33

Fv =

X

i

( ~VbW (d, R))i (3.28)

where i is the index of path, ~Vbis the direction along the path, W (d, R) is a smoothing function

(e.g. gaussian function) with influence range R defined on the path.

The final force for each control particle can be get as a weighted sum of Faand Fvas shown

in Equation 3.29.

Fpath = (Wa∗ Fa+ Wv ∗ Fv) (3.29)

In addition, we apply a feedback control scheme to adjust the force strength during simulation. Considering the fluid force acting on the control particles, which is computed by the fluid solver in 3.2, if this fluid force direction is away from the control direction, we need a larger control force to correct the flowing direction. The strength of path control is proportional to the angle θp between the force direction provided by the fluid solver and the force direction provided by

the paths. Equation 3.30 is the force provided by a path.

FP athf = Wp∗ Fpath (3.30)

Equation 3.31 shows the weighting coefficient Wp computation as a feedback scheme, where

kp is a constant coefficient and θp in Equation 3.32 is the angle between the force direction

provided by the fluid solver and the force direction provided by the paths.

Wp = kp(−θp+ 1.0) (3.31)

θp =

Ff luid· Fpath

kFf luidkkFpathk

(3.32)

where Ff luidis the force direction provided by the fluid solver and Fpath is the force direction

provided by the paths.

Except for the force control we have mentioned above, our path control also controls the temperature absorbtion rate to achieve our goal. Since the fluid will go toward low to flow, we slightly change the heat adoption coefficient δ at starting point of the path, so the melting may

3.5 Flowing Control System 34

happen faster than other regions, causing the melting flow toward our control path. The heat adoption coefficient is simply modified by a Gaussian function, making the adoption coefficient smoothly varying according to the distance to the center of path point. Figure 3.12 shows the adoption coefficients varies in a color-coded particles diagram.

Figure 3.12: Modify the adoption coefficient around the path starting point

3.5.2

Shape Control

It is quite often that artists have the needs to control the fluid to form a target shape. In our object melting simulation, we provide an easy-design shape control interface for artists to do these jobs. Differ from path control scheme, shape control must control all the SPH particles more precisely to obey the constraints we set. Because user will not expect there are small droplet cross the shape or there are some details do not follow this constraint. Therefore, instead of only applying on control particles, we enforce all the SPH particles in the influence range of control shape to accept the control forces.

As the same way in path creation, we allow user to draw shape contours on the scene surface, and store each contours with discrete nodes along the trajectory. Although the shape contour

3.5 Flowing Control System 35

is 2D on the surface, the influence is still on a 3D volume around each node. We think this is adequate for a melting process control since the fluid is usually flowing close to the surface. A shape contour has several shape nodes on the trajectory. And for each shape node, it stores a direction for the contour repulsion force and two influence range as shown in Figure 3.14. The farther range is the attraction range used to attract SPH particles to move close the contour line and the closer range is the repulsion range used to keep the particles stay inside the shape contour region.

Once a shape contour is drawn, we then compute the shape force direction for each node on the shape contour. By finding a nearest surface vertex normal vector and the vector to its next node position, we compute the shape force vector as the cross product of these two vectors, as shown in Figure 3.13.

Figure 3.13: Shape control force for each node.

We divide the influence range for each node into two level: attraction and repulsion range. Figure 3.14 reveals the two level range in our approach. Inside the light red region, we employ the shape force describe above to enforce the fluid inside the shape contour. In the light blue

3.5 Flowing Control System 36

region, we use a attraction force to attract the nearby particles to form our target shape, and the direction is simply the vector from particle to shape node.

Figure 3.14: Shape node two level influence range

Equation 3.33 shows the formulation of shape attraction force, which is used to attract the nearby SPH particles to form the target shape.

FSa= CSaVSaW (d − h1, h2− h1) (3.33)

where FSais the attraction force, CSais the coefficient of attraction strength, VSais the attraction

direction from particle position to the shape node position, and W is a smoothing kernel to smooth interpolate the magnitude of the force with the distance d.

The repulsion force computation can be done as Equation 3.34:

FSr = CSrVSrW (d, h1) (3.34)

where FSris the repulsion force, CSris the coefficient of repulsion strength, VSris the repulsion

direction for the shape node, and W is a smoothing kernel to smooth interpolate the magnitude of the force with the distance d.

3.6 Implementation Details 37

Moreover, we also control the temperature of particles within the influence region. We slightly lower down the temperature per iteration until it is below to the melting point, so the particle will stick to the target shape we design.

3.6

Implementation Details

For the particle moving, we integrate the particle velocity, position by using an explicit Euler scheme (Equation 3.35, 3.36) for a given timestep. Equation 3.35 shows the computation of velocity update in an explicit Euler scheme:

vt+∆t= vt+ ∆t ∗ (

pF

pm

) (3.35)

where vt+∆t is the velocity in next iteration, ∆t is the time interval, pF is the force acting on

particle, and pm is the mass of particle.

Equation 3.36 shows the computation for updating the particle position.

rt+∆t= rt+ ∆t ∗ pv (3.36)

where rt+∆tis the particle position in next iteration, pv is the velocity of the particle.

To speedup the neighborhood searching, we build a k-d tree per iteration to eliminate the farther particles quickly. The obstacle of the scene and the initial fluid volume are all represented as signed distance fields [FPRJ00]. Figure 3.15 shows four different models of the signed distance fields, we only draw the distance larger than zero. The color indicates the distance value to the object surface from black to red.

3.6 Implementation Details 38

Figure 3.15: Signed distance value in four different object

In the real world, some of the materials will change the phases quickly. For instance, the wax of the candle will directly burn out instead of change into the liquid phase. Our system also take into account of this circumstance, we randomly delete the high temperature particles to decrease the fluid generation rate. The possibility of deleting particles is depending on the material properties.

For the surface rendering, we employ the distance based surface tracking method to define the surface induced by particles [APKG07]. Then we extract a triangle mesh by using the marching cube algorithm [LC87]. In some part of our implementation code, we make use of the parallelization of multi-core CPU by applying OpenMP library. This could saves almost

3.6 Implementation Details 39

C H A P T E R

4

Results

In this chapter, we are going to present the experiment results for different simulation scenes. We list the particle number and the computational time for each iteration as a table below, and we also compare our candle melting result with the [PPLT09]. In the end, we show several melting scenes with different materials, and some of them are produced with the effect of control sketches painting by users. All the experiments are performed on Intel core 2 duo CPU E6750 2.66GHz with 4GB ram, and the video card is Nvidia GTX260.

4.1

Simulation Statistics

Our adaptive scheme will dynamically alter the particle number during the simulation, we sum-marize the related statistics in Table 4.1.

4.2 Simulation Results 41

Scene Min particle Avg min time(_iter.sec.) Max particle Avg max time(_iter.sec.) Candle 8017 1.577 16427 6.649 Ice Cream 3724 0.743 15399 7.574 Chocolate Bunny 4356 1.486 17685 11.338

Butter 2640 0.484 2640 0.484

Table 4.1: Simulation statistics of different scenes

Though the particle number in candle, ice cream, chocolate bunny are nearly the same, the average time are not close to each other. This is due to the number of particle in thermal simulation are different, the more particle need to simulate the heat energy transfer, the longer computation time it will cost. Table 4.2 lists the percentage of time consumption in the fluid solver and thermal transfer stages. It explains the reason in Table 4.1, which the particle number in thermal simulation stage plays an critical role in our computation cost.

Scene Fluid Solver Thermal Simulation Candle 25% 75%

Ice Cream 16% 84% Chocolate Bunny 12% 88% Butter 40% 60%

Table 4.2: The percent of time consuming in different stages

4.2

Simulation Results

In this section, we shows the simulation results in several parts, the results are displayed by a sequence of pictures. In subsection 4.2.2, we shows the heat energy propagation between particles by representing the temperature in color. It is clear that the material properties are changing during the heat up process. Subsection 4.2.3 shows the control effect create by users, including the path control and the shape control. And subsection 4.2.4, we compare our results

4.2 Simulation Results 42

to previous methods by simulating the same type of material melting. In the end, subsection 4.2.1 exhibit the rest of melting scene we have simulated.

4.2.1

Object Melting Scenes

We show two of the melting scenes we have simulated here: soft ice-cream melting and a chocolate bunny melting. In these two demo, we use two heat sources for giving heat energy. Figure 4.1 shows the melting process of the ice cream respectively. And Figure 4.2 displays the chocolate bunny melting situation.

4.2 Simulation Results 43

Figure 4.2: Chocolate bunny melting

The results here shows our method is capable to deal with different materials, and produce a reasonable melting behavior compare to the real scene.

4.2.2

Thermal Transfer Simulation Results

The following results shows the thermal simulation temperature in each particle, the color from blue to red represents the temperature from low to high. In the results, we can clearly see the heat energy propagate within the particles from surface particles to inside volume. And the particles close to heat sources may heat up faster than the particles are farther or be occluded.

4.2 Simulation Results 44

4.2 Simulation Results 45

Figure 4.4: Final image in butter scene

4.2.3

Fluid Control Results

Here we show the results of two types control forces in our method. In Figure ??, we make use of the path control force to guide the melting flow movement. We carefully set the control parameter to avoid the unreal fluid flowing and droplet spray. As in the result, by applying the path control, the melting fluid will intend to flow along the direction which provide by users.

4.2 Simulation Results 46

Figure 4.5: Candle melting in path control

The shape force control provides us the ability to design the overall shape of the melting fluid. Figure 4.7 shows a control example in butter scene, where user draws a smiling face as a target shape, and during the melting simulation time, the melting fluid will fill the target shape smoothly.

4.2 Simulation Results 47

Figure 4.6: The shape control line drawn by user

4.2 Simulation Results 48

4.2.4

Comparison with Previous Methods

Figure 4.8, 4.9, 4.10 shows a comparison of a candle melting scene between our result, Carlson et al. [CMVHT02], and Pavia et al.[PPLT09]. Through our thermal transfer works, the heat flow inside the particles can be simulated more naturally. Also, the attraction force of control particles slightly encourage the formation of detail droplet, make our result more realistic than the previous works.

4.2 Simulation Results 49

Figure 4.9: Candle melting in [PPLT09]

C H A P T E R

5

Conclusions

In this chapter, we brief summarize our object melting and flowing control approach algorithm. Furthermore, we analyze the limitations and purpose some future directions to improve our works.

5.1

Summary

In summary, we have purpose a particle-based fluid solver and a physics-inspired thermal trans-fer system to simulate the heat energy flow in fluid particles. By adding a elastic stress tensor to the traditional Navier-Stokes equation, our method can easily handle the high viscosity fluid behavior, which is often seen in melting scene. Also, our method take advantage of the adaptive particle sampling scheme to accelerate the simulation speed. In thermal transfer system, we have simulated three types heat energy transfer manner: thermal radiation, thermal conduction, and thermal dissipation. Through these transfer sways, the heat energy can be propagate to neighboring particles naturally, make us able to trigger the melting event happening based on the temperature in fluid elements.

5.2 Limitations and Future Work 51

Except for the physics-based simulation, we also pay attention to the controllability which is curial for animators. We utilize the control particles to guide the whole fluid volume flowing direction, which may preserve more fluid detail from the fluid solver. And we design two types of control mechanism, the path control can lead the fluid flow along the path, and the shape control has the ability to constraint the melting fluid to form our desire shape automatically. Besides the force control for particle movement, we additional control the temperature variables to achieve our goal during the melting process.

The results demonstrate that our method can melting several different materials naturally. And the animators can easily design the melting scene by the intuitive sketching tools.

5.2

Limitations and Future Work

Our object melting and flowing control framework is able to deal with a wide range of materials, and provide a realistic melting process. But there are still some limitations in our approach. For real objecting melting process, there are plenty of details (e.g. small droplet) flowing on the surface. Current object surface in our framework is reconstructed by marching cube algorithm, it is hard to avoid some small details are smoothing during the polygon creation. Moreover, since our control mechanism is a force-based system, it is difficult to control each particle to the the target position precisely. Finally, we have not try to simulate inhomogeneous materials in this thesis. The inhomogeneous material contains different material types in a single object, which makes the interaction between particles more complicate.

In the future, we would like to enhance the surface detail by applying more recently surface reconstruction method, such as adaptive marching cube or other advance techniques. Besides, currently our computational cost bottleneck is on the thermal simulation part, more particle number in thermal exchanging stage may lower down the simulation speed. To apply a impor-tant sampling scheme in thermal simulation maybe a good idea to improve the efficiency. The region of interest may use more samples to capture the accurate thermal transfer result, and the other region may use less particles for a sufficient simulation. Lastly, we are now focusing on

5.2 Limitations and Future Work 52

extending the path control works to form a 3D target. The path can be done by some mesh information like skeleton, and the nodes can store local volume inside the mesh and the contour around the path node. Hopefully, this control mechanism could help the melting fluid to form any given target shape.

Bibliography

[APKG07] B. Adams, M. Pauly, R. Keiser, and L. J. Guibas. Adaptively sampled particle fluids. In SIGGRAPH ’07: ACM SIGGRAPH 2007 papers, page 48, New York, NY, USA, 2007. ACM.

[CBL+09] Y. Chang, K. Bao, Y. Liu, J. Zhu, and E. Wu. A particle-based method for vis-coelastic fluids animation. In VRST ’09: Proceedings of the 16th ACM Sym-posium on Virtual Reality Software and Technology, pages 111–117, New York, NY, USA, 2009. ACM.

[CMVHT02] M. Carlson, P. J. Mucha, R. B. Van Horn, III, and G. Turk. Melting and flowing. In SCA ’02: Proceedings of the 2002 ACM SIGGRAPH/Eurographics symposium on Computer animation, pages 167–174, New York, NY, USA, 2002. ACM.

[DC96] M. Desbrun and M. Cani. Smoothed particles: A new paradigm for animating highly deformable bodies. In Ronan Boulic and G´erard H´egron, editors, Euro-graphics Workshop on Computer Animation and Simulation, EGCAS ’96, August, 1996, pages 61–76, Poitiers, France, 1996. Springer-Verlag. Published under the name Marie-Paule Gascuel.

[DG96] M. Desbrun and M. Gascuel. Smoothed particles: A new paradigm for ani-mating highly deformable bodies. In In Computer Animation and Simulation

Bibliography 54

96 (Proceedings of EG Workshop on Animation and Simulation, pages 61–76. Springer-Verlag, 1996.

[DKNY08] Y. Dobashi, K. Kusumoto, T. Nishita, and T. Yamamoto. Feedback control of cumuliform cloud formation based on computational fluid dynamics. ACM Trans. Graph., 27(3):1–8, 2008.

[FF01] N. Foster and R. Fedkiw. Practical animation of liquids. In SIGGRAPH ’01: Proceedings of the 28th annual conference on Computer graphics and interactive techniques, pages 23–30, New York, NY, USA, 2001. ACM.

[FL04] R. Fattal and D. Lischinski. Target-driven smoke animation. ACM Trans. Graph., 23(3):441–448, 2004.

[FPRJ00] S. F. Frisken, R. N. Perry, A. P. Rockwood, and T. R. Jones. Adaptively sampled distance fields: a general representation of shape for computer graphics. In SIG-GRAPH ’00: Proceedings of the 27th annual conference on Computer graphics and interactive techniques, pages 249–254, New York, NY, USA, 2000. ACM Press/Addison-Wesley Publishing Co.

[GBO04] T. G. Goktekin, A. W. Bargteil, and J. F. O’Brien. A method for animating viscoelastic fluids. ACM Trans. Graph., 23(3):463–468, 2004.

[HK04] J. M. Hong and C. H. Kim. Controlling fluid animation with geometric potential: Research articles. Comput. Animat. Virtual Worlds, 15(3-4):147–157, 2004.

[KAG+_{05] R. Keiser, B. Adams, D. Gasser, P. Bazzi, P. Dutre, and M. Gross. A unified}

lagrangian approach to solid-fluid animation. Proceedings Eurographics/IEEE VGTC Symposium Point-Based Graphics, 0:125–148, 2005.

[LAD08] T. Lenaerts, B. Adams, and P. Dutr´e. Porous flow in particle-based fluid simula-tions. ACM Trans. Graph., 27(3):1–8, 2008.

Bibliography 55

[LC87] W. E. Lorensen and H. E. Cline. Marching cubes: A high resolution 3d surface construction algorithm. SIGGRAPH Comput. Graph., 21(4):163–169, 1987.

[LD09] T. Lenaerts and P. Dutr´e. Mixing fluids and granular materials. Comput. Graph. Forum, 28(2):213–218, 2009.

[MCG03] M. M¨uller, D. Charypar, and M. Gross. Particle-based fluid simulation for interactive applications. In SCA ’03: Proceedings of the 2003 ACM SIG-GRAPH/Eurographics symposium on Computer animation, pages 154–159, Aire-la-Ville, Switzerland, Switzerland, 2003. Eurographics Association.

[MGG+10] N. Mar´echal, E. Gu´erin, E. Galin, S. M´erillou, and N. M´erillou. Heat transfer simulation for modeling realistic winter sceneries. Computer Graphics Forum, 29(2):449–458, 2010.

[Mon05] J. J. Monaghan. Smoothed particle hydrodynamics. Reports on Progress in Physics, 68(8):1703–1759, August 2005.

[MTPS04] A. McNamara, A. Treuille, Z. Popovi´c, and J. Stam. Fluid control using the adjoint method. ACM Trans. Graph., 23(3):449–456, 2004.

[PPLT06] A. Paiva, F. Petronetto, T. Lewiner, and G. Tavares. Particle-based non-newtonian fluid animation for melting objects. In Sibgrapi 2006 (XIX Brazilian Symposium on Computer Graphics and Image Processing), pages 78–85, Manaus, AM, oc-tober 2006. IEEE.

[PPLT09] A. Paiva, F. Petronetto, T. Lewiner, and G. Tavares. Particle-based viscoplastic fluid/solid simulation. Comput. Aided Des., 41(4):306–314, 2009.

[REN+_{04] N. Rasmussen, D. Enright, D. Nguyen, S. Marino, N. Sumner, W. Geiger,}

S. Hoon, and R. Fedkiw. Directable photorealistic liquids. In SCA ’04: Pro-ceedings of the 2004 ACM SIGGRAPH/Eurographics symposium on Computer

Bibliography 56

animation, pages 193–202, Aire-la-Ville, Switzerland, Switzerland, 2004. Euro-graphics Association.

[SP09] B. Solenthaler and R. Pajarola. Predictive-corrective incompressible sph. In SIGGRAPH ’09: ACM SIGGRAPH 2009 papers, pages 1–6, New York, NY, USA, 2009. ACM.

[SY05] L. Shi and Y. Yu. Taming liquids for rapidly changing targets. In SCA ’05: Pro-ceedings of the 2005 ACM SIGGRAPH/Eurographics symposium on Computer animation, pages 229–236, New York, NY, USA, 2005. ACM.

[TKPR06] N. Th¨urey, R. Keiser, M. Pauly, and U. R¨ude. Detail-preserving fluid control. In SCA ’06: Proceedings of the 2006 ACM SIGGRAPH/Eurographics symposium on Computer animation, pages 7–12, Aire-la-Ville, Switzerland, Switzerland, 2006. Eurographics Association.

[TMPS03] A. Treuille, A. McNamara, Z. Popovi´c, and J. Stam. Keyframe control of smoke simulations. ACM Trans. Graph., 22(3):716–723, 2003.

[WLK03] X. Wei, W. Li, and A. Kaufman. Melting and flowing of viscous volumes. In CASA ’03: Proceedings of the 16th International Conference on Computer Ani-mation and Social Agents (CASA 2003), page 54, Washington, DC, USA, 2003. IEEE Computer Society.