这是一篇关于SPH 3D流水模拟的论文，看的关于这方面文章中最赞的一篇！希望给大家带来帮助
Miller et al, Particle-Based Fluid Simulation for Interactive Applications while the navier-Stokes equation 5 formulates conservation Since particles only carry the three quantities mass, posi- of momentum tion and velocity, the pressure at particle locations has to be luated first. This is de tcps. Eqn. 3)yiclds the dt Vv)=-Vp+pg+uVv density at the location of the particle. Then, the pressure can be computed via the ideal gas state equation where g is an external force density field and u the viscosity of the fluid. Many forms of the Navier-Stokes equation ap p=kp (11) pear in the literature. Eqn. (7)represents a simplified version for incompressible fluids where k is a gas constant that depends on the temperature In our simulations we use a modified version of Eqn. (11) The use of particles instead of a stationary grid simplifies suggested by Desbrurr2 these two equations substantially. First, because the number of particles is constant and each particle has a constant mass p=k(p-p0), (12) mass conservation is guaranteed and qn(6) can be omit where po is the rest density. Since pressure forces depend on ted completely. Second, the expression dv/at + v. Vv on the the gradient of the pressure field, the offset mathematically left hand side of Eqn. (7)can be replaced by the substantial has not cffcct on prcssurc forccs. Howcvcr, thc offsct docs derivative Dv/ Dt. Since the particles move with the fluid influence the gradient of a field smoothed by SPh and makes the substantial derivative of the velocity field is simply the the simulation numerically more stable time derivative of the velocity of the particles meaning that the convective term v. Vy is not needed for particle systems 3.2. Viscosity There are three force density fields left on the right hand side of Eqn. (7)modeling pressure(Vp), external forces Application of the SPh rule to the viscosity term ux2xagain (pg) and viscosity (uV-v). The sum of these force density yiclds asymmetric forccs ficlds f=-Vp+pg +uv-v determines thc change of mo- mentum p of the particles on the left hand side. For the y2W(r-r/,h).(13) acceleration of particle i we, thus, get because the velocity field varies from particle to particle. ai (8) Since viscosity forces are only dependent on velocity diffe d not on absolute velocities there is a natural where v; is the velocity of particle i and f; and p; are the force to symmetrize the viscosity forces by using velocity differ density field and the density field evaluated at the location of cncs particle i, repectively. We will now describe how we model the force density terms using SPh ISC oLly 4∑m (r;-r;h)(14) p 3. 1. Pressure A possible interpretation ofEqn. (14)is to look at the neigh bors of particle i from is own moving frame of referenc Application of the SPHrule described in Eqn. (1)to the pres Then particlc i is accclcratcd in the direction of the rclativc sure term-Vp yields speed of its environment ressure vp(r)=-∑ r今,PivW(r:-r,b) 9) 3.3. Surface Tension Unfortunately, this force is not symmetric as can be seen when only two particlcs interact. Since the gradicnt of thc We model surface tension forces(not present in Eqn. (7) kernel is zero at its center, particle i only uses the pressure of explicitly based on ideas of Morris. Molecules in a fuid particle j to compute its pressure force and vice versa. Be- arc subjcct to attractivc forccs from neighboring molccules cause the pressures at the locations of the two particles are Inside the fluid these intermolecular forces are equal in all not equal in general, the pressure forces will not be symmet directions and balance each other In contrast the forces act ric Diffcrent ways of symmetrization of Eqn. (9)havc bccn ing on molecules at the free surface are unbalanced. The net proposed in the literature. We suggest a very simple solution forces (i.e. surface tension forces)act in the direction of the which we found to be best suited for our purposes of speed surface normal towards the fluid. they also tend to mini- and stability mize the curvature of the surface The larger the curvature. the higher the force Surface tension also depends on a ten- -∑m/+P 3;VW(r;-,b (10) sion coefficient o which depends on the two fluids that form The so computcd prcssurc forcc is symmetric bccausc it uscs The surfacc of thc fluid can bc found by using an addi- the arithmetic mean of the pressures of interacting particles. tional field quantity which is I at particle locations and 0 ()The Eurographics Associa ion 2003 Muller et al, Particle-Based Fluid Simulation for interactive Applications every where else. This field is called color hield in the litera- cure For the smoothed color field we get r)=∑m/W(r-r,b The gradient field of the smoothed color field (16) Figure 2: The three smoothing kernels Wpoly6. Spiky and yields the surface normal field pointing into the fiuid and the Wwiscosiy (from left to right) we use in our simulations.The divergence of n measures the curvature of the surface thick lines show the kernels, the thin lines their gradients the directie ds the center and the dashed lines the K (17) Laplacian Nole that the diagrams ure di/ferenily scaled.The curves show 3-d kernels along one axis through the cente The minus is necessary to get positive curvature for con- for smoothing length h=I ex fluid volumes. Putting it all together, we get for the sur face traction urface oK (18) To distribute the surfacc traction among particles ncar the that is can be evaluated without computing square roots in surface and to get a force density we multiply by a normal ized scalar field 8s=n which is non-zero only near the distance computations. However, if this kernel is used for the surface. For the force density acting near the surface we get computation of the pressure forces, particles tend to build igh pressure. As particles get ver =OK oF.n (19) each other, the repulsion force vanishes because the gradient of the kernel approaches zero at the center. Desbrun solves Evaluating ni n at locations where n is small causes nu- this problem by using a spiky kernel with a non vanishing merical problcms. We only cvaluatc thc forcc if n cxcccds gradicnt ncar the centcr. For prcssurc computations we usc certain threshold Debrun's spiky kernel 3. 4. External forces 15J(h 0<r<h (r,h) 兀610 (21) Our simulator supports cxtcrnal forccs such as gravity, col- lision forces and forces caused by user interaction. These that generates the necessary repulsion forces. At the bound forces are applied directly to the particles without the use of ary where it vanishes it also has zero first and second deriva- SPH. When particles collide with solid objects such as the tIcs glass in our examples, we simply push them out of the ob- jcct and reflect thc vclocity componcnt that is perpendicular Viscosity is a phenomenon that is caused by friction and to the object's surface hus, decreases the fluid's kinetic energy by converting it into heat. Therefore, viscosity should only have a smoothing ef- 3.5. Smoothing Kernels fect on the velocity field. However, if a standard kernel is used for viscosity, the resulting viscosity forces do not al Stability, accuracy and speed of the SPh method highly de- ways have this property. For two particles that get close to pend on the choice of the smoothing kernels. The kernels each other, the laplacian of the smoothed velocity field(on we use have second order interpolation errors because they which viscosity forces depend) can get negative resulting in are all even and normalized(see Fig. 2). In addition, kernels forces that increase their relative velocity. the artifact al that are zero with vanishing derivatives at the boundary are pears in coarsely sampled velocity fields. In real-time appli- conducive to stability. Apart from those constraints, one is cations where the number of particles is relatively low, this free to design kernels for special purposes. We designed the effect can cause stability problems. For the computation of following kernel viscosity forces we, thus designed a third kernel Wholy(r, h) 315「(h2 0<r<h 0 otherwise Viscosity(r, h=-1: 一2++2-10<r<h 2h310 otherwise and usc it in all but two cascs. An important fcaturc of this simple kernel is that r only appears squared which means whose Laplacian is positive everywhere with the following Miller et al, Particle-Based Fiuid Simulation for Interactive Applications additional properties We start searches from all the cells that contain surface par- 45 ticles and from there recursively traverse the grid along the VW(r,h)=-6(h-r) 兀h0 surfacc. With the use of a hash tablc wc makc surc that thc W(r=h,h)=0 ells are not visited more than once for each cell identified to contain the surface, the triangles are generated via a fast VW(r=h, h)=0 table lookup The use of this kernel for viscosity computations increased the stability of thc simulation significantly allowing to omit 3. Implementation any kind of additional damping h, a common way to reduce the computational complexi . t Since the smoothing kernels used in SPH have finite suppe 3. 6. Simulation to usc a grid of cclls of sizc h. Then potentially interacting For the integration of the Eqn(8)we use the Leap-Frog partners of a particle i only need to be searched in is own schemel6. As a second order scheme for which the forces cell and all the neighboring cells. This technique reduces the nced to bc cvaluation only oncc, it bcst fits our purposes and time complexity of the force computation step from o(n in our examples allows time steps up to 10 milliseconds. For to o(nm), m bcing the avcragc numbcr of particles pcr grid the examples we used constant time steps. We expect even better performance if adaptive time steps are used based on With a simple additional trick we were able to speed the courant-Friedrichs-Lewy condition2 up the simulation by an additional factor of 10. Instead of storing references to particles in the grid, we store copies 4. Surface Tracking and visualization of the particlc objects in thc grid cclls(doubling mcmor consumption ) The reason for the speed up is the prox The color field cs and its gradient field n= VCs defined in imity in memory of the information needed for interpola section 3.3 can be used to identify surface particles and to lion which dramatically increases the cash hit rate. Further compute surface normals. We identify a particle i as a surface speedup might be possible through even better clustering us particle if ing Hilbcrt spacc filling curves". The data structurc for fast n(r)|>l, neighbor searches is also used for surface tracking and ren dering where l is a threshold parameter. The direction of the surface normal at the location of particle i is given by 6。 Results n(r;) (24) The water in the glass shown in Fig 3 is sampled with 2200 particles. An external rotational force field causes the fuid 4.1. Point Splatting to swirl. The first imagc (a) shows thc individual particles We now have a set of points with normals but without For the second image(b), point splatting was used to ren connectivity information. This is cxactly thc typc of infor- der the free surface only. In both modes, the animation runs mation needed for point splatting techniques24. However, at 20 frames per second on a 1. 8 GHz Pentium IV PC with rk with point clouds ob a giForce 4 graphics card. The most convincing results are tained from scanners that typically contain at least 10, 000 produccd whcn thc iso surfacc of thc color ficld is visual to 100,000 points. We only use a few thousand particles a ized using the marching cubes algorithm as in image(c) fraction of which arc idcntificd as bcing on the surfacc. Still However, in this mode the frame rate drops to 5 frames per surface splatting yields plausible results as shown in the re- second. Still this frame rate is significantly higher than the sults section one of most off-line fuid simulation techniques and with the next generation of graphics hardware, real-time performance We are currently working on ways to upsample the sur- will be possible face of the fluid. Ilereby, the color field information of sur- facc particlcs is interpolated to find locations for additional The image sequence shown in Fig 4 demonstrates inter- particles on the surface only used for rendering action with the fluid. Through mouse motion, the user gen- crates an cxtcrnal forcc ficld that causc the watcr to splash The free surface is rendered using point splatting while 4.2. Marching Cubes lated particles are drawn as single droplets. The simulation Another way to visualize the free surface is by rendering an with 1300 particles runs at 25 frames per second iso surface of the color field cs. We use the marching cubes For the animation shown in Fig. 5 we used 3000 particle algorithm" to triangulate the iso surfacc. In a grid fixcd in and rendered thc surfacc with the marching cubes tcchniquc space the cells that contain the surface are first identified at 5 frames per second ()The Eurographics Associa ion 2003 Miiller et al/ Particle-Based Fluid Simulation for interactive Applications 7. Conclusions and Future work 8. William E. Lorensen and Harvey E. Cline. Marching cubes high resolution 3d surface construction algorithm. In pro- We have presented a particle-based method for interactive eedings of the 14th annual conference on Computer graphics fuid simulation and rendering. The physical model is based and interactive techniques, pages 163 169. ACM Press, 198 on Smoothed Particlc Hydrodynamics and uscs spccial pur pose kernels to increase stability and speed. We have pre 9. L. B. Lucy. A numerical approach to the testing of the fission sented techniques to track and render the free surface of flu- hypothesis. The Astronomical ournal, 82: 1013-1024, 1977 ids. The results are not as photorealistic yet as animations 10. J.J. Monaghan. Smoothed particlc hydrodynamics. Annual computed off-line. I lowever, given that the simulation runs Review of Astronomy and A.sirophysics, 30: 543-574, 1992 at interactive rates instead of taking minutes or hours per 11. B Moon, H.V. Jagadish, C Faloutsos, and J H Saltz Analy frame as in today's off-line methods, the results are quite sis of the clustering properties of hilbert space-filling curve promising IEEE Transactions on Knowledge and Data Engineering, 13(1):124-141,2001 While we are quite content with the physical model, track- ing and rendering of the fluid surface in real time certainly 12. J. P. Morris. Simulating surface tension with smoothed particle remains an opcn rescarch problcm. In the future we will il hydrodynamics. International Journal for Numerical Methods vestigate upsampling techniques as well as ways to increase Fluids,33(3):333-353,2000 the performance of the marching cubes-based algorithm 13. S.A. Munzel Smoothed particle hydrodynamics und ihre An wendung auf Akkretionsscheiben. Phd thesis, Eberhard-Karls Universitat Thiibingen, 1996 Acknowledgements 14. D. Nixon and R. Lobb. A Huid-bascd soft-objcct modcl The authors would like to thank Simone Hieber for her help- aphics and Applications, pages ful comments and Rolf Bruderer. Simon Schirm and Thomas July/ august 200 Rusterholz. for their contributions to the real-time system 15. D. Pnueli and c. gutfinger. luid Mechanics. Cambridge Univ Press. NY.1992 16. C. Pozrikidis. Numerical Computation in Science and Engi Refe neering. Oxford Univ Press, NY, 1998 1. Mark Carlson. Peter Mucha III R. Brooks van Horn, and 17. W. T Reeves. Particle systems-a technique for modeling a Greg Turk. Melting and nowing. In Proceedings of the ACM class of fuzzy objects. ACM Transactions on Graphics 2(2) SIGGRAPH symposium on Computer animation, pages 167 pages91-108,1983 174. ACM Press. 2002 18. Jos Stam. Stable fluids. In Proceedings of the 26th annual 2. M. Desbrun and M. P. Cani. Smoothed particles: Anew on/erence on Computer graphics and interactive techniques paradign for animating highly deformable bodies. Ir pages 121-128. ACM Press/Addison-Wesley Publishing Co Computer Animation and Simulation 96(Proceedings of 1999 EG Workshop on Animation and Simulation), pages 61-76 Springer-Verlag, Aug 1996 19. Jos Stam and Eugene Fiume. Depicting fire and other gaseous phcnomcna using diffusion proccsscs. Computer Graphics, 3. D Enright, S Marschner, and R. Fedkiw. Animation and ren 29(Annual Conference Series): 129-136, 1995 dering of complex water surfaces. In Proceedings of the 29th annual conference on Computer graphics and interactive tech- 20. Dan Stora, Pierre-Olivier Agliati, Marie-Paule Cani, Fabrice niques, pages 736-744. ACM Press, 2002 Neyret, and Jean-Dominique Gascuel. Animating lava flows In Graphics Interface, pages 203-210, 1999 4. N. Foster and R. Fedkiw. Practical animation of liquid In Proceedings of the 28th annual conference on Computer 21. T. Takahashi, U. Hcihachi, A. Kunimatsu, and H. Fujii. The graphics and interactive techniques, pages 23-30. ACM Press simulatiun uf nuid-rigid body interaction. ACM Siggraph 2001 Sketches Applications, July 2002. 5. R. A. Gingold and JJ. Monaghan. Smoothed particle hy 2. D. Tonnesen Dynamically Coupled Particle Systems for Gen drodynamics: theory and application to non-spherical stars metric Modeling, Reconstruction, and Animation. PhD thesis University of Toronto. November 1998 Monthly Notices of the royal Astronomical Society, 181: 375- 398,1977 23. Andrew Witkin and Paul Heckbert. Using particles to sample and control implicit surfaces. In Computer Graphics(Proc 6. Damicn Hinsingcr, Fabricc Ncyrct, and Maric-Paulc Cani In- SIGGRAPH 94), volume 28. 1994 teractive animation of ocean waves. In Proceedings of the ACM SIGGRAPH symposium on Computer animation, pages 24. Matthias Zwicker, Hanspeter Pfister, Jeroen van Baar, and 161-166. ACM Press.2002. Markus gross. Surface splatting. In Proceedings of the 28th annual conference on Computer graphics and interactive tech 7. Jcan-Christophe Lombardo and claude Pucch Oricnted parti cles: A tool for shape memory objects modelling In Graphics niques, pages 371-378. ACM Press, 2001 Interface 25, pages 255 262, mai 1995. Quebec city, Canada Miller et al/ Particie-Based Fiuid Similation for Ineractive Applications (b) splatting and (c) the iso-surface triangulated via marching cuber mage (a) shows the particles, (b) the surface using point Figure 3: A swirl in a glass induced by a rotational force field Figure 4: The user interacts with the fluid causing it to splash Figure 5: Pouring water into a glass at 5 frames per second. ()The Eurographics Associa ion 2003