Patent Grant
7479963
This invention relates generally to computer graphic animation, and more particularly to computer graphic simulation of fluid flow.
The following non-patent references provide relevant background to the present invention and will be referenced in the discussion of prior art and description:
Animated natural phenomena, such as smoke, fire, and liquids, have become an indispensable component in contemporary computer generated animation, special effects, training simulators, and computer games. Currently, realistic and believable animations of these elements are most effectively achieved through physically-based numerical simulations of fluid dynamics. Efficient and stable methods for performing such simulations have been recently introduced to the computer graphics community [12,3,4,2,10]. However, providing animators with means for controlling such animations in an intuitive yet precise manner remains a challenging open problem. The present invention addresses this problem by introducing a new tool for controlling animations of smoke.
To illustrate just how challenging smoke animation control can be, consider the shot from the recent feature film The Fellowship of the Ring, where the smoke that Gandalf blows forms a beautiful galleon that sails right through Bilbo's expanding smoke ring. It would have been extremely difficult, if not impossible, for an animator to achieve such detailed and precise smoke formations and motion by manually adjusting wind fields and other smoke simulation parameters.
Many researchers in the field of computer graphics have used physically-based models and computational fluid dynamics (CFD) to simulate fluid flows. The interested reader can find good surveys in the previous work sections of [12], [3], and [13]. Different mechanisms for controlling fluid flows are disclosed in [6], [8], [13] and [15].
Foster and Metaxas [6] describe a mechanism referred to as “embedded controllers” via which an animator may introduce various effects into fluid animations by adapting boundary and fluid properties, and pressure and velocity fields. Thus, animators need not understand the fluid dynamics equations, or the low-level details of the solver. However, they must understand the dynamics of the effect they are interested in achieving.
Witting [15] describes the CFD tool that was used at DreamWorks Feature Animation to produce smoke and fluid effects for the feature film The Prince of Egypt. DreamWorks is a trademark of DreamWorks Animation LLC, California, USA. With this tool animators drive fluid simulations by using images and animated sequences to specify inputs such as initial temperature fields, heat sources, and forcing functions. Another production tool that provides animators both artistic and behavioral control for animation of flames is described by Lamorlette and Foster [8]. However, neither of these tools supports direct control over the desired results.
In a recent pioneering work Treuille et al. [13] introduced a method for keyframe control of smoke simulations. They employ a shooting technique to solve for the optimal wind forces that would cause a smoke simulation to approximate a set of user-specified keyframes as closely as possible. An animator controls a smoke simulation by specifying smoke density keyframes. Continuous quasi-Newton optimization is then used to solve for the control forces that would cause the simulated smoke to approximate the keyframes while minimizing the amount of force. The control force field is defined by a collection of parametric wind and vortex forces, and the optimization process searches for the appropriate parameter values.
While this approach has produced some very impressive results, it suffers from several drawbacks. First, the optimization framework is very expensive, as it requires multiple evaluations of the objective function and of its gradient. Each evaluation performs an entire simulation and, in addition to computing the velocity and the smoke density, their derivatives with respect to each of the control parameters are computed as well. Furthermore, the dimensionality of the optimization problem grows with the length of the simulation, so it is necessary to split the original problem into several smaller optimization problems, organized in two overlapping schedules, and to alternate between them. The dimensionality and cost concerns also dictate using a small number of control forces and a relatively coarse grid, making it difficult to specify and match highly detailed keyframes.
U.S. Pat. No. 6,266,071 (Stam) entitled “Method of producing fluid-like animations using a rapid and stable solver for the Navier-Stokes equations” and published Jul. 24, 2001 discloses a method for performing computer graphic simulation of a fluid in motion. The method takes input data and calculates the velocity of the fluid at a plurality of points at successive time intervals. The velocity values are sent to an animator module which produces a geometrical description of the scene. This geometrical description is then sent to a renderer module, which calculates an image using the geometrical description. The animation is then displayed on an output device. Scalar quantities such as temperature and density may also be calculated and sent to the renderer module, where they are used in calculating the image.
U.S. Pat. No. 5,999,194 (Brunelle) entitled “Texture controlled and color synthesized animation process” and published Jul. 12, 1999 discloses a process for producing an animation having fluid color, texture and consistency throughout the entire sequence. The process comprises the steps of creating key frames containing objects and characters having substantial color and texture and which correspond to an animated sequence. The key frames are digitized into a computer system for storage in a memory space. Two consecutive key frames are then defined as a source key frame and a target key frame. Corresponding features in both the source key frame and the target key frame are then outlined and in-between frames are generated by linear interpolation of each outlined Fig. in the source key frame and a corresponding outline in the target key frame.
The contents of both the above-referenced patents are incorporated herein by reference.
It would be desirable to provide a method and system for generating fluid-based animations, which allows complex simulations or animations to be controlled with very little additional cost compared to hitherto-proposed simulations and which avoids the need for optimization as required by Treuille et al. [13].
It is an object of the invention to provide a method and system for generating fluid-based animations, which allows complex simulations or animations to be controlled with very little additional cost compared to hitherto-proposed simulations and which avoids the need for optimization.
This object is realized in accordance with a first aspect of the invention by method for performing computer graphic simulation of an incompressible fluid in motion, the method comprising:
providing to a simulator module input data that includes a scalar field that defines a physical property of the fluid or of matter suspended therein and a sequence of target scalar fields that define a desired temporal evolution of the scalar field;
simulating fluid motion by calculating a velocity vector u of the fluid at a plurality of points in the simulated fluid, at a plurality of time intervals, wherein said calculation is performed by solving an equation:
ut=−u·∇u−∇p+f
where:
sending some or all of the scalar fields to a renderer module for producing a sequence of animation frames each relating to a respective scalar field; and
storing the sequence of animation frames for subsequent display.
According to a second aspect of the invention there is provided a system for performing computer graphic simulation of an incompressible fluid in motion, the system comprising:
a simulator module having an input for receiving data that includes a scalar field that defines a physical property of the fluid or of matter suspended therein and a sequence of target scalar fields that define a desired temporal evolution of the scalar field;
said simulator module being adapted to calculate a velocity vector u of the fluid at a plurality of points in the simulated fluid, at a plurality of time intervals, wherein said calculation is performed by solving an equation:
ut=−u·∇u−∇p+f
where:
an advection unit coupled to the simulator module for and being responsive to said velocity values at each time interval to advect the scalar field;
a renderer module coupled to the advection unit and being responsive to some or all of the scalar fields for producing a sequence of animation frames each relating to a respective scalar field; and
a storage unit coupled to the renderer module for storing the sequence of animation frames for subsequent display.
The invention provides improved animation control of fluid-borne scalar fields, using target-driven animation as the control paradigm. A particular application of the invention relates to simulation of a suspension of smoke particles in air as is done by Treuille et al. [13]. However, rather than attempting to optimally approximate a set of keyframes as done by Treuille et al. [13], the invention controls the simulation by a sequence of target smoke states, each serving in turn as an attractor for the simulated smoke. This is achieved by augmenting the standard fluid dynamics equations with closed-form terms designed specifically to drive the simulation from any given state towards a user-specified smoke density target.
The control paradigm according to the invention does not guarantee that each target is approximated in an optimal manner before switching to the next one, but it eliminates the need for expensive non-linear optimization. The present invention does not require keeping track of derivatives, and may thus use non-differentiable numerical schemes, such as flux limiter based hyperbolic solvers [9], which improve simulation accuracy. Thus, complex and interesting smoke animations, such as the one shown in
Specifically, the invention provides: (i) a new driving force term in the momentum equation, designed to induce the fluid flow that will carry the current scalar field towards a user-specified scalar field target; (ii) a new gathering term in the scalar field advection equation, designed to counteract diffusion of the scalar field due to numerical dissipation, thereby improving the ability of the simulation to match arbitrary targets; (iii) the ability for animators to independently control several scalar fields sharing a fluid.
It should be noted that the target-driven approach according to the invention does not provide direct control over the quality with which targets are matched. Instead, the animator is provided with several intuitive parameters for controlling the rate at which smoke evolves, as well as other characteristics of the simulation.
In order to understand the invention and to see how it may be carried out in practice, a preferred embodiment will now be described, by way of non-limiting example only, with reference to the accompanying drawings, in which:
a and 6b show different animation sequences defined by a dense sequence of similar targets where an animated 3D object is voxelized and the resulting volumes are used to drive a smoke animation;
The principles of the invention will now be described with regard to smoke simulation. It will be understood that smoke simulation requires the application of fluid dynamics to the flow of air in which there are suspended particles of smoke (carbon). However, it will be understood that the same principles are equally applicable to the suspension of other particles (e.g., dust, pigment) in air or other fluid carriers such as, for example, the flow of pigments through liquids. They are likewise applicable to the flow of temperature through a fluid and can be extended to model the propagation of heat and matter through the fluid, in addition to changes in velocity.
Let us denote by ρ=ρ(x,t) the time-dependent scalar field that specifies the density of smoke at location x and time t. Given an initial smoke density, ρ0=ρ(x, 0), and a sequence of target densities ρi*=ρ*(x,ti), our goal is to produce an animation in which the smoke is driven towards each of these targets in turn, while maintaining smoke-like dynamics and appearance. More specifically, during each time interval (ti−1,ti) the smoke should evolve towards the target ρi*. Note that unlike in traditional keyframe animation our problem statement does not require the i-th target to be matched precisely at time ti, but rather that at time ti the i-th target ceases to attract the smoke and ρi+1* becomes the attractor instead.
Realistic smoke animations are typically generated by numerically approximating either the Navier-Stokes or the Euler equations. The Navier-Stokes equations contain a viscosity term, not present in the inviscid Euler equations. These equations govern the mechanics of a medium fluid (thin air) in which smoke is present. Denoting the fluid velocity vector field by u=u(x, t) and its temporal derivative by ut, the inviscid Euler equations for incompressible flow are:
ut=−u·∇u−∇p+f (1)
∇·u=0 (2)
where p=p(x,t) is the hydrostatic pressure and f accounts for external forces affecting fluid flow, such as gravity and buoyancy. The term f may also account for other forces that are not strictly external such as viscosity, which as noted above is used by the Navier-Stokes equations. Therefore, for the purpose of the description and claims no distinction will be made whether the forces affecting fluid flow and encompassed by the force term f are external or intrinsic, thus allowing the Navier-Stokes and the Euler equations to be represented identically. By doing so, it will be readily understood that the invention may be applied to different fluid flow equations. Equation 1 is the balance of momentum equation (Newton's second law) for a fluid with unit density. Equation 2 is the continuity equation, which enforces conservation of mass. An additional equation describes the transport of smoke by the fluid flow (advection):
ρt=−u·∇ρ (3)
Here ρt is the temporal derivative of the smoke density field ρ.
The force term f in equation 1 provides an important means of control over the smoke simulation. For example, Stam [12] used this term to allow a user to control the flow by dragging a mouse. Treuille et al. [13] include in this term a set of parameterized wind and vortex forces. Their optimization process then searches for parameter values that cause the simulation to match the specified keyframes. Thus, in their approach f implicitly depends on the initial smoke density ρ0 and the target density ρ*. The approach according to the invention also utilizes the force term in order to drive the simulation from a given initial state to a target state, but this is achieved by introducing a driving force term that depends only on the instantaneous state of the system, rather than being defined by some global optimization condition; it is an explicit yet simple function F(ρ,ρ*).
In accordance with the invention, two modifications are made to the fluid flow equations presented above. First, as already mentioned, the force term f in equation 1 includes a driving force F(ρ,ρ*) and a momentum attenuation term −vdu is added:
ut=−u·∇u−∇p+f*+vfF(ρ,ρ*)−vdu (4)
In the above equation, the force vector f is replaced by a modified force vector f*, equivalent to the force vector f in equation 1, but without the driving force term F(ρ,ρ*) which is now shown discretely. The additional force vector F is referred to as the driving force term, since it is designed to exert forces on the fluid in such a manner that the resulting flow will carry the smoke to the prespecified target density ρ*. This term is derived below. The momentum attenuation term is added as a means for limiting the accumulation of momentum in the system.
Second, we add a smoke gathering term G(ρ, ρ*) to equation 3:
ρt=−u·∇ρ+vgG(ρ,ρ*) (5)
The purpose of this term is to counteract numerical dissipation that causes smoke to diffuse. It enables our simulations to closely approximate a target density field, even if it is “sharper” than the initial density field. The gathering term is derived below.
These new terms are controlled by the non-negative parameters vf, vd, and vg, whose effect is discussed below.
In this section we develop an appropriate expression for the driving force term F. Recall that the purpose of this force, at any given time t, is to cause the fluid to advect the current smoke density ρ(x,t) towards the next target ρ*(x). Furthermore, we must also take care to define F in a manner that enables the fluid to reach a rest state, u=ut=0, once the target has been reached. Failure to enforce this requirement will cause undesirable fluctuations even in regions where ρ=ρ*.
The rest state requirement may be satisfied by ensuring that the amount of momentum generated by F decreases as p approaches ρ*, allowing the momentum attenuation term to bring the system to u=u, =0 when ρ=ρ*. In this state, equation 4 reduces to:
F(ρ*,ρ*)=∇p.
This means that rather than forcing F to vanish, we can simply allow the hydrostatic pressure to cancel it out. Thus, it suffices to ensure that F(ρ*,ρ*) is a gradient of some potential field, rather than an arbitrary vector field.
Note that the gradient of the target density ∇ρ* always points uphill towards higher concentrations of ρ*. This means that the gradient is the desired direction of flow, and to induce it we align the driving force at any point x with ∇ρ*(x):
F(ρ,ρ*)∝∇ρ*
In practice, however, ρ* may be constant in some regions of the domain, causing ∇ρ*=0 there. In order to avoid this problem we use a blurred version of ρ*, which we denote {tilde over (ρ)}*. In order to ensure that ∇{tilde over (ρ)}*≠0 everywhere, the blurring filter must have sufficiently large support. The filter should also fall off rapidly in order to avoid unnecessarily smoothing the directions of the gradients. In our implementation we obtain {tilde over (ρ)}* by convolving ρ* with a Gaussian kernel given by:
g(x)=e−x
with a small value for σ. Although the gradient ∇{tilde over (ρ)}* is non-zero everywhere, its magnitude is very small away from regions where the target is defined, because of the rapid falloff of the blurring filter. Therefore, we use a “normalized” gradient instead:
Note that we are guaranteed that {tilde over (ρ)}*>0, since it is obtained by convolving the non-negative function ρ* with a Gaussian. As we shall see shortly, this form of normalization enables the driving force to comply with the rest state requirement.
In order to avoid applying forces unnecessarily, the magnitude of the driving force should be roughly proportional to the actual smoke density ρ. This could have been achieved by setting
However, this expression does not comply with the rest state requirement. In order to fix this it suffices to replace ρ with its blurred version {tilde over (ρ)}:
Note that for ρ=ρ* this definition gives us:
and thus F(ρ*,ρ*) is indeed a gradient of a potential field.
In general, it is not possible to match a given target density field solely by advection. Consider, for example, a case where the target field ρ* contains higher values and higher gradients than those present in the initial field ρ0. In the course of a simulation these initial values and gradients only decrease because of numerical dissipation, and therefore will never achieve the target values. Even if we were to somehow avoid numerical dissipation, it would still be impossible to match such a target: from the mathematical standpoint, due to incompressibility of the fluid, any advection applied to ρ0 is merely a spatial remapping of this density distribution, incapable of generating any new density values. Tackling this issue by introducing an inverse diffusion operator, such as −α∇2ρ, is physically unstable [11]. Instead, we introduce an alternative mechanism, which we refer to as smoke gathering. This mechanism is enforced by the G(ρ,ρ*) term in equation 5.
Let us assume that the total mass of the simulated smoke (the integral of ρ over the entire domain) is equal to that of the target smoke. In practice, this assumption can be realized by scaling the source or the target to make sure that the two masses are equal to each other. Now consider the time varying “density error” field e(x,t)=ρ(x,t)−ρ*(x). Obviously, any large differences in the values and/or gradients between ρ and ρ* give rise to high values and high gradients in e. The essence of our idea is to apply diffusion to e gradually smoothing it as time goes by. This is a stable process that will eventually make e constant across the entire domain, which implies ρ=ρ*, because of mass conservation.
Formally, we would like the density error e to satisfy the following diffusion equation:
et=∇2e (8)
with the von Neumann boundary conditions
which essentially state that the equation must be resolved inside the domain, as no flux of e across the boundary is allowed. As time advances, e becomes zero across the domain, which implies that
e=ρ−ρ*=0 and ∇e=∇ρ−∇ρ*=0.
Thus, diffusing the error enables ρ to evolve exactly into ρ*, given enough time.
Substituting the definition of e into equation 8 and noting that ρ* is constant with respect to time, we obtain the following equation in ρ:
ρt=∇2(ρ−ρ*). (9)
Thus, in order to achieve diffusion of the error, we could simply set G(ρ,ρ*)=∇2(ρ−ρ*) in equation 5. In practice, however, we'd like the gathering effect to occur only in the close vicinity of the target ρ*, otherwise we end up simply diffusing ρ:
ρt=∇2(ρ−ρ*)≈∇2ρ.
In addition, gathering should only be allowed to occur where some smoke ρ is actually present, otherwise we might obtain negative values of ρ. To account for these requirements we rewrite equation 9 as:
ρt=∇·∇(ρ−ρ*).
where ∇(ρ−ρ*) is the error flux, and modulate the error flux by both {tilde over (ρ)}* and ρ. Thus, our gathering term becomes
G(ρ,ρ*)=∇·[ρ{tilde over (ρ)}*∇(ρ−ρ*)]. (10)
Our driving force term attracts any smoke present in the system towards the nearest “peak” of the target density. To provide the animator with more precise control over the affinity between the smoke and the target it is sometimes useful to split the smoke in the system to a number of independent density fields, each with its own set of initial and target states. Although each smoke field is controlled independently, there is an implicit interaction between them, as they are advected by, and exert driving forces on, the same fluid. Given n smoke fields ρ1, . . . , ρn and their corresponding targets ρ*,1, . . . , ρ*,n the only change in equation 4 is that the driving force is now the sum of the n driving forces F(ρi,ρ*,i) and we now have n advection equations for the smoke fields, instead of one:
ρti=−u·∇ρi+vgG(ρi,ρ*,i)
Note, however, that there is still only one velocity field u.
The simulations produced by our methods are controlled via four parameters: the force smoothing parameter σ (eq. 6), the force coefficient vf, the momentum attenuation coefficient vd (eq. 4), and the gathering coefficient vg (eq. 10):
Given an initial velocity field u0, initial smoke density ρ0, and a target smoke density ρ* we simulate the fluid flow and the resulting motion of smoke by numerically solving equations 2, 4, and 5 over a series of discrete time steps.
Equations 4 and 5 express the rate of change of u and ρ as a sum of terms, each representing a different contribution. Such equations are often solved by the method of fractional steps, also known as operator splitting [11]. The same approach is used by Stam [12] and by Fedkiw et al. [3]. The main idea behind this method is to solve the original equation by integrating a succession of simpler equations, each with a single term on its right hand side. Each of the simpler equations can be treated more effectively by choosing the most appropriate method for integrating it. Specifically, we perform the following sequence of fractional steps in each time step:
In our implementation, steps 3, 5 and 6 are restricted by a CFL condition resulting from the explicit hyperbolic scheme used at these steps (described in appendix A). Formally, step 2 imposes another restriction on the maximal time step to be less than 1/vd. In practice, this is not a limitation as the typical values of vd are quite small.
We use a regular 3D grid to solve our equations numerically using finite differences. We use a staggered arrangement of the different variables [7]: the velocity variables are defined at the centers of the cell faces, while the pressure and density variables are defined at the center of each cell (see
and
respectively, while the density variables are denoted ρi,j,k. Thus,
is located at ((i+½)Δx,jΔy,kΔz), where Δx, Δy, and Δz are the cell spacings along the x, y, and z axes, respectively. When necessary, we denote the number of the fractional step using a parenthesized superscript.
For compactness and readability we introduce notation for several discrete operators. First derivatives are approximated halfway between the positions of the corresponding variable. For example, the first derivative of the velocity along the x axis is approximated as
and the first derivative of the density along the y axis is
Second derivatives are approximated at the positions where the variables are defined. For example,
is the second derivative of the velocity along the z axis.
Occasionally, we require the value of a field halfway between the locations where its variables are defined. In such cases, we approximate the value by using simple averaging. For example,
is the velocity along the y axis at the center of a cell.
Components of the driving force must be evaluated for each updated velocity variable at its corresponding location. We denote these components
and
Each of these components is evaluated by approximating equation 7 at the corresponding grid point, using a discrete derivative operator and averaging:
and similarly for Fv and Fw. The velocity variables are then updated as follows:
Note that when ρ=ρ*, the driving force is “projected out” by the projection step in the discrete sense, allowing the discrete velocity variables to be at rest. This is due to the fact that in this case the discrete force components reduce to
that is, the gradient of a discrete potential function, which is exactly what is being eliminated by the projection step.
The advection term appearing in fractional steps 2 and 5 can be written in conservation form
qt=−u·∇q=−(uq)x−(vq)y−(wq)z
where q stands for either u or ρ. We solve this equation as three separate 1D hyperbolic equations
qt=−(uq)x
qt=−(vq)y
qt=−(wq)z
using the high resolution hyperbolic solver described in Appendix A. This is a second-order numerical scheme that is considerably less prone to numerical dissipation compared to first-order schemes. Thus, for a given grid resolution more turbulent flows are obtained, and the fine features of the evolving smoke are better preserved. These advantages are demonstrated in
The gathering term defined by equation 10 is already given in conservation form, and is also solved by the same solver, independently for each axis:
ρt=[ρ{tilde over (ρ)}*(ρ−ρ*)x]x
ρt=[ρ{tilde over (ρ)}*(ρ−ρ*)y]y
ρt=[ρ{tilde over (ρ)}*(ρ−ρ*)z]z
The discrete velocity field (U,V,W)(3) after fractional step 3, is not divergence free. To impose equation 2 discretely, we use Chorin's projection technique [1]. We recover a divergence free field (U,V,W)(4) in the discrete sense, i.e.
(Dx(U(4))+Dy(V)(4)+Dz(W)(4))i,j,k=0,
by subtracting the gradient of the pressure field P:
Combining the above equations, we get the following set of equations for the pressure variables Pi,j,k
(Dxx(P)+Dyy(P)+Dzz(P))i,j,k=(Dx(U(3))+Dy(V(3))+Dz(W(3)))i,j,k, (11)
which is a discretization of the Poisson equation ∇2p=∇·u. This is a large sparse set of linear equations, which we solve using a standard multigrid solver [14].
The values of the velocity variables on the boundary are prescribed by the boundary conditions, so there we have (U,V,W)(4)=(U,V,W)(3). Consequently, for the left boundary cells, equation 11 reads:
This is also the case for all the other boundaries.
The invention has been implemented in 2D and in 3D using the C++ programming language. All of the timings reported below were measured on a 2.4 GHz Pentium IV with 1 GB of RAM running Linux. For 2D simulations on a 2562 grid a single time step takes 0.25 seconds. For 3D simulations on a 1283 grid each time step takes 10.6 seconds. This is only about 15 percent slower than the computation time for an ordinary smoke simulation (without our driving force and gathering terms) using the same implementation.
In all of the following examples between 5 and 15 adaptive timesteps were used for each animation frame, so each second of a 20 fps animation took about 50 seconds of simulation time in 2D, and about 35 minutes in 3D. This is faster by at least two orders of magnitude than the running times in [13], where they report that it took 2-5 hours to optimize a 50×50 2D simulation.
The invention has been applied to two different kinds of smoke animations. The first uses target states that are very different from each other, in essence producing a morphing sequence between pairs of states. Examples are shown in
Another kind of smoke animation is that defined by a dense sequence of similar targets. The idea is to voxelize an animated 3D object and use the resulting volumes to drive a smoke animation. Frames from two such animations are shown in
As discussed above, the invention allows animators to control several smoke fields independently.
Finally,
Smoke diffusion is unavoidable in numerical simulations. As pointed out above, it presents a serious problem in the context of controlled smoke animations. The gathering term we have introduced does provide a partial remedy to this problem, as it enables the simulation to converge to an arbitrary target. However, the resulting smoke transition does not always come across as a natural evolution of smoke, since it is based on a diffusion process and not on kinetics. For example, the animation might sometimes look like the target is simply emerging from an amorphous static cloud of smoke. Thus, the gathering mechanism should be employed to the least extent necessary, leaving the task of forming the target, as much as possible, to the advection induced by the driving force.
As pointed out earlier, the target-driven paradigm does not guarantee optimal approximation of the targets, and the animator cannot directly control how well a particular target is approximated at a specific instant in time. However, the animator may control the simulation progress in a less direct manner via the vf parameter. Finally, just as in traditional keyframe animation, skill and experience are still required in order to choose the right target states and the appropriate control parameters to achieve a good animation.
The invention presents a new method and system for direct control of smoke animation by introducing new terms into the standard flow equations. Each of these terms has a simple closed form, and consequently they allow complex smoke animations to be controlled with little additional expense on top of the cost of an ordinary simulation.
Similarly to previous methods for simulation of smoke [3] the simulator according to the invention is easily extended to support the usual external forces (gravity, buoyancy, etc.) in addition to the driving force. It is also possible to add obstacles represented by internal boundary conditions.
It will be appreciated that the gathering term acts as a counter mechanism to the diffusion mechanism which is inherent to fluid flow. The invention contemplates other terms for driving a scalar object, such as smoke, through a sequence of target states, as well as other anti-diffusion mechanisms. A multi-resolution version of the gathering term may be employed to make gathering faster, more responsive, and more natural looking.
Although in the preferred embodiment, the target scalar fields define a desired temporal evolution of the scalar field, the invention also contemplates additional mechanisms for controlling smoke animations. For example, the animator may want to sketch out preferred paths for the smoke to flow along, without the need to generate many intermediate keyframes between the starting and the ending positions.
Referring now to
The system 10 may in practice be realized by a suitably programmed computer having a graphics card or workstation and suitable peripherals, all as are well known in the art. A detailed description of a suitable system is given in above-referenced U.S. Pat. No. 6,266,071 whose contents are incorporated herein by reference and to which reference should be made.
For the sake of completeness,
The velocity values at each time interval are used by the advection unit 14 to advect the scalar fields and some or all of the scalar fields are sent to the renderer module 15 for producing a sequence of animation frames each relating to a respective scalar field, which is then stored in the storage unit 16 for subsequent display by the display unit 17.
As noted above, the system 10 and the simulator module 11 may be suitably programmed computers. Likewise, the invention contemplates a computer program being readable by a computer for executing the method of the invention. The invention further contemplates a machine-readable memory tangibly embodying a program of instructions executable by the machine for executing the method of the invention.
Likewise, the invention contemplates an animation sequence that is generated according to the method of the invention. The invention further contemplates a machine-readable memory tangibly embodying data representative of an animation sequence that is generated according to the method of the invention.
The solution for hyperbolic conservation laws of the form
qt+(H(q))x=0
can be approximated numerically by using the Lax-Wendroff formula, corrected by a “flux limiter”, described in detail in [9]. The idea is improve the accuracy of 1st order upwinding by adding an anti-diffusive correction term. The correction term is modulated by aflux limiter in order to prevent oscillations.
More formally, the scheme uses the following update equation
Here Qi approximates the average density of the conserved quantity q in the cell [(i−½)Δx,(i+½)Δx]. The numerical flux function Hi+1/2 approximates the flux H(q) across the boundary between two neighboring cells (i,i+1). The numerical flux function defined by the Lax-Wendroff formula is given by
HLW(Qn)=HU(Qn)+C(Qn),
where HU is the 1st order upwinding flux function defined by
Here si+1/2 approximates the flux speed ∂H/∂q at Δx(i+½), and C(Qn) is an anti-diffusive correction term given, in its limited form, by
where mml is a flux limiter (similar to the commonly used minmod limiter [9]) given by:
where:
As this scheme is based on the Lax-Wendroff formula, it is subject to a time step restriction: Δt<Δx/s, where s is the maximal flux speed.
We use this scheme to approximate the advection terms in the smoke and momentum equations, as well as the smoke gathering term. In the case of advection, the flux speed is simply the velocity along an axis. Therefore, when advecting a quantity along a particular axis, we need the proper velocity component, approximated on the faces perpendicular to that axis. In the case of smoke advection, the cells around the smoke variables coincide with the computational grid cells, so the flux speed is simply given by the velocity variables. When advecting the momentum components U,V,W, we obtain the flux speed on the boundaries between the cells around these variables by the proper averaging.
In the gathering case the flux speed should be evaluated on the faces of computational grid cells:
| Number | Name | Date | Kind |
|---|---|---|---|
| 5999194 | Brunelle | Dec 1999 | A |
| 6266071 | Stam et al. | Jul 2001 | B1 |
| 6959269 | Welterlen | Oct 2005 | B1 |
| 7184051 | Matsumoto et al. | Feb 2007 | B1 |
| 7243057 | Houston et al. | Jul 2007 | B2 |
| Number | Date | Country | |
|---|---|---|---|
| 20050253854 A1 | Nov 2005 | US |
| Number | Date | Country | |
|---|---|---|---|
| 60571619 | May 2004 | US |
