Patent Application
20150161810
The present disclosure relates generally to the field of computational dynamics, and, more specifically, to the field of computational fluid dynamics.
Fluids, such as liquids, smoke, fire, explosions and related phenomena, are responsible for many visually rich image, and simulating them has been an area of long-standing interest and challenges in computer graphics. In fluid simulation, enforcing incompressibility or constant density of a fluid object is crucial for realism and yet computationally expensive.
There are a variety of techniques available for fluid dynamics simulation, including grid based methods and particle based methods. Particle based methods are popular for their simplicity and flexibility. Smoothed particle hydrodynamics (SPH) is a particle based method for fluid simulation and has many attractive properties: mass-conservation; Lagrangian discretization (particularly useful in games where the simulation domain is often not known in advance); and conceptual simplicity. However, SPH is sensitive to density fluctuations from neighborhood deficiencies, and enforcing incompressibility is computationally costly due to the unstructured nature of the model. Recent work has improved efficiency by an order of magnitude, but unfortunately it requires small time steps, thereby limiting real-time applications.
For interactive application environments, computation robustness is a key issue; the simulation needs to handle degenerate situations gracefully. SPH algorithms often become unstable if the particles do not have enough neighbors for valid density estimates. A typical solution to avoid these situations is by taking sufficiently small time steps, or by using sufficiently many particles, at the cost of increased computation.
SPH can be used for interactive fluid simulation by incorporating viscosity and surface tension by using a low stiffness equation of state (EOS). However to maintain incompressibility, standard SPH or weakly compressible SPH (WCSPH) require stiff equations, resulting in large forces that limit the time-step size. Predictive-corrective incompressible SPH (PCISPH) relaxes this time-step restriction by using an iterative Jacobi-style method that accumulates pressure changes and applies forces until convergence. It has the advantage of not requiring a user-specified stiffness value and of amortizing the cost of neighbor-finding over many density corrections. Based on this framework, related work has been developed to solve incompressibility iteratively by setting up density constraints as a linear complementarity problem through a Gauss-Seidel iteration.
A hybrid approach that combines particle based and grid based methods has also been used to simulate fluid dynamics. For example, a grid is used to solve for pressure, and then velocity changes are transferred back to particles. A fluid implicit particle (FLIP) method can simulate incompressible flow with free surfaces. However, in an interactive setting, where the domain is not known as a priori, grid based methods are generally not applicable.
Position based dynamics (PBD) provides a simulation approach that enables direct control over positions of objects or vertices of a mesh. PBD has been previously used for simulating dynamics in electronic games based on Vertlet integration, where a system of non-linear constraints are solved using Gauss-Seidel iteration by updating particle positions directly. By eschewing forces and deriving momentum changes implicitly from the position updates, the typical instability associated with explicit methods can be avoided.
Therefore, it would be advantageous to provide a mechanism of simulating fluid dynamics while maintaining the incompressibility of a fluid object with high computational efficiency and stability. Embodiments of the present disclosure employ a computer implemented method of simulation of incompressible flow based on a position based dynamics (PBD) method. A set of constraints that enforce constant density of the particles in a fluid object are formulated in terms of neighbor particle positions. The formulated constraint equations can be solved iteratively in a Jacobi method to obtain a new position and new velocity of each particle in large time steps, as offered by the position based dynamics (PBD) framework, which makes the method suitable for real-time applications. Therefore, such a method can advantageously simulate incompressibility of a fluid object with computational convergence, and stability. Further, vorticity confinement may be introduced to simulate turbulent motion of the fluid object based on an unnormalized curl of the particle velocities. Moreover, a positive artificial pressure term can be incorporated in particle position updates to reduce particle clustering or clumping effects caused by negative pressure when a particle has too few neighbors.
In accordance with an embodiment of the present disclosure, a computer implemented method of fluid simulation of a fluid object comprises: (1) determining first positions of a plurality of particles based on external forces applied thereon; (2) identifying neighbor particles for a respective particle based on the first positions; and (3) determining a second position of the respective particle based on a constraint function representing that a density associated with the respective particle at the second position is approximately equal to a rest density of the fluid object, wherein the constraint function is a function of positions of the neighbor particles. The second position may be determined by: computing a position correction for the respective particle based on the constraint function; adding the position correction with a first position of the respective particle to generate the second position; and designating the second position as a new position of the respective particle.
The position correction may be proportional to a gradient of the constraint function. The position correction can be computed by iteratively computing the gradient, a scaling factor, and the position correction in accordance with a substantially Newton-Raphson method and a substantially Jacobi iteration method. The density may be defined based on a smoothing kernel function having a vanishing gradient at a kernel boundary, and further comprising pre-computing a corrective scale based on a reference particle configuration with a filled neighborhood in accordance with a predictive-corrective incompressibility smoothed-particle hydrodynamics (PCISPH) method. The method may further comprise: defining a positive artificial pressure in terms of the smoothing kernel function; and adding the positive artificial pressure to the scaling factor.
In another embodiment of present disclosure, a system comprises: a processing unit; and a non-transient computer-readable storage medium storing a fluid dynamics simulation program embodying instructions that cause the processing unit to perform: (a) accessing first positions of a plurality of particles in a fluid object based on an external force applied on the plurality of particles; (b) identifying neighbor particles for a respective particle based on the first positions; and (c) determining a position correction of the respective particle based on a constraint function representing that a density associated with the respective particle at the second position is approximately equal to a rest density of the fluid object, wherein the constraint function is a function of positions of the neighbor particles; and (d) determining a second position of the respective particle by adding the position correction to a first position of the respective particle; (e) designating the second position as a new position of the respective particle; and (f) repeating the (b) through (e) for all remaining particles of the plurality of particles.
In another embodiment of present disclosure, a computer implemented method of simulating fluid dynamics of an incompressible fluid object comprises: (1) determining first velocities and first positions of a plurality of particles in the incompressible fluid object based on an external force applied on the incompressible fluid object; (2) identifying a set of neighbor particles for a respective particle of the plurality of particles based on the first positions; (3) determining a second position of the respective particle based on a constraint that a density associated with the respective particle at the second position is approximately equal to a rest density of the incompressible fluid object; and (4) determining a second velocity of the respective particle based on a difference between a first position of the respective particle and the second position.
This summary contains, by necessity, simplifications, generalizations and omissions of detail; consequently, those skilled in the art will appreciate that the summary is illustrative only and is not intended to be in any way limiting. Other aspects, inventive features, and advantages of the present invention, as defined solely by the claims, will become apparent in the non-limiting detailed description set forth below.
Embodiments of the present invention will be better understood from a reading of the following detailed description, taken in conjunction with the accompanying drawing figures in which like reference characters designate like elements and in which:
Reference will now be made in detail to the preferred embodiments of the present invention, examples of which are illustrated in the accompanying drawings. While the invention will be described in conjunction with the preferred embodiments, it will be understood that they are not intended to limit the invention to these embodiments. On the contrary, the invention is intended to cover alternatives, modifications and equivalents, which may be included within the spirit and scope of the invention as defined by the appended claims. Furthermore, in the following detailed description of embodiments of the present invention, numerous specific details are set forth in order to provide a thorough understanding of the present invention. However, it will be recognized by one of ordinary skill in the art that the present invention may be practiced without these specific details. In other instances, well-known methods, procedures, components, and circuits have not been described in detail so as not to unnecessarily obscure aspects of the embodiments of the present invention. The drawings showing embodiments of the invention are semi-diagrammatic and not to scale and, particularly, some of the dimensions are for the clarity of presentation and are shown exaggerated in the drawing Figures. Similarly, although the views in the drawings for the ease of description generally show similar orientations, this depiction in the Figures is arbitrary for the most part. Generally, the invention can be operated in any orientation.
It should be borne in mind, however, that all of these and similar terms are to be associated with the appropriate physical quantities and are merely convenient labels applied to these quantities. Unless specifically stated otherwise as apparent from the following discussions, it is appreciated that throughout the present invention, discussions utilizing terms such as “processing” or “accessing” or “executing” or “storing” or “rendering” or the like, refer to the action and processes of a computer system, or similar electronic computing device, that manipulates and transforms data represented as physical (electronic) quantities within the computer system's registers and memories and other computer readable media into other data similarly represented as physical quantities within the computer system memories or registers or other such information storage, transmission or display devices. When a component appears in several embodiments, the use of the same reference numeral signifies that the component is the same component as illustrated in the original embodiment.
In general, position based dynamics (PBD) provides a mechanism for simulating dynamics by solving a system of non-linear constraints by updating particle positions directly. By eschewing forces and deriving momentum changes implicitly from the position updates, the typical instabilities associated with explicit methods can be avoided.
According to the present disclosure, incompressibility of a simulated fluid object can be achieved by projecting uniform density constraints on the particles in the fluid object. Such a constraint can be defined using an equation of state in which density associated with a respective particle is related directly to spatial distances between the respective particle and its neighbor particles. As a result, a new position and new velocity of a respective particle can be derived by solving the equation.
Then, based on a density constraint projected on the center particle, the estimated new position can be advantageously adjusted, or updated, at 103. The density constraint can be expressed as an equation of state that correlates a density associated with the center particle to positions of the identified neighbor particles. The new position of the center particle can thus be updated by iteratively solving the equation of state. At 104, the velocity of the center particle can be derived based on a difference between the current position and the updated new position.
It will be appreciated that the new position and velocity of the center particle can be further adjusted to achieve other simulated properties, effects, and/or phenomena related to the fluid object. For example, at 105, voracity confinement and viscosity can be optionally taken into account to derive a new velocity representing a spinning motion of the particle to achieve a turbulence motion of the fluid object for example. The viscosity may be defined in accordance with a smoothed particle hydrodynamics method (SPH), e.g., XSPH. In some embodiments, collision detection and response are also performed to further adjust the new position of the center particle, e.g., by projection-relevant constraints.
In some embodiments, a system of non-linear constraints can be solved to enforce constant density, with one constraint projected on each particle of a fluid object. Each constraint can be expressed as a function of the center particle's position and the positions of its neighbors, herein referred to collectively as p1, . . . , pn.
Consistent with a PBD framework, the density constraint on the ith particle can be defined using an equation of state:
where ρ0 is the rest density of the fluid object.
It will be appreciated that, the present disclosure is not limited to any specific constraint formula and/or equation that can enforce incompressibility or uniform density on a simulated fluid object. In some embodiments, ρi is given by the standard SPH density estimator:
where mj represents the mass of the jth neighbor particle, W is a smoothing kernel function corresponding to the contribution of a neighbor particle to the ρi and h is the corresponding smoothing length. Any suitable form of kernel function that are well known in the art can be used to correlate ρi to the particle positions p1, . . . , pn, such as a Gaussian function. In some embodiments, all particles can be treated as having equal mass. For purposes of simplification, the mass term is dropped from subsequent equations.
Provided with a set of estimated particle positions, p1, . . . , pn, that do not necessarily satisfy the density constraints, a particle position correction that satisfies the density constraint can be derived by solving
C(p+Δp)=0. (3)
This is found by a series of Newton steps along the constraint gradient:
Δp≈∇C(p)λ (4)
C(p+Δp)≈C(p)+∇CTΔp=0 (5)
C(p+Δp)≈C(p)+∇CT∇Cλ=0 (6)
Applying a recipe for the gradient of a function on the particles in accordance with a smoothed particle hydrodynamics (SPH) method, the gradient of the constraint function (1) with respect to a particle k is given by:
Placing (7) into (6) and solving for λ gives:
It will be appreciated that the present disclosure is not limited to any specific numerical or mathematical method of solving the constraint formula or equation to derive positions of the particles.
While “Iter<solver_iteration,” as at 201, a gradient of the constraint C and scaling factor s are calculated iteratively and independently for each particle at 202, for example based on equation (7) and (8). The calculation may be based on a Newton-Raphson method and any suitable iteration method, e.g., a Jacobi method or a Gauss-Seidel iteration method. In some embodiments, collision detection against an external object, e.g., a solid object, can also be included in the constraint solving loop. At 203, a position correction Δpi can be calculated for each particle, for example based on equation (7). At 204, the position correction Δpi is added with the estimated position pi to generate the new position that satisfies the density constraint.
Because the constraint function (1) can be non-linear, with a vanishing gradient at the smoothing kernel boundary, the denominator in equation (2) may cause instability when particles are close to separating. In some embodiments, this problem can be addressed by pre-computing a conservative corrective scale based on a reference particle configuration with a filled neighborhood, as used in a conventional predictive corrective incompressible SPH (PCISPH) method.
In some other embodiments, a constraint force mixing (CFM) method can be used to regularize a constraint, in which the constraint can be softened by mixing in some of the constraint force back into the constraint function, in the case of PBD, this changes (6) to
C(p+Δp)≈C(p)+∇CT∇Cλ+ελ=0 (9)
Where ε is a user specified relaxation parameter which can be a constant over the simulation. The scaling factor is now
and the total position update Δpi including corrections from neighbor particles density constraint (λj) is
Conventional SPH simulation methods typically suffer from the problem of particle clustering or clumping caused by negative pressure when a particle has only a few neighbors and is unable to satisfy the rest density. In embodiments of the present disclosure, an artificial pressure specified in terms of the smoothing kernel itself can be incorporated to avoid this problem. For example, the artificial pressure may be formulated as a correction of the scaling factor
where Δq is a point some fixed distance inside the smoothing kernel radius and k is a small positive constant. For instance, |Δq|=0.1h . . . 0.3h, k=0.1 and n=4. This term can then be included in the particle position update as
This purely repulsive term can be used to keep particle density slightly lower than the rest density. Consequently, particles pull their neighbors inwards causing surface tension-like effects. In some embodiments, the artificial pressure term is dependent on the spatial resolution and time-step. In some other embodiments, the artificial pressure, spatial resolution and time-step parameters are decoupled, which can advantageously offer anti-clustering effect independent from surface tension effects.
Conventional position based methods may introduce additional velocity or kinetic energy damping which is often undesirable. In some embodiments of the present disclosure, vorticity confinement can be used to replace the lost energy. This involves calculating the vorticity at a particle's location, e.g., by using an estimator
where vij=vj−vi. Then a corrective force can be derived by using the location vector
with η=∇|ω|i, e.g.,
f
i
vorticity=ε(n×ωi). (15)
In some embodiments, an unnormalized curl vector ω can be used such that vorticity only increases where it already exists. In some other embodiments, a normalized ω value is used such that viscosity is increased indiscriminately. Further, viscosity can be incorporated for purposes of achieving coherent motion. For example, the viscosity may be defined as typically used in a conventional XSPH method.
vi + Δltext(xi)
xi + Δlvj
xi* + Δpi
xi*
Table 1 provides pseudo-code describing an exemplary PBD simulation loop for simulating fluid object dynamics while maintaining incompressibility in accordance with an embodiment of the present disclosure. Lines 1-4 represent that an estimated velocity is first determined for each particle based on external force fext applied thereon and a time step Δt, and an estimated position is determined based on the estimated velocity. Lines 5-7 represent that a respective set of neighbor particles are identified for each particle based on the estimated positions resulted from Lines 1-4. The neighbor particles can be identified in any suitable method that is well known in the art.
Lines 8-19 represent the computation loop to iteratively compute the constraint error C, the scaling factor, and the position correction for each particle. The constraint error C, the scaling factor, and the position correction can be computed based on above-described equations (7) and (8) respectively for instance. Collision detection and response are also performed in the computation loop. Lines 16-19 represent that the position of each particle is updated by adding the position correction with the estimated position resulted from Lines 1-4. In some embodiments, signed distance fields stored as textures can be used for the particle-solid collision detection. Lines 20-24 represent that the particle velocities are updated based on the updated positions. In turn, particle velocity may be further adjusted by applying vorticity confinement and viscosity, e.g, based on equations (14)-(15).
In this example, distance and constraint values are recalculated in each solver iteration step, and particle neighborhoods are recomputed in each time step. This possibly leads to density underestimates, for example if a particle separates from the initial set of neighbors. In the conventional PCISPH approach, this can cause serious problems because once a particle becomes isolated, each iteration makes its pressure increasingly negative. If it then comes back into contact on a subsequent iteration, large erroneous pressure forces can be accumulated and applied. In contrast, the process according to the present disclosure can process the computation without being affected by such a situation because only the current positions, rather than accumulated pressure, are considered.
The data shown in
Table 2 summarizes the performance results of sample simulation processes in accordance with embodiment of the present disclosure. A frame time of 16 ms is used in all cases listed in Table 2. Table 3 summarizes a breakdown of a frame (percentages) for two examples that generated accordance with embodiment of the present disclosure. In the examples listed in Table 3, constraint solving includes collision handing processes.
A process according to the present disclosure, e.g., as described in Table 1, can be performed on any suitable processing unit, such as a central processing unit (CPU), or a graphic processing unit (GPU). Each stage of the algorithm is parallelizable, thus making the algorithm suitable for parallel architectures, such as a multi-processor GPU.
The present disclosure can be applied to simulate a fluid object including a liquid, smoke, fire, explosion and related phenomena.
Although certain preferred embodiments and methods have been disclosed herein, it will be apparent from the foregoing disclosure to those skilled in the art that variations and modifications of such embodiments and methods may be made without departing from the spirit and scope of the invention. It is intended that the invention shall be limited only to the extent required by the appended claims and the rules and principles of applicable law.
