conduction, and thermal dissipation. Once the temperature of a particle reaches the melting point, the phase transition from solid to liquid takes place and the viscosity of melting liquid also varies according to the temperature. Moreover, to simulate the combustion, the liquid transits to the air when it reaches the flash point.
To reveal the realistic cohesion between melting liquids and solids and prevent liquid particles from penetrating the solid, we consider solid particles as ghost particles, and assign some elaborate attributes to them. We also compensate density deficiency near the free surface by sampling air particles that are viewed as ghost particles. The fluid particles thus can distribute evenly on the solid surface instead of being clustered to a shell. After solving Navier-Stokes equation using SPH, we cloud determine the velocity of each particle. However, in order to damp off the unphysical particle oscillation and allow the melting liquids flow smoothly on the solid surface, we further smooth the ve-locity of each liquid particle by applying XSPH to nearby particles. Besides, to simulate the thread formed by melting liquids, we also split the liquid particle dynamically.
We use screen space rendering approach to display the simulation result. To model the shape of viscous melting liquids flowing on the surface, we stretch the particle’s sphere to the ellipsoid. The stretchiness and the orientation of the ellipsoid are different according to the shape we want to model.
3.2 Fluid Simulation Framework
In this section, we first describe the framework of the fluid solver, and how we formulate the force exerted on the fluids by using SPH. Then, we explain the concept of using ghost particle in the object melting simulation.
3.2 Fluid Simulation Framework 13
3.2.1 SPH formulations and Navier-Stokes equation
SPH is an interpolation method defined on the particle system. With SPH, fluid at-tributes could be evaluated at discrete particle positions by using a symmetric radial function. According to SPH, a scalar quantity A at location r could be interpolated by a weighted sum of contributions from the neighbor particles:
A(r) =∑
j
A(rj)mj
ρj W (r− rj, h), (3.1) where j iterates over neighbor particles, mjis the mass of the particle j, rjis the location of the particle j, ρj is the density define on the particle j, A(rj) is the quantity at rj. The function W is a symmetric radial function with the kernel size h. We compute the density at r by using
which is also a derivation of Equation 3.1.
With SPH, the derivatives of field quantities only affect the smoothing kernel, so the gradient and Laplacian of A at location r could be derived easily by Equation 3.3 and 3.4 Fluids are often described by a velocity field v, density field ρ, and pressure field p.
The evolution of these quantities over time is governed by two equations. The first is mass conservation equation
∂ρ
∂t +∇ · (ρv) = 0. (3.5)
3.2 Fluid Simulation Framework 14
The second is Navier-Stokes equation which formulates the conservation of momentum ρ(∂v
∂t + v· ∇v) = −∇p + µ∇2v + ρg, (3.6) where µ is the viscosity of the fluid, and g is external force field. Because we use the particle system to model the entire fluid, both of these two equations could be further simplified. First, because the number of particles is constant and each particle has a constant mass during the whole simulation process, mass conservation is guaranteed naturally. Second, since the particles move with the fluid, the convective term v· ∇ v is not needed for the particle system. That is,
ρ(∂v
∂t) = −∇p + µ∇2v + ρg. (3.7)
There are three forces on the right hand side of the Equation 3.7 modeling pressure, viscosity, and external forces. We could formulate these force density terms acting on a particle i using SPH as follows [MCG03]:
fipressure =−∇p = −∑ To compute the pressure term in Equation 3.8 for each particle, we use ideal gas state equation
pi = k(ρi− ρo), (3.10)
where k determines the incompressibility of the fluid, ρi is the density defined on the particle i, and ρois the rest density of the fluid [DC96].
After computing all the forces exerted to the particle i, the acceleration of it can be computed by using Equation 3.11.
3.2 Fluid Simulation Framework 15
where ai is the acceleration of particle i. The velocity and the position of the particle i is then updated using Equation 3.12 and 3.13.
v∗i = vi+ ai∆t, (3.12)
where v∗i is the updated velocity and ∆t is the time interval, and
r∗i = ri+ v∗i∆t, (3.13)
where r∗i is the updated position.
3.2.2 Smooth flow of melting liquids
To reveal the cohesion phenomenon between solid and melting liquid and let the melting liquid flows smoothly on the solid surface, we apply the concept proposed by Schechter et al.[SB12] in our object melting simulation. First, by using ghost particles near fluid boundary, we prevent the density of a liquid particle from being under-estimated, which introduces obvious artifacts in SPH simulation. Second, by making the density field continuous across the solid boundary, we could reach no-penetration and no-separation solid-liquid boundary condition easily. Finally, we apply an artificial viscosity to damp off the non-physical oscillations and obtain smooth fluid behavior when the melting liquid flow on the solid surface.
Ghost particle Obviously, the fluid density distribution computed by Equation 3.2 is not even because the amount of neighbors between interior and boundary liquid particles differs. The pressure equation (Equation 3.10) leads to a case that particles near the boundary are clustered to a shell, aiming to rebalance the density; see the Figure 3.2 .
To correct the density summation for a fluid particle near the boundary of the fluid, we first consider solid particles as ghost particles with mass equal to the liquid particle and also contribute to the density summation for liquid particles. As a result, we pre-vent liquid particles near the solid boundary from being clustered to a shell due to the
3.2 Fluid Simulation Framework 16
Figure 3.2: The fluid particles clustered into a shell due to neighbor deficiency [SB12].
neighbor deficiency. Also, we add air ghost particles above the free surface of the fluid and within the kernel size [SB12]. However, for simplicity, the position of an air ghost particle for the liquid particle is set to be the reflection of some solid particles; as shown in Figure 3.3.
Figure 3.3: Sample air ghost particles.
Besides, to prevent liquid particles from either penetrating or leaving solid surface, the pressure field should be continuous through the boundary. Thus, after computing the density for each liquid particle, we set the density of each solid particle equal to the density of it’s nearest liquid particle. This step naturally enforces the no-penetration and no-separation boundary condition between solid and liquid.
3.2 Fluid Simulation Framework 17
XSPH To damp off the unphysical oscillation when the fluids flow on the solid sur-face, we apply the artificial viscosity to diffuse the velocity with nearby particles by Equation 3.14 :
vnewi = vi∗+ ϵ∑
j
mj
ρj (vj∗− vi∗)W (r− rj, h), (3.14) where j iterates over neighbor particle including solid particles, v∗i and v∗j are the ve-locities of particle i and j after solving the Navier-Stokes equation, ϵ is a parameter that represents the degree of diffusion [SB12].
However, before we apply the artificial viscosity, we set elaborate ghost velocity for solid particles. The velocity of liquids near solid boundary thus tends to converge to that ghost velocity after the computation of Equation 3.14 and achieves specific stickiness.
For the inviscid liquid, it often freely slips on the solid boundary because the normal component of liquid and solid velocities match and the motion of the liquid is only determined by its tangential velocity. Thus, Schechter et al. [SB12] constructed the ghost velocity vghostfor a solid particle by summing the normal component of particle’s true velocity vNsolidand the tangential component of the nearest liquid particle’s velocity vliquidT , as shown in Equation 3.15.
vghost = vNsolid+ vTliquid. (3.15)
Obviously, by constructing the ghost velocity for solid particle as Equation 3.15, the liquid near the solid boundary could achieve no-stick boundary condition because there does not exist relative-velocity in the normal direction and the velocity of the liquid does not suffer any friction in the tangential direction.
However, we modify the Equation 3.15 as the following:
vghost = vsolidN + αvliquidT , (3.16)
where α, 0≤ α ≤ 1, is a tunable tangential friction coefficient that allows us to achieve the desired fluid stickiness more easily. For the very high viscous liquid in our case, we let α approach to zero such that the liquid could almost stick on the solid surface.
3.2 Fluid Simulation Framework 18
3.2.3 Velocity Update
Figure 3.4 shows the detail steps for computing the velocity of the melting liquids. We first compute the density for each liquid particle by using Equation 3.2. In this step, both the solid ghost particles and the additionally sampled air ghost particles contribute to the density summation. Then, we set the density of each solid particle as the density of its nearest liquid particle, expecting that the pressure computed in the next step could be continuous across the boundary between solid and liquid. After computing the pressure using Equation 3.10 for each particle, we could compute the forces exerted on the liquid particle by using Equation 3.8 and Equation 3.9. Then, we solve the velocity for melting liquids by Equation 3.12. Finally, we set the ghost velocity for solid particle using Equation 3.16 and smooth the velocity of liquid particles by applying XSPH (Equation 3.14).
Figure 3.4: Flow for updating the velocity of melting liquids.