Simulation and Control Framework
3.2 Fluid Simulation Framework
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
Ajmj
ρ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
Ajmj
ρj ∇W (r − rj, h) (3.2)
∇2A(r) =X
j
Ajmj
ρj ∇2W (r − rj, 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 + µ∇2v + fsurf ace+ f + felastic (3.7) There are five force field in the right hand side of Equation 3.7. −∇p is the pressure term, µ∇2v is the viscosity, fsurf 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:
fijpressure = −∇p = −mi
3.2 Fluid Simulation Framework 16 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 otalt+∆t = T otalt +dT 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) 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:
fijElastic = mi The κi represents the elastic stress coefficient of particle i.