Real-Time Smoke Simulation based on Vortex Method

In this paper, we use three different experimental methods (particle method, grid method and hybrid method) to model and simulate the smoke from the perspective of fluid dynamics. Through the comparison of different methods, we conclude: The particle method can avoid the numerical dissipation problem caused by grid calculation, but it also brings problems such as the distortion of the trajectory of the example. The grid method is accurate in calculation, but it is prone to numerical dissipation and loss of details. Finally, we choose the hybrid method to store the vorticity in the form of particles in vortex particles, avoiding the numerical dissipation problem caused by the use of grids, and including rich turbulence, which perfectly shows the simulation effect of smoke.


Introduction
Physics based smoke animation has been simulated for many years. In the early days, researchers used the Navier-Stocks equation to speed the motion in the graph. For example, Foster and Metaxes[1] introduced fluid simulation using the finite difference method into the graph for the first time, Stam [2] introduced a semi-Lagrangian advection scheme to achieve unconditional stability of the fluid, Fedkiw et al. [3] applied fluid simulation to similar smoke and added vortex control to enhance vortex detail. Based on this, Kim et al. [4] and Selle et al. [5] gave a higher order interpolation algorithm to improve its accuracy later. Zhu and Bridson [6] introduced FLIP to concentrate the advection term on the particles and discard the previous numerical dissipation.
Compared with other fluid animation, the simulation of smoke animation pays more attention to the display of surface vorticity. Later, some scholars introduced vortex particle methods, which can avoid the need to solve the elliptic equation separately in the velocitypressure equation. The vorticity-based fluid simulation mainly includes the brute force direct summation of the N-body problem of the particle, and a series of treecode [7], FMM [8], etc. These methods all show wonderful effects, but also bring problems such as flow field distortion. At the same time, at each time step, the new position of the N vortices is calculated, and O(N 2 ) operations are required. Due to limitations of computer resources and costs, it is impossible to make N too large. In order to improve the effect of vortex methods to simulate the fluid, a larger number of vortices are needed, and when the number of vortices is large, the chaotic state is prone to occur.
To overcome these shortcomings, we choose the Vortex-In-Cell method (VIC), which was first proposed by Christiansen [9]. We know that the vortex particle method is a no-grid Lagrange method. The VIC method is a hybrid method of semi-Lagrange and semi-Euler. In the numerical simulation of the vorticity field, the particle tracking Lagrange method is still maintained, but in the processing of the velocity field, we adopt the Euler method with grid. First, the velocity field is calculated by the Poisson equation of the flow function in the Euler frame, and then the vortex element is tracked along the trace in the Lagrange frame. This idea can be traced back to the MAC method of Harlow and Welch [10].

Related Work
In this paper, we introduce three methods, including particle method, grid method and hybrid method: Particle method: Selle et al. [11] used vortex particles to augment a more traditional Eulerian velocity-pressure simulation. Pfaff et al. [12] used vortex particles from preprocessed boundary layers to synthesize turbulence. Weiβmann and Pinkall [13] used Biot-Savart summation for vortex rings. Brochu et al. [14] captured vorticity with a Lagrangian mesh, including its generation from buoyancy, and used FMM for lineartime Biot-Savart summation. Lagrangian vorticity-based fluid solvers [15] evaluate the velocity field by summing over vortex sources with the Biot-Savart law. Berlin noise [16] etc.
Grid method: In fluid simulation, the common advection is the advection of velocity, such as Foster and Metaxas [1] used the Euler method to simulate the advection of velocity, and Stam [2] introduced the semi-Lagrangian method, by using linear backtracking, make it unconditionally stable. Fedkiw et al. [3] introduced the semi-Lagrangian method into the smoke simulation for solving the advection of temperature and density. Based on this, Kim et al., [4], Selle et al. [5] gave a higher order interpolation algorithm to improve its accuracy by using BFECC scheme and MacComark method respectively. With regard to the semi-Lagrangian method for vorticity advection, mainly for solving atmospheric models [17], contour [18,19], etc.
Hybrid method: Zhang and Bridson [20] used vortex segments to simulate vortex stretching. Marshall and Grant [21], Liu and Doorly [22] used the VIC method to directly calculate vortex stretching on a temporary Euler grid. Selle et al. [11] determined the stretching term by given velocity fields and grids. Zhang et al. [23] simplified the item and no longer used complex Jacobian matrices. Graham [24] first used the explicit Euler method to solve the diffusion term in the viscous solution 2D VIC method. Cheng et al. [25] used the "diffusion-vortex" scheme proposed by Lu and Ross [26] to solve the flow on a circular cylinder in a 2D VIC solver. In a system, if there are no other external factors, assuming that the temperature and density changes of the fluid element are represented by the material derivative DT/Dt and Dρ/Dt respectively, we can obtain the simplest equation Lagrangian form: Euler form: where k T and k ρ are non-negative diffusion coefficients.

Vortex equation
According to Boussinesq approximation [27]: A small change in density has no effect on the inertia and viscosity terms in the fluid motion equation, and the only effect is reflected in the buoyancy generated. Therefore, according to this kind of movement mechanism, combined with Navier-Stocks equation[1], we construct the fluid movement model as follows: Lagrangian form: where u is the ascending velocity of the fluid, ν is the viscosity coefficient, p is the pressure, and F bv is the buoyancy of the fluid. Then we find the curl (∇×) on both sides of the Equation ( Lagrangian form:

External forces
Compared with the surrounding air, the smoke fluid has the characteristics of high temperature and low density. Assuming that the unit volume is 1, the buoyancy of the smoke is: Where ρ 0 is the density of the air, ρ is the density of the smoke fluid, and g is the acceleration of gravity. At the same time, according to the thermal expansion performance of the fluid, the buoyancy of the smoke rising process can also be expressed as where temperature of smoke is T, the temperature of air is T 0 , g is the acceleration of gravity and β is the coefficient of thermal expansion. Therefore, the external force exerted by the smoke rising process is The vorticity is mainly generated by viscous forces and baroclinic fluids. Viscous force is mainly the interaction between fluid and solid, while baroclinic fluid forces is mainly caused by uneven density distribution between different fluids, as shown in Figure 1(a): In fluid A (smoke), the density is ρ, the pressure is p, the flow velocity is u upward due to buoyancy; In fluid B (air), the density is ρ 0 , the pressure is p 0 . Fluid A has a higher density and pressure due to its higher heating temperature, and the density ρ < ρ 0 , pressure p < p 0 . So the lower fluid B flows to the left while the upper portion starts to move to the right due to the increase of the density of fluid A, thereby causing the generation of vorticity [28]. At the same time, during the ascent, all vortices at the boundary of the smoke fluid and the surrounding fluid form a vortex ring (shown by the dotted blue line in Figure 1(b)). Meanwhile, Figure 1(b) shows the vorticity distribution produced by rising smoke. The value of vorticity is the curl of buoyancy of smoke rising (∇ × F bv ).

Application of semi-Lagrangian method in fluid field
The implicit semi-Lagrangian method can make the step size larger or the grid spacing smaller, which can achieve unconditional stability.  Take the advection term of the temperature field as an example, According to Equation (2), Similarly, to derive the density field and vorticity field by using the linear backtracking method, we can get Figure 2 shows a process in which the vorticity at time t and at position x can be linearly traced back to the vorticity at time t − ∆t and at position x − ∆tu n (x).
According to Equation (13), as long as the linear interpolation algorithm is used, there must be |T n+1 max | = |T n max |. The advection method of this item is unconditionally stable, that is, no matter how large ∆t is, or how fast u is, the final interpolation range must be within the range of the maximum value of the previous step. Therefore, the vorticity field will only decrease without increasing the external force.

Diffusion term
Meanwhile, in the diffusion term of the temperature field, for larger viscosity or smaller grid size, the implicit semi-Lagrangian method is also used, and the following form is changed: Then Where I is the identity matrix. Similarly, in the density field and vorticity field,

Induced velocity field of vorticity field
For the incompressible vector velocity field u, we define the stream function vector ψ, and u = ∇ × ψ. The relationship between vorticity ω, velocity u and stream function ψ is shown in Figure 3. Then By the limitation of ∇ · ω = 0, it can be assumed that the divergence of the stream function in the vorticity field is also 0, that is ∇ · ψ = 0, then The equation is a Poisson equation, which can be solved iteratively by the finite difference method. After ψ is obtained, the value of u can be obtained from equation When we use the particle method, in 3D space, according to the Green function method, we can get the stream function vector Where r is the distance and dV is the volume microelement, so u can be expressed as Preprints (www.preprints.org) | NOT PEER-REVIEWED | Posted: 31 May 2021 doi:10.20944/preprints202105.0762.v1 Then discretize it to get the velocity solution, namely Where r ij is the distance between the ith vortex particle and the jth vortex particle. Equation (24) and Equation (25) are the analytical and discrete forms of the Biot-Savart law, respectively. It can be seen from the Equation (25) that the velocity of each particle is obtained by summing the vorticity of all other particles in the space. This problem has the same form as solving the universal gravitation in the universe, which is also the N-Body problem.

Vortex-In-Cell (VIC) method
The Vortex-In-Cell (VIC) method is to use the Poisson Equation (21) of the Euler coordinate system flow function to calculate the velocity field, and then track the vortex particles along the trace in the Lagrange coordinate system.
First, we divide the grid on the computing plane. In order to solve the Poisson equation, we distribute the vorticity carried by vortices to the corresponding grid nodes. When we use the area weighting method, as shown in Figure 4.  Figure 4, it can be seen that the vorticity ω p of the vortex particles located in the grid is distributed to the surrounding grid points by means of area weighting (Figure 4 in 1, 2, 3, 4). The proportion of each node can be allocated according to the size of its diagonal area, respectively, S 1 , S 2 , S 3 , S 4 , extended to three dimensions, the vorticity on each node After the vorticity of all vortex particles is distributed to each node of the corresponding grid according to the above method, if n particles of one grid node are adjacent to each other, then the vorticity of each grid node ω is At the same time, the sum of the area weight W at this node is After obtaining the vorticity ω on the grid nodes, Equation (21) can be solved by the finite difference method (FDM), including Jacobian iteration, Gauss-Seidel iteration, conjugate gradient method, and multigrid method, etc. After obtaining ψ, the value of u can be obtained according to u = ∇ × ψ. The next step is to track the motion of vortex particles in the Lagrangian framework. According to the velocity on the grid node, the velocity of each vortex particle can be obtained by area weighted interpolation From this, the position of each vortex particle at the next moment can be calculated. The biggest advantage of the Vortex-In-Cell method is the high computational efficiency of the velocity field. Each time step only needs to calculate O(N) + O(Mlog 2 M) times, where M is the number of grids, which is one order of magnitude smaller than O(N 2 ) in the direct summation method. In this way, we can introduce a much larger number of vortex particles and improve the resolution of the flow field, especially in the area where the vorticity is concentrated, and can give full play to the advantages of Lagrangian and Euler methods. Because the grid vortex method is a method that combines grids and particles, compared with the pure finite difference grid method, this method is not constrained by the stability determination conditions when dealing with vorticity convection. That is to say, no matter what the time step is, it can be unconditionally stable.

Vorticity stretching
In 3D fluids, the interaction between the local vorticity vector and the velocity gradient results in vorticity stretching.
(1) If adopt the particle method, we can use the vortex segment method instead [20]. The vortex segment can be regarded as a cylinder rotating around the central axis. The direction of the central axis is consistent with the direction of the vortex. It has two endpoints x a and x b , length h, and curl κ. A vortex particle with a vortex intensity of ω i can have a length of h i and a unit direction oft = ω i / ω i 2 , curl κ = ω/(h iti ) represents a vortex segment. At the beginning of each time step, the vortex particles are converted into such vortex segments, and then the speed at both ends of the vortex segment is calculated, and the positions of both ends are moved.
When both ends of the vortex segment are updated, the vortex segment is stretched and converted back to vortex particles with vorticity It can be found that the process is cyclic. The whole process can be represented by Figure 5 as (2) If we adopt grid method, when the vortex element passes through the advection, it is simultaneously stretched by the velocity field, changing the vorticity field. Through the separation of the Euler equation, we need to solve From this, we can obtain a new vorticity In the calculation, if the grid scale is selected as h, then ε = h ω 2 . In this way, the velocity can be sampled at both ends of the grid to better show how the vorticity field is stretched by the velocity field.  It can be seen from Equation (33) that the velocity increment generated by the stretching term at a point in the vortex ring coincides with the direction of the vorticity of the point. As shown in Figure 6, the red arrow in the figure indicates the velocity increment caused by the stretching, which coincides with the direction of the vorticity at that point. At the same time, due to the stretching of the vorticity, the trajectory of the smoke flow is also bent (shown by the solid black line in Figure 6).

Implementation process
In this article, we have established three smoke models, using the particle method based on N-body summation, the semi-Lagrangian vortex method based on the grid, and the Vortex-In-Cell method based on the hybrid of particles and grids. The following are the flowcharts, algorithms and detailed steps of the three methods. In the three flowcharts, the round box represents the calculation based on particles, and the rectangular box represents grid-based computing. In order to satisfy the above vector field at the same time, we have select particles that contain temperature T p , density ρ p , vorticity ω p and tracking particles with velocity u p . At the same time, the variables on the grid include: temperature T, density ρ, vorticity ω, velocity u and flow function vector ψ.

Particle method (N-body method)
In the particle method, the model uses the N-body summation method based on the Biot-Sarvat law, which all calculations of temperature, density, and vorticity are based on particles. This method can avoid the numerical dissipation problem caused by grid calculation, but at the same time it brings about the distortion of particle trajectory and other problems. The flowchart is shown in Figure 7. The implementation process is shown in Algorithm 1: Algorithm 1 Particle Method 1: Emit(ω p , T p , ρ p , τ) 2:T p ,ρ p ← Diffuse(T p , ρ p , k, ∆t) 3: F bv = −(ρ 0 /ρ p )g − β(T p − T 0 )g 4: u p ← ω p 5:ω p ←Stretch(ω p , u p , ∆t) 6:ω p =ω p + ∇ × F bv 7: ω n+1 p ←Diffuse(ω p , k, ∆t) The specific steps are as follows: (1) The system emits a group of particles including vorticity ω p , temperature T p , density ρ p and particle lifetime τ.
(2) Temperature and density diffusion: It is very difficult to accurately simulate the viscous diffusion of particles. But a simple "trick" can be used to demonstrate particle diffusion [29]. Taking temperature as an example, according to the upward ejection of hydrothermal fluid, it can be assumed that the temperature of the smoke ejected from the nozzle is the highest, and when the smoke reaches the neutral buoyancy surface, the temperature is the lowest. Therefore, a function can be constructed that makes the temperature T p and the temperature decrease continuously over time, that is,T p = T p e −kt , where k is the coefficient that controls temperature diffusion. As shown in Figure 8, the temperature of each particle T p and at the initial moment of eruption, its temperature is the maximum value T p , and then as the particles continue to rise, their temperature continues to decrease. After a certain period of time, its temperature reaches its minimum value. In the same way, a density function can be constructed, that is,ρ p = ρ p e −kt , where k is the coefficient that controls the density diffusion. (3) According to the obtained new temperatureT p and densityρ p on the particle, the buoyancy on the grid can be obtained as (4) We use the Biot-Savart rule to perform the N-body summation. According to the particle velocity equation (25) obtained by analytic discretization, the velocity of any vortex particle can be obtained as Where r p is the position of the vortex particle, and r j is the position of any other particle. (5) According to the equation (31) that uses the particle method to solve the vortex of the stretching term, we convert the vortex particles into the position and velocity of the two ends, respectively, x p,a , x p,b and u p,a , u p,b , the new position of the vortex segment is After the vortex segment is updated, the vortex segment is stretched, and then converted to the vorticity of the new vortex particle as (6) According to the buoyancy of the particle, the new vorticity of the vortex particle can be obtained asω (7) Vortex diffusion: Similar to step (2), according to the process of smoke rising, it can be assumed that the vorticity of the first sprayed smoke is the largest. When the smoke rises to a certain position, its velocity and vorticity All tend to 0. Therefore, a function can be constructed that makes the vorticity ω p decrease continuously over time, that is, ω n+1 p =ω p e −kt , where k is the coefficient that controls the vorticity diffusion. (8) we use the third-order Runge-Kutta to find the position p n+1 of the next stage of tracking particles. First, get the particle velocity as Iterate this velocity to the next step to obtain the particle velocity Iterate again, get The position of the final particle is p n+1 = p n + 2 9 ∆tk 1 + 3 9 ∆tk 2 + 4 9 ∆tk 3 (42) (9) If the particle reaches its lifetime τ, remove the particle.

Grid method (semi-Lagrangian method)
This method uses the semi-Lagrangian method on the grid to calculate the convection and diffusion of temperature, density and vorticity. Using the grid method can make the calculation results more accurate, and can achieve unconditional stability, no matter how large the step size is. At the same time, it is not necessary to determine the stability condition of the step length. However, due to the existence of grid interpolation, numerical dissipation will occur. The flowchart is shown in Figure 9: The implementation steps are , Figure 9. Flow chart of smoke (grid method) shown in the algorithm 2: The specific steps are as follows: (1) The system emits a group of particles including vorticity ω p , temperature T p , density ρ p and particle lifetime τ. (2) Assign the particle temperature, density and vorticity to the corresponding grid nodes. Divide the grid in the calculation plane. According to the area weighting method (Figure 4), the vorticity ω, temperature T and density ρ at each grid node can be obtained.
(3) Solve the stream function vector ψ and velocity field u. First set the initial value ψ (0) = 0, and then solve the Poisson equation (21). Using the finite difference method in the grid, we can get: Two-dimensional: Three-dimensional: Most of vortex methods use fast and effective numerical solutions to elliptic equations. Using the iterative method, starting from an initial "guess" value ψ (0) , each step n produces an improved solution ψ (n) , gradually approaching the exact value. The superscript symbol indicates the number of iterations. Common solution methods include Jacobian iteration, Gauss-Seidel iteration, conjugate gradient, multigrid method [30], etc. In order to achieve rapid convergence, we adopt the multigrid method.
After obtaining ψ, the value of u can be obtained according to u = ∇ × ψ.
Where ψ 1 , ψ 2 , and ψ 3 are the three components of ψ. In 2D, because ψ 1 and ψ 2 are 0, then (4) This step uses the induced velocity u obtained in step (3) and combines the semi-Lagrangian method to advection temperature, density and vorticity. Taking the temperature solution as an example, Figure 10 shows the temperature field, the temperature on each grid (represented by pure black dots). To calculate the temperature at time t (marked with a red circle), the temperature field needs to be traced back to the red cross ×. Then we perform bilinear interpolation on the 4 temperature scalars (such as red box) at the grid closest to the red cross ×.
The new temperatureT can be solved by advection temperature. According to the Equation (13), the position of red × is defined as p, and any grid point is defined as p i,j,k , then p = p i,j,k − ∆tu/∆h, where ∆h is the grid size in this direction. Then according to the position of p, find the eight grid points adjacent to the position, and then calculate it based on bilinear interpolation (area weighting) (Figure 4) Similarly, according to Equation (14) and Equation (15), the density and vorticity are respectivelyρ The temperature diffusion is shown in Figure 11: Figure 11. Temperature diffusion According to Figure 11 and Equation (17), combined with the semi-Lagrangian method, which can be solved iteratively by referring to Equation (21).
Similarly, according to the Equation (18) and the Equation (19), the diffusion of the density field and the vortex field can be obtained respectivelỹ ρ i,j,k =ρ i,j,k + k ρ ∆t ∆h 2 ρ i−1,j,k +ρ i+1,j,k +ρ i,j−1,k +ρ i,j+1,k +ρ i,j,k−1 +ρ i,j,k+1 1 + 6k ρ ∆t ∆h 2 (53) ω i,j,k =ω i,j,k + ν∆t ∆h 2 ω i−1,j,k +ω i+1,j,k +ω i,j−1,k +ω i,j+1,k +ω i,j,k−1 +ω i,j,k+1 1 + 6ν∆t ∆h 2 (54) (6) The interaction between the local vorticity vector and the velocity gradient leads to vorticity stretching. In a 3D fluid, when the vortex element passes through advection, it is simultaneously stretched by the velocity field, changing the vorticity field. By separating Euler's equations, need to solve ∂ω ∂t = ω · ∇u (55) According to Equation (33), define two positions p1 = p ijk + ω∆t/(2∆h), p2 = p ijk − ω∆t/(2∆h), p1 and p2 are shown in figure 12 in the vorticity field: 1 p 2 p Figure 12. The position distribution of p1 and p2 in vorticity field in the stretch term Extended to three-dimensional grid coordinate representatioñ According to Figure 12, combined with the area weighting method, we can get The new vorticityω (7) According to the obtained temperature and density on the new grid, the external influence on the vorticity field can be obtained. Therefore, the new vorticitỹ (8) In this step, we need to return to the Lagrangian framework to obtain the motion of the vortex particles, and interpolate the vorticity on the grid to the vortex particles. According to Figure 4, we first divide the vorticityω(k) on the grid node by the assigned node weight W(k), Then use the area weighting method to sum the values on the obtained nodes to obtain the vorticity of the vortex particle At the same time, in order to describe the motion trajectory of the hydrothermal fluid, it is necessary to interpolate the velocity on the grid to the tracking particles, and then use the area weighting again according to the velocityũ(k) on the grid node that has been calculated Method interpolation to find the velocity of the trajectory particle (9) Same as step 8 of Algorithm 1. (10) If the particle reaches its lifetime τ, remove the particle.

Hybrid methods (Vortex-In-Cell method)
The processing of temperature and density is mainly on the grid, and the semi-Lagrangian method is used. As mentioned above, this method can achieve unconditional stability. No matter how large the step length is, it does not need to be stepped Stable condition judgment. At the same time, when dealing with vorticity advection, we use the Vortex-In-Cell method. This method saves the vorticity in the form of particles in the vortex particles, avoiding the numerical dissipation caused by the use of grids, and when dealing with induced velocity When the vortex is stretched, it is processed on the grid. This method is simple, efficient, and can fully show the details of fluid changes. The flowchart is shown in Figure 13. Based on this, the specific implementation process is shown in Algorithm 3: (1) The system emits a group of particles including vorticity ω p , temperature T p , density ρ p and particle lifetime τ.
(4) The advection of temperature and density is the same as step 4 of Algorithm 2.

Experiments and conclusions
This experiment uses three modeling methods (particle method, grid method and hybrid method) to simulate smoke. In this experiment, the three methods all use 1 600 vortex particles to provide vorticity and 409 600 tracking particles to track the trajectory; at the same time, in the grid method and the hybrid method, the grid scale is 32 × 64 × 32. The final real scene and the rendering scene of the three methods are shown in Figure 14. In Figure 14(a), although the particle method avoids the problem of numerical dissipation caused by grid calculation, the fluid trajectory is distorted and appears unreal. In Figure 14(b), the approximate outline of the smoke is drawn using the grid method, but due to the existence of numerical dissipation, numerical viscosity is generated, making the surface relatively smooth. In Figure 14(c), the vorticity is stored in the form of particles using a hybrid method, which avoids the problem of numerical dissipation caused by the grid, and when dealing with other variables (such as velocity, vortex stretching, etc.), the grid method is used, which fully shows the details of turbulence.