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.1 SPH formulations and Navier-Stokes equation
SPH is an interpolation method defined on the particle system. With SPH, fluid attributes 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), (1)
where j iterates over neighbor particles, mj is the mass of the particle j, rj is the location of the particle j, ρjis 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 for each particle i by using
ρ(r) =∑
which is also a derivation of Equation 1.
Besides, the gradient and Laplacian of A at location r could be derived easily by Equation 3 and 4
∇A(r) =∑
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. (5)
The second is Navier-Stokes equation which formulates the conservation of momentum ρ(∂v
∂t + v· ∇v) = −∇p + µ∇2v + ρg, (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. (7)
There are three forces on the right hand side of the Equation 7 modeling pressure, viscosity, and external forces. We could formulate these force density terms acting on a particle i using SPH as follows [MCG03]:
fpressurei =−∇p = −∑
To compute the pressure term in Equation 8 for each particle, we use ideal gas state equation
pi = k(ρi− ρo), (10)
where k determines the incompressibility of the fluid, ρiis 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 11.
ai = fi
ρi = fpressurei + fviscosityi + ρig
ρi , (11)
where ai is the acceleration of particle i. The velocity and the position of the particle i is then updated using Equation 12 and 13.
v∗i = vi+ ai∆t, (12)
where v∗i is the updated velocity and ∆t is the time interval, and
r∗i = ri+ v∗i∆t, (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 simula-tion. 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 The fluid density distribution computed by Equation 2 is not even because the amount of neigh-bors between interior and boundary liquid particles differs. The pressure equation (Equation 10) leads to a case that particles near the boundary are clustered to a shell, aiming to rebalance the density; see the Figure 7 .
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 prevent liquid particles near the solid boundary from being clustered to a shell due to the neighbor deficiency. Also, we add air ghost particles above the free surface of the fluid and within the kernel size [SB12]. 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 8.
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.
Figure 7: The fluid clusters into a shell due to neighbor deficiency [SB12].
Figure 8: Sample air ghost particles.
XSPH To damp off the unphysical oscillation when the fluids flow on the solid surface, we apply the artificial viscosity to diffuse the velocity with nearby particles by Equation 14 :
vnewi = v∗i + ϵ∑
j
mj
ρj (v∗j − v∗i)W (r− rj, h), (14) where j iterates over neighbor particle including solid particles, v∗i and v∗j are the velocities of particle i and j after solving the Navier-Stokes equation, ϵ is a parameter that represents the degree of diffusion [SB12].
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 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 vghost for a solid particle by summing the normal component of particle's true velocity vsolidN and the tangential component of the nearest liquid particle's velocity vliquidT , as shown in Equation 15.
vghost = vsolidN + vliquidT . (15)
By constructing the ghost velocity for solid particle as Equation 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.
For achieving the desired fluid stickiness more easily, we modify the Equation 15 as the following:
vghost = vsolidN + αvliquidT , (16)
where α, 0≤ α ≤ 1, is a tunable tangential friction coefficient. 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.3 Velocity Update
Figure 9 shows the detail steps for computing the velocity of the melting liquids. We first compute the density for each liquid particle by using Equation 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 10 for each particle, we could compute the forces on the liquid particle by using Equation 8 and Equation 9. Then, we solve the velocity for melting liquids by Equation 12. Finally, we set the ghost velocity for solid particle using Equation 16 and smooth the velocity of liquid particles by applying XSPH (Equation 14).
Figure 9: Flow of updating the velocity of melting liquids.