CN101484922A - Method of simulating detailed movements of fluids using derivative particles - Google Patents

Method of simulating detailed movements of fluids using derivative particles Download PDF

Info

Publication number
CN101484922A
CN101484922A CNA2006800548768A CN200680054876A CN101484922A CN 101484922 A CN101484922 A CN 101484922A CN A2006800548768 A CNA2006800548768 A CN A2006800548768A CN 200680054876 A CN200680054876 A CN 200680054876A CN 101484922 A CN101484922 A CN 101484922A
Authority
CN
China
Prior art keywords
mrow
msub
math
particle
advection
Prior art date
Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
Pending
Application number
CNA2006800548768A
Other languages
Chinese (zh)
Inventor
宋五泳
金度烨
高亨锡
Current Assignee (The listed assignees may be inaccurate. Google has not performed a legal analysis and makes no representation or warranty as to the accuracy of the list.)
Seoul National University Hospital
Original Assignee
Seoul National University Hospital
Priority date (The priority date is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the date listed.)
Filing date
Publication date
Application filed by Seoul National University Hospital filed Critical Seoul National University Hospital
Publication of CN101484922A publication Critical patent/CN101484922A/en
Pending legal-status Critical Current

Links

Images

Classifications

    • GPHYSICS
    • G06COMPUTING; CALCULATING OR COUNTING
    • G06FELECTRIC DIGITAL DATA PROCESSING
    • G06F30/00Computer-aided design [CAD]
    • G06F30/20Design optimisation, verification or simulation
    • G06F30/23Design optimisation, verification or simulation using finite element methods [FEM] or finite difference methods [FDM]
    • GPHYSICS
    • G06COMPUTING; CALCULATING OR COUNTING
    • G06FELECTRIC DIGITAL DATA PROCESSING
    • G06F2111/00Details relating to CAD techniques
    • G06F2111/10Numerical modelling

Landscapes

  • Engineering & Computer Science (AREA)
  • Physics & Mathematics (AREA)
  • Theoretical Computer Science (AREA)
  • Computer Hardware Design (AREA)
  • Evolutionary Computation (AREA)
  • Geometry (AREA)
  • General Engineering & Computer Science (AREA)
  • General Physics & Mathematics (AREA)
  • Management, Administration, Business Operations System, And Electronic Commerce (AREA)

Abstract

A method of using derivative particles to simulate the detailed movement of fluid. The present invention provides a method, which provides a new fluid simulation technique that significantly reduces the non-physical dissipation of velocity using particles and derivative information. In solving the conventional Navier-Stokes equations, the method replacing the advection part with a particle simulation. When swapping between the grid-based and particle-based simulators, the physical quantities such as the level set and velocity must be converted. A novel dissipation-suppressing conversion procedure that utilizes the derivative information stored in the particles as well as in the grid points is developed. Through several experiments, the proposed technique can reproduce the detailed movements of high-Reynolds-number fluids, such as droplets/bubbles, thin water sheets, and whirlpools. The increased accuracy in the advection, which forms the basis of the proposed technique, can also be used to produce better results in larger scale fluid simulations.

Description

Method for simulating detailed fluid motion using derivative particles
Technical Field
The present invention relates to a method of simulating detailed movement of a fluid using derivative particles.
Background
1 introduction
When water interacts violently with solids, air, or the water itself, it assumes a variety of configurations, including droplets/bubbles, thin films of water, and eddies, as shown in fig. 7. Such a characteristic can occur when the fluid undergoes motion characterized by a high Reynolds number, where the Reynolds number represents the relative magnitude of the ratio of inertia to the viscosity of the fluid. Water moving at high speed is a typical high reynolds number fluid. The present invention explores a physics-based simulation of such phenomena.
Assuming that the Navier-Stokes equations correctly model the physical motion of a fluid, a plausible simulation should be able to reproduce the high Reynolds number behavior of water. Although this work primarily discusses water, it is applicable to any fluid. However, these phenomena have not been satisfactorily reproduced so far. The present invention recognizes the factors that cause this failure and presents a method that allows for a realistic simulation of high Reynolds number liquid motion.
Simulation of seemingly unrealistic high Reynolds number fluids is related to numerical dissipation. In particular, discretized modeling of the navier-stokes equations inevitably requires estimation of physical quantities at non-grid points. In most methods, such values are calculated by interpolation of values of physical quantities at grid points. However, the error introduced by this approximation acts like a non-physical viscosity or numerical dissipation. This dissipation can be reduced by using a thinner grid; however, increasing the grid resolution increases computation time and memory requirements to impractical levels. Over the past few years, a number of excellent techniques have been proposed to address this problem. Introducing particles into the euler scheme can help capture inertial motion of the fluid and increase the effective resolution of the simulation. Enright et al [2002b ] introduced a particle level set method that increases the accuracy of surface tracking by spreading particles around the interface. Losasso et al [2004] propose an octree-based multi-level fluid solver that allows finer resolution simulations in more interesting areas such as water surfaces.
Unfortunately, the above techniques do not produce high Reynolds number liquids with sufficient detail and realism. The simulated fluid exhibits a greater viscosity than in real physics; fluids often look thick/viscous and do not typically exhibit movement of small scale features in complex flows. The applicant has found that the reason why the above dissipation suppression techniques produce such artifacts is that they are not able to effectively suppress velocity dissipation, even if they significantly reduce mass dissipation.
In the present invention, the applicant introduced a new concept called derivative particle (derivative particle) and developed a fluid simulator based on the concept. The applicant utilized the non-dissipative nature of the simulated raleigh's scheme for the advection part; for fluid regions that require simulation details, applicants use particles to solve the advection step. Switching between grid-based and particle-based simulators may introduce numerical dissipation. The present invention develops a new conversion procedure that allows the reproduction of detailed fluid movements. An innovative aspect of this work is that in addition to storing the physical quantities (velocity and level set values), the derivative particles also store the spatial derivatives of those quantities, which enables a more accurate assessment of the physical quantities at non-grid/non-particle locations. The use of particles for this work differs from the particle level set approach [ Enright et al 2002b ] in that the derivative particles carry the fluid velocity as well as the level set value.
Experiments show that the proposed simulator accurately tracks the interface. More interestingly, the proposed method proved to significantly reduce non-physical damping, allowing the reproduction of macroscopic motions of small scale features that occur in real high reynolds number fluids.
The remainder of the description is organized as follows: section 2 reviews previous work; section 3 gives an overview of the simulator; section 4 introduces an octree-based Constrained Interpolation Profile (CIP) solver; section 5 introduces derivative particle models; section 6 reports our experimental results; section 7 summarizes the description.
2 previous work
Since Foster and Metaxas [ 1996; 1997] first introduced a fluid animation technique based on a full 3D navier-stokes simulation, so this approach has prevailed in the computer graphics world. Jos Stam [1999] introduced a Stable fluid solver known as a Stable fluid. The advection step of the solver is implemented using the semi-Lagrangian method [ Stanifforth and C ^ ot' e 1991], which remains stable even when large time steps are used. Since then, there has been active development of fast fluid animation techniques based on the semi-lagrange method in computer graphics. Rasmussen et al [2003] propose a technique for generating 3D large-scale animations of gases using a 2D semi-lagrange solver, and Feldman et al [2003] propose an explosion model incorporating a particle-based combustion model into a stable fluid solver. Treuille et al [ 2003; 2004, to et al, introduces an optimization technique that generates fluid flows that satisfy specified key frame constraints. Additionally, Feldman et al propose the use of a hybrid mesh combining interlaced meshes and unstructured meshes [ Feldman et al 2005a ], but also the use of deformable meshes [ Feldman et al 2005b ], for which they must modify the semi-Lagrangian method.
To simulate a fluid, in addition to a stable Navier-Stokes solver, a model for tracking the surface of the liquid is required. For this purpose, Foster and Fedkiw [2001] propose a mixed surface model that introduces mass-free particles into the level set field. This model motivated the development of a particle level set method [ Enright et al 2002b ] that can capture the dynamic motion of a fluid surface with significant accuracy. Enright et al [2005] demonstrated that this particle level set approach allows the use of a large time step in the semi-lagrange scheme. This method has been used to model the interaction between fluids and rigid bodies [ Carlson et al 2004] and between fluids and deformable shell bodies [ guillelman et al 2005] and to simulate viscoelastic fluids [ Goktekin et al 2004], sand [ Zhu and Bridson 2005], multiphase fluids [ Hong and Kim 2005], and water droplets [ Wang et al 2005 ].
Instead of grid-based methods, purely particle-based methods have also been investigated. Stam and Fiume [1995] used Smooth Particle Hydrodynamics (SPH) to model the gaseous phenomenon. In SPH, the fluid is represented by a collection of particles, and the simulator calculates their motion by calculating each term of the navier-stokes equation. M. Uller et al [2003] use SPH models to simulate fluids, and in [ M1 Uller et al 2005] they use this technique to simulate multiphase fluid interactions. Premo ˇ ze et al [2003] introduced the moving-particle semi-implicit (MPS) method to the graphical community, a technique that showed better performance than SPH in simulating incompressible motion of a fluid. However, a fundamental drawback of purely particle-based methods is their lack of surface tracking capability; if an inappropriate number of particles are used, they may produce a grainy or excessively smooth surface.
The main factor that impairs the visual realism (and physical accuracy) of the simulation results is the numerical dissipation of the velocity field. To reduce dissipation when simulating the gaseous phenomenon, Fedkiw et al [2001] used cubic interpolation. They also include an additional step called vorticity constraint (vorticity constraint) that amplifies the vorticity (curl) of the velocity field, producing a true vortex-shaped portion in the flue gas motion. A more physical-based prevention of vorticity dissipation by employing the vortex (vortex) particle method [ Selle et al 2005; park and Kim 2005 ]. This method works with the vorticity-type navier-stokes equation and solves for the (solve) advection term with the particles, which results in an efficient preservation of vorticity. However, in the case of liquids, the viscosity shows a more pronounced effect; in particular, the swirling motion is short lived and less observed. Therefore, modeling the behavior of high Reynolds number liquids requires velocity dissipation suppression techniques that work in more general cases. To reduce dissipation in liquid simulation, Song et al [2005] proposed a technique based on the CIP advection method. Their technique solves the velocity advection with third order accuracy. However, with this method, the numerical viscosity from the mesh-based advection cannot be avoided. Zhu and Bridson [2005] proposed a FLuid simulation technique based on the FLuid-Implicit-Particle (FLIP) method. The FLIP method [ Brackbill and Ruppel 1986] uses particles to solve the advection step, and does not include all other steps of the grid. Although the particle advection in this approach has no numerical dissipation, interpolation errors are introduced when passing velocity values back and forth between the grid and the particles.
Therefore, there is a long felt need for a method of simulating detailed motion of a fluid using derivative particles. The present invention is intended to solve these problems and meet the long-felt need.
Disclosure of Invention
The present invention seeks to address the disadvantages of the prior art.
It is an object of the present invention to provide a method for simulating detailed movements of a fluid using derivative particles.
It is another object of the present invention to provide a method for simulating detailed motion of a fluid using derivative particles, wherein the advection part is implemented with the lagrangian scheme.
It is a further object of the present invention to provide a method for simulating detailed movements of fluids using derivative particles, which provides a procedure for converting physical quantities at the switching of advection and non-advection parts, suppressing the intervention of unnecessary numerical dissipation.
Overview of the invention 3
In the development of multiphase fluid solvers, the applicant assumed that both air and water were incompressible. The Navier-Stokes equation for incompressible fluids can be written as
<math> <mrow> <mfrac> <mrow> <mo>&PartialD;</mo> <mi>u</mi> </mrow> <mrow> <mo>&PartialD;</mo> <mi>t</mi> </mrow> </mfrac> <mo>=</mo> <mo>-</mo> <mrow> <mo>(</mo> <mi>u</mi> <mo>&CenterDot;</mo> <mo>&dtri;</mo> <mo>)</mo> </mrow> <mi>u</mi> <mo>-</mo> <mfrac> <mrow> <mo>&dtri;</mo> <mi>p</mi> </mrow> <mi>&rho;</mi> </mfrac> <mo>+</mo> <mo>&dtri;</mo> <mo>&CenterDot;</mo> <mrow> <mo>(</mo> <mi>&mu;</mi> <mo>&dtri;</mo> <mi>u</mi> <mo>)</mo> </mrow> <mo>+</mo> <mfrac> <mi>f</mi> <mi>&rho;</mi> </mfrac> <mo>-</mo> <mo>-</mo> <mo>-</mo> <mrow> <mo>(</mo> <mn>1</mn> <mo>)</mo> </mrow> </mrow></math>
<math> <mrow> <mo>&dtri;</mo> <mo>&CenterDot;</mo> <mi>u</mi> <mo>=</mo> <mn>0</mn> <mo>-</mo> <mo>-</mo> <mo>-</mo> <mrow> <mo>(</mo> <mn>2</mn> <mo>)</mo> </mrow> </mrow></math>
Where u is velocity, p is pressure, f is external force, μ is viscosity coefficient, and ρ is fluid density. To accurately model density and viscosity discontinuities across an interface between media, applicants employed a virtual fluid method with variable density [ Liu et al 2000; kang et al 2000; hong and Kim 2005 ]. In our solver, surface tracking is based on a level set method. Updating the level set field phi according to
<math> <mrow> <mfrac> <mrow> <mo>&PartialD;</mo> <mi>&phi;</mi> </mrow> <mrow> <mo>&PartialD;</mo> <mi>t</mi> </mrow> </mfrac> <mo>+</mo> <mi>u</mi> <mo>&CenterDot;</mo> <mo>&dtri;</mo> <mi>&phi;</mi> <mo>=</mo> <mn>0</mn> <mo>-</mo> <mo>-</mo> <mo>-</mo> <mrow> <mo>(</mo> <mn>3</mn> <mo>)</mo> </mrow> </mrow></math>
Symbol distance condition
<math> <mrow> <mo>|</mo> <mo>&dtri;</mo> <mi>&phi;</mi> <mo>|</mo> <mo>=</mo> <mn>1</mn> <mo>-</mo> <mo>-</mo> <mo>-</mo> <mrow> <mo>(</mo> <mn>4</mn> <mo>)</mo> </mrow> </mrow></math>
For the time integration of the Naviger-Stokes equation, applicants employ a step-by-step method [ Stam1999], which incrementally calculates (accounts for) the effect of the term in equation (1) and the mass conservation condition in equation (2). The techniques developed by the applicant are focused here on the increase of the advection part. Thus, as shown in FIG. 8, applicants have split (fractional steps) into two parts, a non-advection part and an advection part. The time integration of equation (3) involves only the advection part.
The advection part consists of a grid advection part and a particle advection part. The mesh advection part advects the velocity and level set fields computed from the non-advection part according to the current velocity field. Our method differs from pure euler simulation in that our simulator includes a particle advection part that is used to simulate the area where detail needs to be generated. To this end, applicants have introduced a large number of particles in the air-water interface.
In the mesh advection part, advection transport speed and advection transport level set are euler advection steps. The precise processing of these steps forms the basis for successful simulation of high Reynolds number behaviour. For euler advection, the applicant developed an octree-based CIP solver (part 4) that reduces the speed and mass dissipation to a significant level.
When the advection of a fluid region needs to be simulated using the particle advection part, the physical quantities (e.g., level set and velocity) currently defined on the grid are delivered to the particles. After this transport, the particles can be transported in parallel flow directly. The results of the particle advection are then transmitted back to the grid. At this point in the program, the final velocity and level set values are stored at the grid points, regardless of whether the simulation was conducted via the grid advection part or the particle advection part. Applicants describe the grid-to-particle and particle-to-grid speed/level set conversion procedures in section 5. The use of particle-based advection and the development of conversion programs are essential to reproduce the details of fluid motion.
A method of simulating detailed movement of a fluid using derivative particles comprising the steps of: a) modeling the fluid with an adaptive mesh based on an octree data structure; b) solving the Navier-Stokes equations for the incompressible fluid for the fluid velocities at the grid points; c) using a step-wise method in the time integration of the Navier-Stokes equation; d) dividing the steps into a non-advection part and an advection part; e) selecting a grid advection part or a particle advection part according to the detail degree of the simulation; f) using octree-based CIP method for the mesh advection part; g) derivative particle models are used for the particle advection part.
When the level of detail of the simulation is high, the step of using the derivative particle model for the particle advection part is used.
Adapting the mesh to the details of the region of high interest increases the number of meshes.
The octree data structure comprises a tree data structure with up to eight children per internal node.
The Navier-Stokes equation for incompressible fluids can be written as
<math> <mrow> <mfrac> <mrow> <mo>&PartialD;</mo> <mi>u</mi> </mrow> <mrow> <mo>&PartialD;</mo> <mi>t</mi> </mrow> </mfrac> <mo>=</mo> <mo>-</mo> <mrow> <mo>(</mo> <mi>u</mi> <mo>&CenterDot;</mo> <mo>&dtri;</mo> <mo>)</mo> </mrow> <mi>u</mi> <mo>-</mo> <mfrac> <mrow> <mo>&dtri;</mo> <mi>p</mi> </mrow> <mi>&rho;</mi> </mfrac> <mo>+</mo> <mo>&dtri;</mo> <mo>&CenterDot;</mo> <mrow> <mo>(</mo> <mi>&mu;</mi> <mo>&dtri;</mo> <mi>u</mi> <mo>)</mo> </mrow> <mo>+</mo> <mfrac> <mi>f</mi> <mi>&rho;</mi> </mfrac> <mo>,</mo> </mrow></math>
<math> <mrow> <mo>&dtri;</mo> <mo>&CenterDot;</mo> <mi>u</mi> <mo>=</mo> <mn>0</mn> <mo>,</mo> </mrow></math>
Where u is velocity, p is pressure, f is external force, μ is viscosity coefficient, and ρ is fluid density.
The method may further comprise the step of tracking the surface of the fluid using a level set method. The level set method includes a level set field phi. Updating the level set field φ according to:
<math> <mrow> <mfrac> <mrow> <mo>&PartialD;</mo> <mi>&phi;</mi> </mrow> <mrow> <mo>&PartialD;</mo> <mi>t</mi> </mrow> </mfrac> <mo>+</mo> <mi>u</mi> <mo>&CenterDot;</mo> <mo>&dtri;</mo> <mi>&phi;</mi> <mo>=</mo> <mn>0</mn> </mrow></math>
distance of symbol <math> <mrow> <mo>|</mo> <mo>&dtri;</mo> <mi>&phi;</mi> <mo>|</mo> <mo>=</mo> <mn>1</mn> <mo>.</mo> </mrow></math>
<math> <mrow> <mfrac> <mrow> <mo>&PartialD;</mo> <mi>&phi;</mi> </mrow> <mrow> <mo>&PartialD;</mo> <mi>t</mi> </mrow> </mfrac> <mo>+</mo> <mi>u</mi> <mo>&CenterDot;</mo> <mo>&dtri;</mo> <mi>&phi;</mi> <mo>=</mo> <mn>0</mn> </mrow></math> The time integral of (a) involves only the advection part. The mesh advection part advects the velocity and level set fields computed from the non-advection part according to the current velocity field.
The grid advection part advects the delivery speed and level set using an Euler advection step. Euler advection is performed by the octree-based CIP (constrained interpolation profile) method.
The step of using the octree-based CIP method for the grid advection part comprises the step of advecting the velocity and level sets and their derivatives with the octree-based CIP method. Direct from equation <math> <mrow> <mfrac> <mrow> <mo>&PartialD;</mo> <mi>&phi;</mi> </mrow> <mrow> <mo>&PartialD;</mo> <mi>t</mi> </mrow> </mfrac> <mo>+</mo> <mi>u</mi> <mo>&CenterDot;</mo> <mo>&dtri;</mo> <mi>&phi;</mi> <mo>=</mo> <mn>0</mn> </mrow></math> Obtaining an equation for the advection transport derivative, where φ is the physical quantity advected, and obtaining an advection transport equation for the derivative of the spatial variable x by differentiating the equation for x, i.e.:
<math> <mrow> <mfrac> <mrow> <mo>&PartialD;</mo> <msub> <mi>&phi;</mi> <mi>x</mi> </msub> </mrow> <mrow> <mo>&PartialD;</mo> <mi>t</mi> </mrow> </mfrac> <mo>+</mo> <mi>u</mi> <mo>&CenterDot;</mo> <mo>&dtri;</mo> <msub> <mi>&phi;</mi> <mi>x</mi> </msub> <mo>=</mo> <mo>-</mo> <msub> <mi>u</mi> <mi>x</mi> </msub> <mo>&CenterDot;</mo> <mo>&dtri;</mo> <mi>&phi;</mi> <mo>.</mo> </mrow></math>
the step of solving the advection equation for the derivative of x comprises the step of using a step-by-step method, which comprises the steps of: a) solving for non-advective parts using finite difference method <math> <mrow> <mo>&PartialD;</mo> <msub> <mi>&phi;</mi> <mi>x</mi> </msub> <mo>/</mo> <mo>&PartialD;</mo> <mi>t</mi> <mo>=</mo> <mo>-</mo> <msub> <mi>u</mi> <mi>x</mi> </msub> <mo>&CenterDot;</mo> <mo>&dtri;</mo> <mi>&phi;</mi> <mo>;</mo> </mrow></math> b) Advective conveying device <math> <mrow> <mfrac> <mrow> <mo>&PartialD;</mo> <msub> <mi>&phi;</mi> <mi>x</mi> </msub> </mrow> <mrow> <mo>&PartialD;</mo> <mi>t</mi> </mrow> </mfrac> <mo>+</mo> <mi>u</mi> <mo>&CenterDot;</mo> <mo>&dtri;</mo> <msub> <mi>&phi;</mi> <mi>x</mi> </msub> <mo>=</mo> <mn>0</mn> </mrow></math> The result of (1); and c) sampling pressure values at the center of cells (cells) defined by the grid and sampling other values including velocity, level sets and their derivatives at nodes.
The method may further comprise the steps of: a) combining an octree-based CIP method with an adaptive grid; and b) minimizing octree artifacts caused by regional variations in grid resolution of the octree data structure by simulating advection with the octree-based CIP method.
The derivative particle model uses a particle-based lagrangian scheme. By using <math> <mrow> <mi>&phi;</mi> <mrow> <mo>(</mo> <mi>x</mi> <mo>)</mo> </mrow> <mo>=</mo> <msub> <mi>&phi;</mi> <mi>p</mi> </msub> <mo>+</mo> <mfrac> <mrow> <mo>&dtri;</mo> <msub> <mi>&phi;</mi> <mi>p</mi> </msub> </mrow> <mrow> <mo>|</mo> <mo>&dtri;</mo> <msub> <mi>&phi;</mi> <mi>p</mi> </msub> <mo>|</mo> </mrow> </mfrac> <mo>&CenterDot;</mo> <mrow> <mo>(</mo> <mi>x</mi> <mo>-</mo> <msub> <mi>x</mi> <mi>p</mi> </msub> <mo>)</mo> </mrow> </mrow></math> And <math> <mrow> <mo>&dtri;</mo> <mi>&phi;</mi> <mrow> <mo>(</mo> <mi>x</mi> <mo>)</mo> </mrow> <mo>=</mo> <mo>&dtri;</mo> <msub> <mi>&phi;</mi> <mi>p</mi> </msub> </mrow></math> the level set and its derivative at the grid points are computed, where p is the gradient carried by the derivative particles.
The method may further comprise the steps of: a) converting the defined velocity and level sets on the grid to defined velocity and level sets on the particles (grid to particles) prior to the step of using the derivative particle model for the particle advection part; and b) converting the defined velocity and level sets on the particles to defined velocity and level sets on the grid (particle to grid) before returning to the step of using the octree-based CIP model for the grid advection part.
The step of converting the velocity and level set grid to particles comprises the step of, for velocity uPUsing the steps of the FLIP method with monotonic CIP interpolation replacing linear interpolation, and <math> <mrow> <msup> <mi>u</mi> <mi>P</mi> </msup> <mo>=</mo> <msubsup> <mi>u</mi> <mo>-</mo> <mi>P</mi> </msubsup> <mo>+</mo> <msubsup> <mi>&Sigma;</mi> <mrow> <mi>i</mi> <mo>=</mo> <mn>1</mn> </mrow> <mn>4</mn> </msubsup> <msub> <mi>w</mi> <mi>i</mi> </msub> <mi>&Delta;</mi> <msubsup> <mi>u</mi> <mi>i</mi> <mi>G</mi> </msubsup> <mo>,</mo> </mrow></math> wherein
Figure A200680054876D00142
Is the particle velocity just before the start of the non-advection part, and
Figure A200680054876D00143
is due to a non-advective part of GiThe speed of (c) is changed.
For speed uPGrid points G, nearest particle P selected from each quadrant of G1、P2、P3And P4、PiVelocity of and
Figure A200680054876D00144
and
Figure A200680054876D00145
the step of converting the velocity and level set particles to a grid comprises the steps of: a) will be provided with
Figure A200680054876D00146
Derivative of (2)Coordinate transformation to
Figure A200680054876D00148
Thereby to obtainIndicates along the direction
Figure A200680054876D00149
Component (a) ofIs shown perpendicular to
Figure A200680054876D001411
Wherein similarly obtained
Figure A200680054876D001412
Derivative of the coordinate transformation of
Figure A200680054876D001413
b) Along the directionPerforming one-dimensional monotonic CIP interpolation on the result obtained in step (a) to calculate u at position AAAnd it
Figure A200680054876D001415
A directional derivative M; c) by applying the inverse (inverse) of the above-mentioned coordinate transformation
Figure A200680054876D001416
Obtaining (u)Ax,uAy) (ii) a d) To particle P3And P4Performing steps (a) - (c) to calculate u of BBAnd (u)Bx,uBy) (ii) a And e) performing a y-direction monotonic CIP interpolation on the results obtained in steps (c) and (d) which gives the velocity and derivative at G.
The step of using a derivative particle model for the particle advection part comprises the steps of: a) adjusting a level set value of the particle; and b) performing particle reseeding.
The advantages of the present invention are (1) the method of simulating detailed motion of a fluid using derivative particles reduces dissipation of both mass and velocity, resulting in a true reproduction of the dynamic fluid; (2) the method of simulating detailed movement of a fluid using derivative particles implements the advection part with the lagrangian scheme; and (3) the method of simulating detailed motion of a fluid using derivative particles provides a particle-to-grid and grid-to-particle procedure.
While the present invention has been briefly summarized, a more complete understanding of the invention can be obtained from the following drawings, detailed description, and claims.
Drawings
These and other features, aspects, and advantages of the present invention will become better understood with reference to the accompanying drawings.
FIG. 1 is a flow chart illustrating a method of simulating detailed movement of a fluid using derivative particles in accordance with the present invention;
FIG. 2 is another flow chart illustrating a method according to the present invention;
FIG. 3 is a further flow chart illustrating a method according to the present invention;
FIG. 4 is a partial flow diagram showing grid-to-particle and particle-to-grid conversion;
FIG. 5 is another partial flow diagram showing grid-to-particle and particle-to-grid conversion;
FIG. 6 is a partial flow diagram showing details of a step of using a derivative particle model;
FIG. 7 shows simulation results of water produced by the impact of a ball;
FIG. 8 illustrates the architecture of a method according to the invention;
FIG. 9 shows the calculation from particle velocity to grid velocity;
FIG. 10 shows a screenshot taken from a 2D dam break (breaking-dam) simulation: the upper and lower sequences show the results produced with a conventional linear model and derivative particle model of the particle level set method, respectively;
FIG. 11 shows a screenshot taken from a water drop simulation: the left and right images show the results produced with a conventional linear model and derivative particle model of the particle level set method, respectively;
FIG. 12 shows a screenshot taken from a simulation of a rotating water tank; and
FIG. 13 shows a screenshot taken from a flood simulation: the three images of the bottom row show the progression of the flood (top view); the top image shows a side view.
Detailed Description
4 development octree-based CIP solver
This section describes the development of the mesh advection part of the simulator by combining the CIP method with the octree data structure.
4.1 CIP method introduction
Advection terms in solving the Navell-Stokes equations (part 3) in a step-wise manner
Figure A200680054876D00151
And
Figure A200680054876D00152
special attention is required due to its hyperbolic nature. The semi-Lagrangian method provides a stable solution in which to handle the above hyperbolic equation [ Stam1999]]. Unfortunately, however, this approach suffers from severe numerical dissipation resulting from linear interpolation of the physical quantities at non-grid points used to determine a back-pull (backtracking) from the physical quantities at adjacent grid points.
While the CIP process was originally developed by Yabe and Aoki [1991 a; 1991b ] A third order technique was proposed and subsequently improved by Yabe et al [2001 ]. The key idea of this method is to advect not only its physical quantity but also its derivative. Here, the problem would be how to advect the derivative. An interesting observation of Yabe and Aoki [1991a ] is that the equation for advective transport derivatives can be obtained directly from the initial hyperbolic equation
<math> <mrow> <mfrac> <mrow> <mo>&PartialD;</mo> <mi>&phi;</mi> </mrow> <mrow> <mo>&PartialD;</mo> <mi>t</mi> </mrow> </mfrac> <mo>+</mo> <mi>u</mi> <mo>&CenterDot;</mo> <mo>&dtri;</mo> <mi>&phi;</mi> <mo>=</mo> <mn>0</mn> <mo>-</mo> <mo>-</mo> <mo>-</mo> <mrow> <mo>(</mo> <mn>5</mn> <mo>)</mo> </mrow> </mrow></math>
Where φ is the physical quantity being advected. For the spatial variable x differential equation (5), applicants have obtained
<math> <mrow> <mfrac> <mrow> <mo>&PartialD;</mo> <msub> <mi>&phi;</mi> <mi>x</mi> </msub> </mrow> <mrow> <mo>&PartialD;</mo> <mi>t</mi> </mrow> </mfrac> <mo>+</mo> <mi>u</mi> <mo>&CenterDot;</mo> <mo>&dtri;</mo> <msub> <mi>&phi;</mi> <mi>x</mi> </msub> <mo>=</mo> <mo>-</mo> <msub> <mi>u</mi> <mi>x</mi> </msub> <mo>&CenterDot;</mo> <mo>&dtri;</mo> <mi>&phi;</mi> <mo>-</mo> <mo>-</mo> <mo>-</mo> <mrow> <mo>(</mo> <mn>6</mn> <mo>)</mo> </mrow> </mrow></math>
Which can be used to predict the evolution of x over time. In solving equation (6), applicants again apply the step-wise approach: first, applicants use finite differences to solve for non-advection parts <math> <mrow> <mo>&PartialD;</mo> <msub> <mi>&phi;</mi> <mi>x</mi> </msub> <mo>/</mo> <mo>&PartialD;</mo> <mi>t</mi> <mo>=</mo> <mo>-</mo> <msub> <mi>u</mi> <mi>x</mi> </msub> <mo>&CenterDot;</mo> <mo>&dtri;</mo> <mi>&phi;</mi> <mo>;</mo> </mrow></math> Applicant then advected the results according to the following equation
<math> <mrow> <mfrac> <mrow> <mo>&PartialD;</mo> <msub> <mi>&phi;</mi> <mi>x</mi> </msub> </mrow> <mrow> <mo>&PartialD;</mo> <mi>t</mi> </mrow> </mfrac> <mo>+</mo> <mi>u</mi> <mo>&CenterDot;</mo> <mo>&dtri;</mo> <msub> <mi>&phi;</mi> <mi>x</mi> </msub> <mo>=</mo> <mn>0</mn> <mo>-</mo> <mo>-</mo> <mo>-</mo> <mrow> <mo>(</mo> <mn>7</mn> <mo>)</mo> </mrow> </mrow></math>
A detailed description of CIP advection can be found in Song et al 2005. An extension of the two-dimensional and three-dimensional aspects of the above methods exists [ Yabe and Aoki 1991b ]. Song et al [2005] found that the generalization given in [1991b ] may cause instability and proposed a modification to address this problem. In this work, to solve the advection part, applicants employed the CIP method with a second order Runge-Kutta time integral.
4.2 combining CIP with adaptive octree Trees
The use of an adaptive mesh based on an octree data structure introduced by Losasso et al [2004] allows simulations to be performed with non-homogeneous accuracy in all fluids. From a practical point of view, this approach is useful because it allows simulating more interesting areas, such as surfaces, with higher accuracy by introducing a small amount of additional computation and memory. Here, applicants take this octree data structure and propose that the actual values of this adaptive approach can be further improved by combining it with the CIP approach.
Losasso et al [2004] use linear interpolation for the half lagrange advection step. However, if applicants replace linear interpolation with CIP interpolation, then the advection will have third order accuracy. As in Guendelman et al [2005], applicants sampled pressure values at the cell center, but sampled all other quantities-velocity, level set, and derivatives thereof-at the junction. When the cell is refined, CIP interpolation is performed on the values at the new grid points by referring to all the derivative values.
We note that the CIP method fits very well with the octree data structure because, although it has third-order spatial precision, it is confined to a single grid cell template (Tencel) rather than multiple templates. Due to this compactness, applicants are able to use the CIP method for simulating an adaptive mesh without any major modifications. In contrast, for higher order schemes defined on multiple templates, it is difficult to extend the method from first order to third order: since the simulation adaptively refines the mesh, cells with different sizes will be generated, which makes the development of multi-template high-order schemes a difficult problem.
We will also note a feature that applicant calls for octree artifacts. Regional variations in the grid resolution of the octree data structure produce variations in the amount of dissipation. The regional variation in mass dissipation is not significant. However, for velocity dissipation, regional differences can be significant, especially for fast fluid motion. For example, in the simulation of a dam break, water undergoes rapid lateral movement. Fluids near high resolution surfaces undergo a small amount of numerical diffusion and thus undergo rapid movement. Instead, the fluid at the bottom moves in a pattern characteristic of a more viscous fluid due to the low mesh resolution. When these two types of motion occur together, the upper portion of the water appears to crawl over the lower portion. High-order schemes, including CIP, do not avoid such artifacts. Then, when advection is simulated using the octree-based CIP solver proposed in this work, the artifacts are less noticeable; applicants attribute this improvement to an overall reduction in numerical dissipation, especially in the low resolution region.
5 derivative particle model
The particle-based Lagrangian scheme is more accurate than the mesh-based Euler scheme in terms of momentum/mass conservation. Therefore, the applicant applies the lagrangian scheme to regions that need to be simulated in more detail. However, this method has its own limitations in surface tracking and pressure/viscosity calculations. Therefore, the applicant applies this method only to the simulation of the advection part; all other parts of the simulation are based on the euler scheme. Switching between the lagrange and euler schemes requires grid-to-particle and particle-to-grid conversion of physical quantities (velocity and level set values). The procedures for performing these conversions should be carefully developed so that they do not introduce unnecessary numerical diffusion. In this section, applicants introduced a novel conversion procedure that does not impair the desirable properties of Lagrangian advection, which is a necessary component to enable the reproduction of small scale features in high Reynolds number fluids. For simplicity, the description in this section is made for 2D.
5.1 mesh to particle velocity conversion
At the end of the advection part shown in fig. 8, (1) the velocities that need to be advected are stored at the grid points, and (2) a large number of particles are scattered over the grid. The particle advection part must find the velocity of the particles so that their velocity and level set values produce lagrangian motion. Suppose particle P is at four grid points G1、G2、G3And G4In defined units. For i ∈ {1, 2, 3, 4}, let uG i=(uG i,vG i) Is GiThe grid speed of (d). The applicant is looking for a speed u that can be used for PP=(uP,vP) The formula (2).
In-cell mass point method [ Harlow 1963]One possible method used in (1) is to use bilinear interpolation <math> <mrow> <msup> <mi>u</mi> <mi>P</mi> </msup> <mo>=</mo> <msubsup> <mi>&Sigma;</mi> <mrow> <mi>i</mi> <mo>=</mo> <mn>1</mn> </mrow> <mn>4</mn> </msubsup> <msub> <mi>w</mi> <mi>i</mi> </msub> <msubsup> <mi>u</mi> <mi>i</mi> <mi>G</mi> </msubsup> <mo>,</mo> </mrow></math> Wherein, wiIs a bilinear weighting determined based on the particle position P. Unfortunately, this conversion introduces severe numerical diffusion. In the FLIP method [ Brackbill and Ruppel 1986]Only the increasing portion of the grid velocity is transmitted to the site velocity. That is to say that the first and second electrodes, <math> <mrow> <msup> <mi>u</mi> <mi>P</mi> </msup> <mo>=</mo> <msubsup> <mrow> <msubsup> <mi>u</mi> <mo>-</mo> <mi>P</mi> </msubsup> <mo>+</mo> <mi>&Sigma;</mi> </mrow> <mrow> <mi>i</mi> <mo>=</mo> <mn>1</mn> </mrow> <mn>4</mn> </msubsup> <msub> <mi>w</mi> <mi>i</mi> </msub> <msubsup> <mi>&Delta;u</mi> <mi>i</mi> <mi>G</mi> </msubsup> <mo>,</mo> </mrow></math> wherein,is the particle velocity just before the start of the non-advection part, and
Figure A200680054876D00184
is due to a non-advective part of GiThe speed of (c) is changed. This approach has been shown to significantly reduce numerical diffusion [ Zhu and Bridson 2005]. Thus, in the development of the grid-to-particle-velocity conversion program, applicants used the FLIP method, but interpolated with monotonic CIP [ Song et al 2005]Instead of linear interpolation. This gives not only uPAnd give
<math> <mrow> <mrow> <mo>(</mo> <mfrac> <mo>&PartialD;</mo> <mrow> <mo>&PartialD;</mo> <mi>x</mi> </mrow> </mfrac> <msup> <mi>u</mi> <mi>P</mi> </msup> <mo>,</mo> <mfrac> <mo>&PartialD;</mo> <mrow> <mo>&PartialD;</mo> <mi>y</mi> </mrow> </mfrac> <msup> <mi>u</mi> <mi>P</mi> </msup> <mo>)</mo> </mrow> <mi>and</mi> <mo>(</mo> <mfrac> <mo>&PartialD;</mo> <mrow> <mo>&PartialD;</mo> <mi>x</mi> </mrow> </mfrac> <msup> <mi>v</mi> <mi>P</mi> </msup> <mo>,</mo> <mfrac> <mo>&PartialD;</mo> <mrow> <mo>&PartialD;</mo> <mi>y</mi> </mrow> </mfrac> <msup> <mi>v</mi> <mi>P</mi> </msup> <mo>)</mo> </mrow></math>
5.2 particle-to-grid speed conversion
Suppose each particle is based on velocity uPMovement, the spatial derivative of the velocity carried, and the velocity itself. Applicants must now develop a program that communicates the results of particle advection to the grid. Particle-to-grid speed conversion versus grid-to-particle conversionAnd is more complex.
Referring to fig. 9, let G be a grid point. Applicant selects the closest particle, P, from each quadrant of G1、P2、P3And P4. Make itIs PiAnd is at a speed of
Figure A200680054876D00193
And
Figure A200680054876D00194
are respectively as
Figure A200680054876D00195
The derivatives of the x and y components of (a). The applicant has to find the speed u at GG=(uG,vG) And the formula of its spatial derivatives. Since the four particles are not placed rectangularly, applicants cannot use traditional grid-based CIP interpolation.
Our particle-to-grid velocity conversion consists of the following steps. For simplicity, applicants show the calculation of only the x component.
(1) The applicant will
Figure A200680054876D00196
Derivative of (2)
Figure A200680054876D00197
Coordinate transformation toAnd, thus,
Figure A200680054876D0019115629QIETU
is shown in the directionAnd u is1⊥I denotes perpendicular to
Figure A200680054876D00199
The component (c). As such, applicants have obtained
Figure A200680054876D001910
Derivative of the coordinate transformation of
Figure A200680054876D001911
(2) Applicant follows the direction
Figure A200680054876D001912
Performing one-dimensional monotonic CIP interpolation on the result obtained in step (1) to compute u at position AAAnd it
Figure A200680054876D001913
Derivative of direction
Figure A200680054876D001914
(3) Applicant derives from the inverse of the above coordinate transformation
Figure A200680054876D001915
Obtaining (u)Ax,uAy)。
(4) Similarly, Applicant pairs particle P3And P4Performing steps (1) - (3) to calculate u of BBAnd (u)Bx,uBy)。
(5) Applicants perform y-direction monotonic CIP interpolation on the results obtained in steps (3) and (4), which gives the velocity and derivative at G.
This monotonic interpolation method can be extended directly to the 3D case. Applicants have noted that this work does not employ the general radial basis interpolation scheme because (1) they do not utilize derivative information, and (2) it has proven difficult to modify the scheme so that they embody the derivatives but remain monotonous.
5.3 level set conversion
Level set conversion should be performed differently than speed conversion; the level set value represents the shortest distance to the surface, so the conversion will not be based on interpolation. Here, applicants use a widely used level set conversion program [ Enright et al 2002 a; enright et al 2002b ].
In the initial particle level set method, the spherical implicit function phi (x) is used as sp(|φp|-|x-xp| where p is the level set of particles and s is for positive (negative) particles, to compute the level set value at grid point xpAnd (1) (+ 1. In this work, using the derivative information stored in the derivative particles, applicants used the following equation to calculate the level set and its derivatives at the grid points
<math> <mrow> <mi>&phi;</mi> <mrow> <mo>(</mo> <mi>x</mi> <mo>)</mo> </mrow> <mo>=</mo> <msub> <mi>&phi;</mi> <mi>p</mi> </msub> <mo>+</mo> <mfrac> <mrow> <mo>&dtri;</mo> <msub> <mi>&phi;</mi> <mi>p</mi> </msub> </mrow> <mrow> <mo>|</mo> <mo>&dtri;</mo> <msub> <mi>&phi;</mi> <mi>p</mi> </msub> <mo>|</mo> </mrow> </mfrac> <mo>&CenterDot;</mo> <mrow> <mo>(</mo> <mi>x</mi> <mo>-</mo> <msub> <mi>x</mi> <mi>p</mi> </msub> <mo>)</mo> </mrow> <mo>-</mo> <mo>-</mo> <mo>-</mo> <mrow> <mo>(</mo> <mn>8</mn> <mo>)</mo> </mrow> </mrow></math>
<math> <mrow> <mo>&dtri;</mo> <mi>&phi;</mi> <mrow> <mo>(</mo> <mi>x</mi> <mo>)</mo> </mrow> <mo>=</mo> <mo>&dtri;</mo> <msub> <mi>&phi;</mi> <mi>p</mi> </msub> <mo>,</mo> <mo>-</mo> <mo>-</mo> <mo>-</mo> <mrow> <mo>(</mo> <mn>9</mn> <mo>)</mo> </mrow> </mrow></math>
Where p is the gradient carried by the derivative particle. Since equation (8) is based on the distance in the gradient direction rather than the euclidean distance, it produces a more accurate result than the spherical implicit function. The conversion procedure is the same as that set forth in Enright et al [2002a ], except for the above modifications.
6 test results
The technique proposed in the present invention is implemented on a Power Mac with dual G52.5 GHz processors and 5.5GB of memory. The applicant used a program to simulate a number of experimental situations in the real world that produce high reynolds number fluid behaviour. In experiments, applicants solved a two-phase fluid consisting of water and air with the following constants: g is 9.8m/sec2、ρwater=1000kg/m3、μwater=1.137×103kg/ms、ρair=1.226kg/m3And muair=1.78×105kg/ms, where g is gravity. Extraction of water surface using marching cubes algorithm, and by smart shooting
Figure A200680054876D00211
Rendering is performed.
Dam break: to compare the numerical diffusion between the derivative particle model and the linear model using the particle level set method, applicants used an effective resolution of 1282To perform a 2D dam break test in which a column of water of 0.2 x 0.4m is released in the gravitational field. The results are shown in FIG. 10, where applicants can see that derivative particles produce less diffusion: the wave breaks more sharply and the vorticity is preserved well.
Bulk water falling onto shallow water: FIG. 11 shows a screenshot taken from the following simulation: i.e. where the bulk of the water falls onto a reservoir 0.05m deep. Applicant uses a no-slip boundary condition for the bottom surface of the reservoir which causes the water to move quickly at the top and slowly at the bottom, resulting in a crown-like splash. 256 for this experiment3Is performed at the effective grid resolution. The same scenario is also simulated with a traditional linear model using the particle level set approach. The comparison demonstrates that the proposed technique yields more realistic results. The derivative particle model uses 60 seconds per time step on average, and the linear model uses 60 seconds per time step on averageAnd 30 seconds each including the file output time.
Impact of solid ball: in this experiment, solid spheres with a radius of 0.15m were impacted into water at a velocity (5.0, -3.0, 0.0) m/sec. This impact produces the complex structure shown in fig. 7. This experiment shows that the proposed technique is able to produce detailed fluid motion and surface features that occur in violent solid-water-air interactions. The effective resolution of this experiment was 192 × 128 × 128.
Rotating the water tank: figure 12 shows a 0.9 x 0.3m spin tank with half full of water. In the experiment, the centrifugal force produced the effect of pushing the water to the tank side. This simulation is performed with an effective grid resolution of 96 x 32. This experiment demonstrates that derivative particle models can simulate complex motion of a fluid even in a coarse resolution grid.
Flood: in this experiment, applicants simulated mesoscale flooding that occurred over a 6m × 2m × 4m domain. Fig. 13 shows a screenshot taken from this experiment. When the dike is removed, the water begins to slide off the steps. The water undergoes violent movements that produce complex geometries, due to obstructions along the path. The limited resolution of this experiment was 384 × 128 × 256.
FIGS. 1-6 show flow charts of the present invention.
1-6, a method for simulating detailed motion of a fluid using derivative particles includes the steps of: a) modeling the fluid with an adaptive mesh based on an octree data structure (S100); b) solving the navier-stokes equation for the incompressible fluid for the fluid velocity at the grid points (S200); c) using a step-by-step method in the time integration of the Navier-Stokes equation (S300); d) dividing the steps into a non-advection part and an advection part (S400); e) selecting a mesh advection part or a particle advection part according to the degree of detail of the simulation (S500); f) using an octree-based CIP method for the mesh advection part (S600); and g) using the derivative particle model for the particle advection part (S700).
The step of using the derivative particle model for the particle advection part (S700) is used when the degree of detail of the simulation is high.
The adaptive mesh increases the number of meshes for the details of the region of high interest.
The octree data structure comprises a tree data structure with up to eight children per internal node.
The Navier-Stokes equation for incompressible fluids can be written as
<math> <mrow> <mfrac> <mrow> <mo>&PartialD;</mo> <mi>u</mi> </mrow> <mrow> <mo>&PartialD;</mo> <mi>t</mi> </mrow> </mfrac> <mo>=</mo> <mo>-</mo> <mrow> <mo>(</mo> <mi>u</mi> <mo>&CenterDot;</mo> <mo>&dtri;</mo> <mo>)</mo> </mrow> <mi>u</mi> <mo>-</mo> <mfrac> <mrow> <mo>&dtri;</mo> <mi>p</mi> </mrow> <mi>&rho;</mi> </mfrac> <mo>+</mo> <mo>&dtri;</mo> <mo>&CenterDot;</mo> <mrow> <mo>(</mo> <mi>&mu;</mi> <mo>&dtri;</mo> <mi>u</mi> <mo>)</mo> </mrow> <mo>+</mo> <mfrac> <mi>f</mi> <mi>&rho;</mi> </mfrac> <mo>,</mo> </mrow></math>
<math> <mrow> <mo>&dtri;</mo> <mo>&CenterDot;</mo> <mi>u</mi> <mo>=</mo> <mn>0</mn> <mo>,</mo> </mrow></math>
Where u is velocity, p is pressure, f is external force, μ is viscosity coefficient, and ρ is fluid density.
The method may further include the step of tracking the surface of the fluid using a level set method (S800).
As shown in fig. 1 to 3, the level set method includes a level set field phi. Updating the level set field φ according to:
<math> <mrow> <mfrac> <mrow> <mo>&PartialD;</mo> <mi>&phi;</mi> </mrow> <mrow> <mo>&PartialD;</mo> <mi>t</mi> </mrow> </mfrac> <mo>+</mo> <mi>u</mi> <mo>&CenterDot;</mo> <mo>&dtri;</mo> <mi>&phi;</mi> <mo>=</mo> <mn>0</mn> </mrow></math>
distance of symbol <math> <mrow> <mo>|</mo> <mo>&dtri;</mo> <mi>&phi;</mi> <mo>|</mo> <mo>=</mo> <mn>1</mn> <mo>.</mo> </mrow></math>
<math> <mrow> <mfrac> <mrow> <mo>&PartialD;</mo> <mi>&phi;</mi> </mrow> <mrow> <mo>&PartialD;</mo> <mi>t</mi> </mrow> </mfrac> <mo>+</mo> <mi>u</mi> <mo>&CenterDot;</mo> <mo>&dtri;</mo> <mi>&phi;</mi> <mo>=</mo> <mn>0</mn> </mrow></math> The time integral of (a) involves only the advection part. The mesh advection part advects the velocity and level set fields computed from the non-advection part according to the current velocity field.
The grid advection part advects the delivery speed and level set using an Euler advection step. The euler advection is performed by an octree-based CIP (constrained interpolation profile) method.
The step of using the octree-based CIP method for the grid advection part (S600) includes the step of advecting the transport speed and level sets and their derivatives with the octree-based CIP method (S610). Direct from equation <math> <mrow> <mfrac> <mrow> <mo>&PartialD;</mo> <mi>&phi;</mi> </mrow> <mrow> <mo>&PartialD;</mo> <mi>t</mi> </mrow> </mfrac> <mo>+</mo> <mi>u</mi> <mo>&CenterDot;</mo> <mo>&dtri;</mo> <mi>&phi;</mi> <mo>=</mo> <mn>0</mn> </mrow></math> Obtaining an equation for advection transport derivatives, where φ is a physical quantity advected, and obtaining an advection transport equation for a derivative of a spatial variable x by differentiating the equation for x, i.e.:
<math> <mrow> <mfrac> <mrow> <mo>&PartialD;</mo> <msub> <mi>&phi;</mi> <mi>x</mi> </msub> </mrow> <mrow> <mo>&PartialD;</mo> <mi>t</mi> </mrow> </mfrac> <mo>+</mo> <mi>u</mi> <mo>&CenterDot;</mo> <mo>&dtri;</mo> <msub> <mi>&phi;</mi> <mi>x</mi> </msub> <mo>=</mo> <mo>-</mo> <msub> <mi>u</mi> <mi>x</mi> </msub> <mo>&CenterDot;</mo> <mo>&dtri;</mo> <mi>&phi;</mi> <mo>.</mo> </mrow></math>
as shown in fig. 2 and 3, the step of solving the advection equation of the derivative of x (S610) includes employing a step using a step method, and the step using the step method includes the steps of: a) solving for non-advective parts using finite difference method <math> <mrow> <mo>&PartialD;</mo> <msub> <mi>&phi;</mi> <mi>x</mi> </msub> <mo>/</mo> <mo>&PartialD;</mo> <mi>t</mi> <mo>=</mo> <mo>-</mo> <msub> <mi>u</mi> <mi>x</mi> </msub> <mo>&CenterDot;</mo> <mo>&dtri;</mo> <mi>&phi;</mi> </mrow></math> (S620); b) advective conveying device <math> <mrow> <mfrac> <mrow> <mo>&PartialD;</mo> <msub> <mi>&phi;</mi> <mi>x</mi> </msub> </mrow> <mrow> <mo>&PartialD;</mo> <mi>t</mi> </mrow> </mfrac> <mo>+</mo> <mi>u</mi> <mo>&CenterDot;</mo> <mo>&dtri;</mo> <msub> <mi>&phi;</mi> <mi>x</mi> </msub> <mo>=</mo> <mn>0</mn> </mrow></math> The result of (S630); and c) sampling pressure values at the centers of the cells defined by the grid and sampling other values including velocity, level sets and their derivatives at the nodes (S640).
As shown in fig. 3, the method may further comprise the steps of: a) combining an octree-based CIP method with an adaptive mesh (S650); and b) minimizing octree artifacts caused by regional variations in grid resolution of the octree data structure by simulating advection with the octree-based CIP method (S660).
The derivative particle model uses a particle-based lagrangian scheme. By using <math> <mrow> <mi>&phi;</mi> <mrow> <mo>(</mo> <mi>x</mi> <mo>)</mo> </mrow> <mo>=</mo> <msub> <mi>&phi;</mi> <mi>p</mi> </msub> <mo>+</mo> <mfrac> <mrow> <mo>&dtri;</mo> <msub> <mi>&phi;</mi> <mi>p</mi> </msub> </mrow> <mrow> <mo>|</mo> <mo>&dtri;</mo> <msub> <mi>&phi;</mi> <mi>p</mi> </msub> <mo>|</mo> </mrow> </mfrac> <mo>&CenterDot;</mo> <mrow> <mo>(</mo> <mi>x</mi> <mo>-</mo> <msub> <mi>x</mi> <mi>p</mi> </msub> <mo>)</mo> </mrow> </mrow></math> And <math> <mrow> <mo>&dtri;</mo> <mi>&phi;</mi> <mrow> <mo>(</mo> <mi>x</mi> <mo>)</mo> </mrow> <mo>=</mo> <mo>&dtri;</mo> <msub> <mi>&phi;</mi> <mi>p</mi> </msub> </mrow></math> to calculate the grid pointsWhere p is the gradient carried by the derivative particle, and its derivative.
The method may further comprise the steps of: a) converting the set of velocities and levels defined on the grid to a set of velocities and levels defined on the particles (S690) (grid to particles) prior to the step of using the derivative particle model for the particle advection part (S700); and b) converting the velocity and level set defined on the particles to a velocity and level set defined on the grid (S710) (particle to grid) before returning to the step of using octree-based CIP method for the grid advection part (S600).
Step of converting the velocity and level sets (S690) grid to particle includes using the FLIP method with monotonic CIP interpolation instead of linear interpolation for velocity <math> <mrow> <msup> <mi>u</mi> <mi>P</mi> </msup> <mo>=</mo> <msubsup> <mi>u</mi> <mo>-</mo> <mi>P</mi> </msubsup> <mo>+</mo> <msubsup> <mi>&Sigma;</mi> <mrow> <mi>i</mi> <mo>=</mo> <mn>1</mn> </mrow> <mn>4</mn> </msubsup> <msub> <mi>w</mi> <mi>i</mi> </msub> <mi>&Delta;</mi> <msubsup> <mi>u</mi> <mi>i</mi> <mi>G</mi> </msubsup> </mrow></math> And a step of conversion of uP (S692), wherein,
Figure A200680054876D00245
is the particle velocity just before the start of the non-advection part, and
Figure A200680054876D00246
is the velocity change at Gi due to the non-advective part.
As shown in fig. 5, for speed uPGrid points G, nearest particle P selected from each quadrant of G1、P2、P3And P4、PiVelocity of, and
Figure A200680054876D00247
Figure A200680054876D00248
and
Figure A200680054876D00249
the step of converting (S710) the velocity and level sets to a grid includes the steps of: a) will be provided withDerivative of (2)
Figure A200680054876D002411
Coordinate transformation to
Figure A200680054876D002412
Thus, it is possible to provide
Figure A200680054876D0024121550QIETU
Indicates along the direction
Figure A200680054876D002413
Component (a) of1⊥Is shown perpendicular to
Figure A200680054876D00251
Wherein similarly obtained
Figure A200680054876D00252
Derivative of the coordinate transformation of
Figure A200680054876D00253
(S712); b) along the direction
Figure A200680054876D00254
Performing one-dimensional monotonic CIP interpolation on the result obtained in step (a) to calculate u at position AAAnd it
Figure A200680054876D00255
Derivative of direction
Figure A200680054876D00256
(S714); c) by applying the inverse of the above coordinate transformation
Figure A200680054876D00257
Obtaining (u)Ax,uAy) (S716); d) to particle P3And P4Performing steps (a) - (c) to calculate u of BBAnd (a)uBx,uBy) (S718); and e) performing y-direction monotonic CIP interpolation on the results obtained in steps (c) and (d), which gives the velocity and derivative at G (S720).
As shown in FIG. 6, the step of using the derivative particle model for the particle advection part comprises the steps of: a) adjusting a level set value of the particle (S702); and b) performing particle reseeding (S704).
7 conclusion
In the present invention, applicants have introduced a new concept called derivative particles and proposed fluid simulation techniques that increase the accuracy of the advection part. The non-advection part of the simulator is implemented using the euler scheme, while the advection part is implemented using the lagrangian scheme. An important issue in developing a simulator that incorporates in this way is how to convert the physical quantities at the switching of the advection and non-advection parts. Applicants have successfully overcome this problem by developing conversion procedures that effectively suppress the unwanted dissipation of values. In fluid regions where particle advection need not be performed, advection was simulated using octree-based CIP solvers, which applicants developed in this work by combining CIP interpolation with an octree data structure. This method also serves to reduce numerical dissipation. The results of multiple experiments demonstrate that the proposed technique significantly reduces velocity dissipation, resulting in a true reproduction of the dynamic fluid.
This work would not be possible without the pioneering work of previous researchers. The idea of using particles entirely to solve the advection part comes from PIC [ Harlow 1963] and FLIP [ Brackbill and Ruppel 1986] methods, but for particle-to-grid and grid-to-particle velocity scaling, applicant uses derivative-based cubic interpolation instead of linear interpolation. In addition, the idea of letting the particles carry a level set value comes from the particle level set method [ Enright et al 2002b ]. This model is extended in this work so that the particles also carry velocity information. In addition, derivative information is used from CIP methods [ Yabe et al 2001 ]. Here, the applicant extended the application of this method to particles, leading to the development of a particle-to-grid-speed conversion program that contains extensions to the CIP method itself: CIP interpolation of particles in a non-rectangular configuration.
In the experiments, the ability to reproduce details was emphasized. However, the proposed technique can also be applied to large-scale fluid movements to produce accurate results. This technique involves an increased amount of memory for storing derivative information. However, applicants have demonstrated that 2563 grids can be simulated on current PCs, which is already a useful resolution for delineating fluid scenarios of interest.
While the invention has been shown and described with reference to various embodiments, those skilled in the art will recognize that changes may be made in form, detail, composition and operation without departing from the spirit and scope of the invention, as defined by the claims.
Reference to the literature
BRACKBILL,J.U.,AND RUPPEL,H.M.1986.Flip:A method foradaptively zoned,particle-in-cell calculations of fluid flows in twodimensions.J.Comp.Phys.65,314-343.
CARLSON,M.,MUCHA,R.J.,AND TURK,G.2004.Rigid fluid:Animating the interplay between rigid bodies and fluid.ACM Transactionson Graphics 23,3,377-384.
ENRIGHT,D.,FEDKIW,R.,FERZIGER,J.,AND MITCHELL,I.2002.A hybrid particle level set method for improved interface capturing.J.Comp.Phys.183,83-116.
ENRIGHT,D.,MARSCHNER,S.,AND FEDKIW,R.2002.Animation and rendering of complex water surfaces.ACM Transactionson Graphics 21,3,736-744.
ENRIGHT,D.,LOSASSO,F.,AND FEDKIW,R.2005.A fast andaccurate semi-lagrangian particle level′set method.Computers andStructures 83,479-490.
FEDKIW,R.,STAM,J.,AND JENSEN,H.W.2001.Visualsimulation of smoke.Computer Graphics(Proc.ACM SIGGRAPH 2001)35,15-22.
FELDMAN,B.E.,O′BRIEN,J.F.,AND ARIKAN,0.2003.Animating suspended particle explosions.ACM Transactions on Graphics22,3,708-715.
FELDMAN,B.E.,O′BRIEN,J.F.,AND KLINGNER,B.M.2005.Animating gases with hybrid meshes.ACM Transactions on Graphics 24,3,904-909.
FELDMAN,B.E.,O′BRIEN,J.F.,KLINGNER,B.M.,ANDGOKTEKIN,T.G.2005.Fluids in deforming meshes.In Proceedings ofEurographics/ACM SIGGRAPH Symposium on Computer Animation2005,255-259.
FOSTER,N.,AND FEDKIW,R.2001.Practical animation of liquids.Computer Graphics(Proc.ACM SIGGRAPH 2001)35,23-30.
FOSTER,N.,AND METAXAS,D.1996.Realistic animation ofliquids.Graphical models and image processing:GMIP 58,5,471-483.
FOSTER,N.,AND METAXAS,D.1997.Controlling fluidanimation.In Computer Graphics International 97,178-188.
GOKTEKIN,T.G.,BARGTEIL,A.W.,AND O′BRIEN,J.F.2004.A method for animating viscoelastic fluids.ACMT ransactions onGraphics 23,3,463-468.
GUENDELMAN,E.,SELLE,A.,LOSASSO,F.,AND FEDKIW,R.2005.Coupling water and smoke to thin deformable and rigid shells.ACMTransactions on Graphics 24,3,973-981.
HARLOW,F.H.1963.The particle-in-cell method for numericalsolution of problems in fluid dynamics.In Experimental arithmetic,high-speed computations and mathematics.
HONG,J.-M.,AND KIM,C-H.2005.Discontinuous fluids.ACMTransactions on Graphics 24,3,915-920.
KANG,M.,FEDKIW,R.,AND LIU,X.-D.2000.A boundary-condition capturing method for multiphase incompressible flow.J.Sci.Comput.15,323-360.
LIU,X.-D.,FEDKIW,R.,AND KANG,M.2000.A boundarycondition capturing method for poisson′s equation on irregular domains.J.Comp.Phys.160,151-178.
LOSASSO,F.,GIBOU,F.,AND FEDKIW,R.2004.Simulatingwater and smoke with an octree data structure.ACM Transactions onGraphics 23,3,457-462.
MCNAMARA,A.,TREUILLE,A.,POPOVI′C,Z.,AND STAM,J.2004.Fluid control using the adjoint method.ACM Transactions onGraphics 23,3,449-456.
M"ULLER,M.,CHARYPAR,D.,AND GROSS,M.2003.Particlebased fluid simulation for interactive applications.In Proceedingsof ACM SIGGRAPH/Eurographics Symposium on Computer Animation2003,154-159.
M"ULLER,M.,SOLENTHALER,B.,KEISER,R.,AND GROSS,M.2005.Particle-based fluid-fluid interaction.In Proceedings ofEurographics/ACM SIGGRAPH Symposium on Computer Animation2005,237-244.
PARK,S.I.,AND KIM,M.J.2005.Vortex fluid for gaseousphenomena.In Proceedings of Eurographics/ACM SIGGRAPHSymposium on Computer Animation 2005,261-270.
PREMCTZ E,S.,TASDIZEN,T.,BIGLER,J.,LEFOHN,A.,ANDWHITAKER,R.T.2003.Particle-based simulation of fluids.InEurographics 2003 proceedings,Blackwell Publishers,401-410.
RASMUSSEN,N.,NGUYEN,D.Q.,GEIGER,W.,AND FEDKIW,R.2003.Smoke simulation for large scale phenomena.ACM Transactionson Graphics 22,3,703-707.
SELLE,A.,RASMUSSEN,N.,AND FEDKIW,R.2005.A vortexparticle method for smoke,water and explosions.ACM Transactions onGraphics 24,3,910-914.
SONG,O.-Y.,SHIN,H.,AND KO,H.-S.2005.Stable butnondissipative water.ACM Transactions on Graphics 24,1,81-97.
STAM,J.,AND FIUME,E.1995.Depicting fire and other gaseousphenomena using diffusion processes.Computer Graphics(Proc.ACMSIGGRAPH′95)29,Annual Conference Series,129-136.
STAM,J.1999.Stable fluids.Computer Graphics(Proc.ACMSIGGRAPH′99)33,Annual Conference Series,121-128.
STANIFORTH,A.,AND CλO TTE,J.1991.Semi-lagrangianintegration scheme for atmospheric model-a review.Mon.Weather Rev.119,12,2206-2223.
TREUILLE,A.,MCNAMARA,A.,POPOVI′C,Z.,AND STAM,J.2003.Keyframe control of smoke simulations.ACM Transactions onGraphics 22,3,716-723.
WANG,H.,MUCHA,P.J.,AND TURK,G.2005.Water drops onsurfaces.ACM Transactions on Graphics 24,3,921-929.
YABE,T.,AND AOKI,T.1991.A universal solver for hyperbolicequations by cubic-polynomial interpolation i.one-dimensional solver.Comp.Phys.Comm.66,219-232.
YABE,T.,AND AOKI,T.1991.A universal solver for hyperbolicequations by cubic-polynomial interpolation ii.two-and three dimensionalsolvers.Comp.Phys.Comm.66,233-242.
YABE,T.,XIAO,F.,AND UTSUMI,T.2001.The constrainedinterpolation profile method for multiphase analysis.J.Comp.Phys.169,556-593.
ZHU,Y.,AND BRIDSON,R.2005.Animating sand as a fluid.ACMTransactions on Graphics 24,3,965-972.

Claims (19)

1. A method of simulating detailed movement of a fluid using derivative particles, comprising the steps of:
a) modeling the fluid with an adaptive mesh based on an octree data structure;
b) solving a Navier-Stokes equation for the incompressible fluid for the fluid velocity at the grid points;
c) using a step-wise method in the time integration of the Navier-Stokes equation;
d) dividing each step into a non-advection part and an advection part;
e) selecting a grid advection part or a particle advection part according to the detail degree of the simulation;
f) using an octree-based CIP method for the mesh advection part; and
g) a derivative particle model is used for the particle advection part,
wherein, when the degree of detail of the simulation is high, the step of using the derivative particle model for the mass point advection part is used.
2. The method of claim 1, wherein the adaptive mesh increases the number of meshes for regions of high interest in detail.
3. The method of claim 1, wherein the octree data structure comprises a tree data structure, wherein each internal node has up to eight children.
4. The method of claim 1, wherein the Navier-Stokes equation for incompressible fluids is written as
<math> <mrow> <mfrac> <mrow> <mo>&PartialD;</mo> <mi>u</mi> </mrow> <mrow> <mo>&PartialD;</mo> <mi>t</mi> </mrow> </mfrac> <mo>=</mo> <mo>-</mo> <mrow> <mo>(</mo> <mi>u</mi> <mo>&CenterDot;</mo> <mo>&dtri;</mo> <mo>)</mo> </mrow> <mi>u</mi> <mo>-</mo> <mfrac> <mrow> <mo>&dtri;</mo> <mi>p</mi> </mrow> <mi>&rho;</mi> </mfrac> <mo>+</mo> <mo>&dtri;</mo> <mo>&CenterDot;</mo> <mrow> <mo>(</mo> <mi>&mu;</mi> <mo>&dtri;</mo> <mi>u</mi> <mo>)</mo> </mrow> <mo>+</mo> <mfrac> <mi>f</mi> <mi>&rho;</mi> </mfrac> <mo>,</mo> </mrow></math>
<math> <mrow> <mo>&dtri;</mo> <mo>&CenterDot;</mo> <mi>u</mi> <mo>=</mo> <mn>0</mn> </mrow></math>
Where u is velocity, p is pressure, f is external force, μ is viscosity coefficient, and ρ is fluid density.
5. The method of claim 1, further comprising the step of tracking the surface of the fluid using a level set method, wherein the level set method comprises a level set field φ.
6. The method of claim 5, wherein the level set field φ | is updated according to
<math> <mrow> <mfrac> <mrow> <mo>&PartialD;</mo> <mi>&phi;</mi> </mrow> <mrow> <mo>&PartialD;</mo> <mi>t</mi> </mrow> </mfrac> <mo>+</mo> <mi>u</mi> <mo>&CenterDot;</mo> <mo>&dtri;</mo> <mi>&phi;</mi> <mo>=</mo> <mn>0</mn> </mrow></math>
With signed distance <math> <mrow> <mrow> <mo>|</mo> <mo>&dtri;</mo> <mi>&phi;</mi> <mo>|</mo> </mrow> <mo>=</mo> <mn>1</mn> <mo>.</mo> </mrow></math>
7. The method of claim 6, wherein the <math> <mrow> <mfrac> <mrow> <mo>&PartialD;</mo> <mi>&phi;</mi> </mrow> <mrow> <mo>&PartialD;</mo> <mi>t</mi> </mrow> </mfrac> <mo>+</mo> <mi>u</mi> <mo>&CenterDot;</mo> <mo>&dtri;</mo> <mi>&phi;</mi> <mo>=</mo> <mn>0</mn> </mrow></math> The time integral of (a) involves only the advection part.
8. The method of claim 6, wherein the mesh advection portion advects the velocity and level set fields computed from the non-advection portion according to a current velocity field.
9. The method of claim 6, wherein the grid advection part advects delivery speed and level sets using an Euler advection step.
10. The method according to claim 9, wherein the euler advection is performed by octree-based CIP (constrained interpolation profile) method.
11. The method of claim 1, wherein the step of using octree-based CIP method for the advection part of the grid comprises the step of advecting the set of velocities and levels and their derivatives with the octree-based CIP method.
12. The method of claim 11, wherein the equation is directly derived from <math> <mrow> <mfrac> <mrow> <mo>&PartialD;</mo> <mi>&phi;</mi> </mrow> <mrow> <mo>&PartialD;</mo> <mi>t</mi> </mrow> </mfrac> <mo>+</mo> <mi>u</mi> <mo>&CenterDot;</mo> <mo>&dtri;</mo> <mi>&phi;</mi> <mo>=</mo> <mn>0</mn> </mrow></math> Obtaining an equation for advection transport derivatives, where φ is a physical quantity advected, and obtaining an advection transport equation for a derivative of a spatial variable x by differentiating the equation with respect to x, i.e.:
<math> <mrow> <mfrac> <msub> <mrow> <mo>&PartialD;</mo> <mi>&phi;</mi> </mrow> <mi>x</mi> </msub> <mrow> <mo>&PartialD;</mo> <mi>t</mi> </mrow> </mfrac> <mo>+</mo> <mi>u</mi> <mo>&CenterDot;</mo> <msub> <mrow> <mo>&dtri;</mo> <mi>&phi;</mi> </mrow> <mi>x</mi> </msub> <mo>=</mo> <msub> <mrow> <mo>-</mo> <mi>u</mi> </mrow> <mi>x</mi> </msub> <mo>&CenterDot;</mo> <mo>&dtri;</mo> <mi>&phi;</mi> <mo>.</mo> </mrow></math>
13. the method of claim 12, wherein the step of solving the advection equation for the derivative of x comprises the step of using a step method, wherein the step of using a step method comprises the steps of:
a) solving for non-advective parts using finite difference method <math> <mrow> <msub> <mrow> <mo>&PartialD;</mo> <mi>&phi;</mi> </mrow> <mi>x</mi> </msub> <mo>/</mo> <mo>&PartialD;</mo> <mi>t</mi> <mo>=</mo> <mo>-</mo> <msub> <mi>u</mi> <mi>x</mi> </msub> <mo>&CenterDot;</mo> <mo>&dtri;</mo> <mi>&phi;</mi> <mo>;</mo> </mrow></math>
b) Advective conveying device <math> <mrow> <mfrac> <msub> <mrow> <mo>&PartialD;</mo> <mi>&phi;</mi> </mrow> <mi>x</mi> </msub> <mrow> <mo>&PartialD;</mo> <mi>t</mi> </mrow> </mfrac> <mo>+</mo> <mi>u</mi> <mo>&CenterDot;</mo> <mo>&dtri;</mo> <msub> <mi>&phi;</mi> <mi>x</mi> </msub> <mo>=</mo> <mn>0</mn> </mrow></math> The result of (1); and
c) the pressure values are sampled at the center of the cells defined by the grid and other values including velocity, level sets and their derivatives are sampled at each node.
14. The method of claim 1, further comprising the step of:
a) combining an octree-based CIP method with an adaptive grid; and
b) OctreeImage artifacts caused by regional variations in the grid resolution of the octree data structure are minimized by simulating advection with an octree-based CIP method.
15. The method as in claim 1, wherein the derivative particle model uses a particle-based Lagrangian scheme, wherein the level set at grid points and their derivatives are calculated with the following equation
<math> <mrow> <mi>&phi;</mi> <mrow> <mo>(</mo> <mi>x</mi> <mo>)</mo> </mrow> <mo>=</mo> <msub> <mi>&phi;</mi> <mi>p</mi> </msub> <mo>+</mo> <mfrac> <msub> <mrow> <mo>&dtri;</mo> <mi>&phi;</mi> </mrow> <mi>p</mi> </msub> <mrow> <mo>(</mo> <msub> <mrow> <mo>&dtri;</mo> <mi>&phi;</mi> </mrow> <mi>p</mi> </msub> <mo>)</mo> </mrow> </mfrac> <mo>&CenterDot;</mo> <mrow> <mo>(</mo> <mi>x</mi> <mo>-</mo> <msub> <mi>x</mi> <mi>p</mi> </msub> <mo>)</mo> </mrow> </mrow></math>
And
<math> <mrow> <mo>&dtri;</mo> <mi>&phi;</mi> <mrow> <mo>(</mo> <mi>x</mi> <mo>)</mo> </mrow> <mo>=</mo> <msub> <mrow> <mo>&dtri;</mo> <mi>&phi;</mi> </mrow> <mi>p</mi> </msub> <mo>,</mo> </mrow></math>
where p is the gradient carried by the derivative particle.
16. The method of claim 1, further comprising the steps of:
a) converting the set of velocities and levels defined on each grid to a set of velocities and levels defined on each particle (grid to particles) prior to said step of using the derivative particle model for the particle advection part; and
b) the velocity and level sets defined on each particle are converted to the velocity and level sets defined on each grid (particle to grid) before returning to the step of using the octree-based CIP method on the grid advection part.
17. The method of claim 16, wherein the step of grid to particle conversion speed and level set comprises for speed uPUsing the steps of the FLIP method with monotonic CIP interpolation instead of linear interpolation, wherein <math> <mrow> <msup> <mi>u</mi> <mi>P</mi> </msup> <mo>=</mo> <msubsup> <mi>u</mi> <mo>-</mo> <mi>P</mi> </msubsup> <mo>+</mo> <msubsup> <mi>&Sigma;</mi> <mrow> <mi>i</mi> <mo>-</mo> <mn>1</mn> </mrow> <mn>4</mn> </msubsup> <msub> <mi>w</mi> <mi>i</mi> </msub> <msubsup> <mi>&Delta;u</mi> <mi>i</mi> <mi>G</mi> </msubsup> <mo>,</mo> </mrow></math> WhereinIs the particle velocity just before the onset of the non-advection part, and
Figure A200680054876C0005135500QIETU
is due to a non-advective part of GiThe speed of (c) is changed.
18. The method of claim 16, wherein for speed u, uPGrid points G, nearest particle P selected from each quadrant of G1、P2、P3And P4The speed of Pi, and
Figure A200680054876C00053
Figure A200680054876C00054
and
Figure A200680054876C00055
the derivatives of the x-component and the y-component of (c),
the step of converting the particle to grid speed and level set comprises the steps of:
a) will be provided with
Figure A200680054876C00056
Derivative of (2)
Figure A200680054876C00057
Coordinate transformation into (u)1‖u1⊥) So that u is1‖Indicates along the direction
Figure A200680054876C00058
Component (a) of1⊥Is shown perpendicular to
Figure A200680054876C00059
Wherein similarly obtained
Figure A200680054876C000510
Derivative (u) of the coordinate transformation of (c)2‖,u2⊥);
b) Along the direction
Figure A200680054876C000511
Performing one-dimensional monotonic CIP interpolation on the result obtained in step (a) to calculate u at position AAAnd it
Figure A200680054876C000512
Derivative of direction (u)A‖,uA⊥);
c) By applying the inverse of the above-described coordinate transformation,from (u)A‖,uA⊥) Obtaining (u)Ax,uAy);
d) To particle P3And P4Performing steps (a) - (c) to calculate u of BBAnd (u)Bx,uBy) B is the grid line and the connection P3And P4The intersection points between the line segments of (a); and
e) performing y-direction monotonic CIP interpolation on the results obtained in steps (c) and (d), giving velocity and derivatives at G.
19. The method of claim 15, wherein the step of using a derivative particle model for the mass advection part comprises the steps of:
a) adjusting a level set value of the particle; and
b) particle reseeding is performed.
CNA2006800548768A 2006-04-05 2006-10-27 Method of simulating detailed movements of fluids using derivative particles Pending CN101484922A (en)

Applications Claiming Priority (2)

Application Number Priority Date Filing Date Title
US11/398,982 2006-04-05
US11/398,982 US7565276B2 (en) 2006-04-05 2006-04-05 Method of simulating detailed movements of fluids using derivative particles

Publications (1)

Publication Number Publication Date
CN101484922A true CN101484922A (en) 2009-07-15

Family

ID=38563801

Family Applications (1)

Application Number Title Priority Date Filing Date
CNA2006800548768A Pending CN101484922A (en) 2006-04-05 2006-10-27 Method of simulating detailed movements of fluids using derivative particles

Country Status (4)

Country Link
US (1) US7565276B2 (en)
EP (1) EP2008250A4 (en)
CN (1) CN101484922A (en)
WO (1) WO2007114548A1 (en)

Cited By (5)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN102147932A (en) * 2011-03-30 2011-08-10 北京航空航天大学 Method for simulating smog driven by movable Euler grid-based model
CN102203782A (en) * 2010-09-09 2011-09-28 天津空中代码工程应用软件开发有限公司 Numerical method for solving eulerian equation of lagrange type
CN103473418A (en) * 2013-09-18 2013-12-25 梧州学院 Simulation method of pressure difference in low Reynolds number incompressible flow at bending boundaries
CN107133418A (en) * 2017-05-26 2017-09-05 刘哲 Fluids material stratospheric transport analogy method based on alternately TVD difference algorithms
CN110457798A (en) * 2019-07-29 2019-11-15 广东工业大学 A kind of adaptive vorticity restraint method based on vorticity loss

Families Citing this family (68)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
GB2442496A (en) * 2006-10-02 2008-04-09 Univ Sheffield Evaluating displacements at discontinuities within a body
US7921003B2 (en) * 2007-01-23 2011-04-05 Adobe Systems Incorporated System and method for simulating shallow water effects on arbitrary surfaces
WO2009035737A2 (en) * 2007-06-13 2009-03-19 United States Of America As Represented By The Administrator Of The National Aeronautics And Space Administration Adaptive refinement tools for tetrahedral unstructured grids
KR100872434B1 (en) * 2007-10-25 2008-12-05 한국전자통신연구원 System and method for simulating particle fluid having multi-resolution
CA2648441A1 (en) * 2007-12-31 2009-06-30 Exocortex Technologies, Inc. Fast characterization of fluid dynamics
FR2930350B1 (en) * 2008-04-17 2011-07-15 Inst Francais Du Petrole PROCESS FOR SEARCHING FOR HYDROCARBONS IN A GEOLOGICALLY COMPLEX BASIN USING BASIN MODELING
CN101329772B (en) * 2008-07-21 2010-06-02 北京理工大学 Emulation modelling method interacted with movable object and water based on SPH
JP5268496B2 (en) * 2008-08-22 2013-08-21 株式会社東芝 Flow analysis method, flow analysis apparatus, and flow analysis program
KR101159163B1 (en) * 2008-12-10 2012-06-22 한국전자통신연구원 Method, recording medium, and apparatus for incompressible fluid simulations
US20100185420A1 (en) * 2009-01-18 2010-07-22 Ejiang Ding Computer system for computing the motion of solid particles in fluid
US8289327B1 (en) * 2009-01-21 2012-10-16 Lucasfilm Entertainment Company Ltd. Multi-stage fire simulation
US8335675B1 (en) 2009-02-27 2012-12-18 Adobe Systems Incorporated Realistic real-time simulation of natural media paints
US8229719B2 (en) * 2009-03-26 2012-07-24 Seiko Epson Corporation Finite element algorithm for solving a fourth order nonlinear lubrication equation for droplet evaporation
US8055490B2 (en) * 2009-03-31 2011-11-08 Seoul National University Semi-Lagrangian CIP fluid solver without dimensional splitting
US8219370B1 (en) * 2009-05-20 2012-07-10 Adobe Systems Incorporated Simulation of shallow viscoelastic flows
US8711140B1 (en) * 2009-06-01 2014-04-29 Paradigm Sciences Ltd. Systems and methods for building axes, co-axes and paleo-geographic coordinates related to a stratified geological volume
US8600708B1 (en) 2009-06-01 2013-12-03 Paradigm Sciences Ltd. Systems and processes for building multiple equiprobable coherent geometrical models of the subsurface
US9418182B2 (en) 2009-06-01 2016-08-16 Paradigm Sciences Ltd. Systems and methods for building axes, co-axes and paleo-geographic coordinates related to a stratified geological volume
US8014986B2 (en) * 2009-06-02 2011-09-06 Seiko Epson Corporation Finite difference algorithm for solving lubrication equations with solute diffusion
US8285530B2 (en) * 2009-10-15 2012-10-09 Seiko Epson Corporation Upwind algorithm for solving lubrication equations
US8743115B1 (en) 2009-10-23 2014-06-03 Paradigm Sciences Ltd. Systems and methods for coordinated editing of seismic data in dual model
US8386224B2 (en) * 2009-11-06 2013-02-26 Doyub Kim Method for simulating stretching and wiggling liquids
US8255194B2 (en) * 2009-12-02 2012-08-28 Seiko Epson Corporation Judiciously retreated finite element method for solving lubrication equation
US8285526B2 (en) * 2009-12-02 2012-10-09 Seiko Epson Corporation Finite difference algorithm for solving slender droplet evaporation with moving contact lines
KR101269553B1 (en) * 2009-12-18 2013-06-04 한국전자통신연구원 Fire simulation method with particle fuel
WO2011093541A1 (en) * 2010-01-28 2011-08-04 (주)에프엑스기어 Fluid simulation shape control system and method
KR101170909B1 (en) 2010-01-28 2012-08-03 (주)에프엑스기어 System and method for fluid simulation using moving grid
US20110196657A1 (en) * 2010-02-11 2011-08-11 Jie Zhang Solving a Solute Lubrication Equation for 3D Droplet Evaporation on a Complicated OLED Bank Structure
US8271238B2 (en) * 2010-03-23 2012-09-18 Seiko Epson Corporation Finite difference scheme for solving droplet evaporation lubrication equations on a time-dependent varying domain
KR101113301B1 (en) * 2010-03-23 2012-02-24 한국과학기술원 Method and recording medium of a hybrid approach to multiple fluid simulation using volume fraction
US8725476B1 (en) * 2010-05-04 2014-05-13 Lucasfilm Entertainment Company Ltd. Applying details in a simulation
US10321892B2 (en) * 2010-09-27 2019-06-18 Siemens Medical Solutions Usa, Inc. Computerized characterization of cardiac motion in medical diagnostic ultrasound
US20120143573A1 (en) * 2010-12-01 2012-06-07 Hyeong-Seok Ko Method for Tracking Detail-Preserving Fully-Eulerian Interface
US8457939B2 (en) * 2010-12-30 2013-06-04 Aerion Corporation Generating inviscid and viscous fluid-flow simulations over an aircraft surface using a fluid-flow mesh
US8538738B2 (en) 2011-03-22 2013-09-17 Aerion Corporation Predicting transition from laminar to turbulent flow over a surface
US8892408B2 (en) 2011-03-23 2014-11-18 Aerion Corporation Generating inviscid and viscous fluid flow simulations over a surface using a quasi-simultaneous technique
US8970592B1 (en) 2011-04-19 2015-03-03 Lucasfilm Entertainment Company LLC Simulating an arbitrary number of particles
US8744812B2 (en) * 2011-05-27 2014-06-03 International Business Machines Corporation Computational fluid dynamics modeling of a bounded domain
US9984489B2 (en) 2011-07-27 2018-05-29 Dreamworks Animation L.L.C. Fluid dynamics framework for animated special effects
JP6220797B2 (en) * 2012-03-08 2017-10-25 エンジン シミュレーション パートナーズ Fluid dynamics system boundaries
CN102800116B (en) * 2012-06-18 2014-11-05 浙江大学 Method for rapidly creating large-scale virtual crowd
KR102205845B1 (en) * 2013-10-28 2021-01-22 삼성전자주식회사 Method and apparatus for modeling based on particles
EP2869096B1 (en) 2013-10-29 2019-12-04 Emerson Paradigm Holding LLC Systems and methods of multi-scale meshing for geologic time modeling
KR101514806B1 (en) * 2014-01-20 2015-04-24 동국대학교 산학협력단 Method and apparatus of fluid simulation
US20160004802A1 (en) * 2014-07-03 2016-01-07 Arizona Board Of Regents On Behalf Of Arizona State University Multiscale Modelling of Growth and Deposition Processes in Fluid Flow
US9934339B2 (en) 2014-08-15 2018-04-03 Wichita State University Apparatus and method for simulating machining and other forming operations
CN104182598A (en) * 2014-09-18 2014-12-03 重庆大学 Constraint damping structure optimizing and designing method based on level set method
KR102263096B1 (en) * 2014-11-13 2021-06-09 삼성전자주식회사 Method and apparatus for modeling objects consisting of particles
US9690002B2 (en) 2015-06-18 2017-06-27 Paradigm Sciences Ltd. Device, system and method for geological-time refinement
CN105279781B (en) * 2015-10-23 2018-06-08 山东师范大学 Fluid animation generation method based on the fusion of more precision
US9881110B1 (en) * 2015-10-29 2018-01-30 Sohrab Mohajerin Apparatus and method for estimating and modeling turbulent flow
CN105760588B (en) * 2016-02-04 2022-02-25 自然资源部第一海洋研究所 SPH fluid surface reconstruction method based on two-layer regular grid
CN105841777B (en) * 2016-03-17 2019-04-30 深圳大学 A kind of multibeam echosounding estimation method and system based on adaptively selected node
CN106019930A (en) * 2016-08-03 2016-10-12 中国人民解放军63821部队 Aerodynamic/control integrated coupling simulating technology in aircraft maneuvering process
CN106327524B (en) * 2016-08-31 2019-04-02 上海交通大学 A kind of rapid fluid imaging surface method for tracing
US10466388B2 (en) 2016-09-07 2019-11-05 Emerson Paradigm Holding Llc System and method for editing geological models by switching between volume-based models and surface-based structural models augmented with stratigraphic fiber bundles
CN108269299B (en) * 2017-01-04 2021-07-27 北京航空航天大学 Viscous fluid modeling method based on approximate solution of SPH (particle-spray-drying) method
CN109359312B (en) * 2018-08-01 2022-11-15 中国科学院软件研究所 Real-time simulation and interaction method for dry particle flow
JP7012944B2 (en) * 2018-10-11 2022-01-31 オムロン株式会社 Simulation equipment, simulation method and simulation program
US10520644B1 (en) 2019-01-10 2019-12-31 Emerson Paradigm Holding Llc Imaging a subsurface geological model at a past intermediate restoration time
US11156744B2 (en) 2019-01-10 2021-10-26 Emerson Paradigm Holding Llc Imaging a subsurface geological model at a past intermediate restoration time
CN110059363A (en) * 2019-03-25 2019-07-26 天津大学 A method of fluid-mixing Phase transition simulation and liquid level reconstruct based on SPH
JP7160752B2 (en) * 2019-04-25 2022-10-25 株式会社日立製作所 Particle behavior simulation method and particle behavior simulation system
CN110457791B (en) * 2019-07-26 2023-03-28 嘉兴学院 Hydromechanical finite element algorithm based on characteristic line of Longge-Kutta time dispersion
WO2021101548A1 (en) * 2019-11-21 2021-05-27 Nvidia Corporation Fluid simulations using one or more neural networks
CN110955998B (en) * 2019-11-28 2023-04-25 青岛科技大学 GIS-based large-scale debris flow numerical simulation and numerical processing method
CN115698672B (en) * 2020-06-01 2024-05-31 住友金属矿山株式会社 Simulation device, simulation method, and storage medium
CN116303483B (en) * 2023-05-23 2023-07-21 北京适创科技有限公司 Compression method and device for structured grid, electronic equipment and storage medium

Family Cites Families (2)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
US7379852B2 (en) * 2004-02-18 2008-05-27 Chevron U.S.A. Inc. N-phase interface tracking method utilizing unique enumeration of microgrid cells
US7479963B2 (en) * 2004-05-14 2009-01-20 Yissum Research Development Company Of The Hebrew University Of Jerusalem Method and system for performing computer graphic simulation of a fluid using target-driven control

Cited By (7)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN102203782A (en) * 2010-09-09 2011-09-28 天津空中代码工程应用软件开发有限公司 Numerical method for solving eulerian equation of lagrange type
CN102147932A (en) * 2011-03-30 2011-08-10 北京航空航天大学 Method for simulating smog driven by movable Euler grid-based model
CN103473418A (en) * 2013-09-18 2013-12-25 梧州学院 Simulation method of pressure difference in low Reynolds number incompressible flow at bending boundaries
CN103473418B (en) * 2013-09-18 2017-02-08 梧州学院 Simulation method of pressure difference in low Reynolds number incompressible flow at bending boundaries
CN107133418A (en) * 2017-05-26 2017-09-05 刘哲 Fluids material stratospheric transport analogy method based on alternately TVD difference algorithms
CN107133418B (en) * 2017-05-26 2020-04-21 刘哲 Earth fluid material advection transportation simulation method based on alternative TVD differential algorithm
CN110457798A (en) * 2019-07-29 2019-11-15 广东工业大学 A kind of adaptive vorticity restraint method based on vorticity loss

Also Published As

Publication number Publication date
US20070239414A1 (en) 2007-10-11
EP2008250A4 (en) 2010-04-21
EP2008250A1 (en) 2008-12-31
US7565276B2 (en) 2009-07-21
WO2007114548A1 (en) 2007-10-11

Similar Documents

Publication Publication Date Title
CN101484922A (en) Method of simulating detailed movements of fluids using derivative particles
Stam Real-time fluid dynamics for games
US7647214B2 (en) Method for simulating stable but non-dissipative water
Kelager Lagrangian fluid dynamics using smoothed particle hydrodynamics
Song et al. Stable but nondissipative water
Selle et al. A vortex particle method for smoke, water and explosions
Losasso et al. Two-way coupled SPH and particle level set fluid simulation
Thürey et al. Free Surface Lattice-Boltzmann fluid simulations with and without level sets.
Thürey et al. Animation of open water phenomena with coupled shallow water and free surface simulations
US8055490B2 (en) Semi-Lagrangian CIP fluid solver without dimensional splitting
CN110717269A (en) Fluid surface detail protection method based on grid and particle coupling
Yang et al. Pairwise force SPH model for real-time multi-interaction applications
Qiu Two-dimensional SPH simulations of landslide-generated water waves
Yang et al. Versatile interactions at interfaces for SPH-based simulations.
Fan et al. Adapted unstructured LBM for flow simulation on curved surfaces
Shao et al. Particle‐based simulation of bubbles in water–solid interaction
Song et al. Derivative particles for simulating detailed movements of fluids
Maljaars et al. A numerical wave tank using a hybrid particle-mesh approach
Pan et al. Wake synthesis for shallow water equation
Zaspel et al. Photorealistic visualization and fluid animation: coupling of Maya with a two-phase Navier-Stokes fluid solver
Shi et al. An advanced hybrid smoothed particle hydrodynamics–fluid implicit particle method on adaptive grid for condensation simulation
Baek et al. Muddy water animation with different details
CN113673140B (en) Real-time interactive visual simulation method for fluid particles under action of air pressure
Coquerelle et al. A vortex method for bi-phasic fluids interacting with rigid bodies
Lee et al. Simulation of swirling bubbly water using bubble particles

Legal Events

Date Code Title Description
C06 Publication
PB01 Publication
C10 Entry into substantive examination
SE01 Entry into force of request for substantive examination
C12 Rejection of a patent application after its publication
RJ01 Rejection of invention patent application after publication

Application publication date: 20090715