Modeling Highly-Deformable Lluid

全文

(1)Modeling Highly-Deformable Lluid Chih-Wei Chiu Jung-Hong Chuang Cheng-Chung Lin Department of Computer Science and Information Engineering National Chiao Tung University Submitted to ICS Workshop on Multimedia Technologies, 2002. Contact:. Jung-Hong Chuang. phone: fax: email:. 886-3-5731829 886-3-5724176 jhchuang@csie.nctu.edu.tw Abstract. Modeling liquid with complex surfaces is a great challenge in computer graphics. We proposed a hybrid approach combining the volume tracking method and smoothed particle hydrodynamics (SPH) to model liquids with highly deformable surfaces. The 3-D NavierStokes equations are solved by the marker-and-cell (MAC) method and a density volume representing the liquid is evolved over time and space by the volume-of-fluid (VOF) method. The VOF method uses a finite-volume discretization and maintaining a value for the liquid volume in each grid cell. The advection of the volumes is calculated by using cell-face fluxes, in which the liquid volume that leaves one cell is exactly the same as the liquid volume that enters the adjacent cells. Only a scalar convective equation needs to be solved to evolve the liquid surface forward in time and space. Potentially under-resolved effects are captured by SPH particles generated in that high-deformable regions near the liquid surface. Dynamics of these particles provides a better fieldity than that of other particle methods due to their capability of approximating the equations of fluid mechanics. Particle volumes are incorporated into the interpolated volume fractions before an unified iso-surface is extracted. As a result the liquid surface is more smooth than those of previous researches which rendered the splash by coating the particles with a field functions or as hard spheres.. Keywords: Liquid simulation, Fluid dynamics, Splashing, Particle dynamics..

(2) 1. Introduction. Producing realistic liquid animation is a great challenge in compute graphics due to the complex dynamics of flows. Various approaches have been attempted. Early graphics work concentrated on modeling just the surface of a body of water as a parametric function that could be animated over time to simulate wave transport. Recently, many researches have simulated various natural objects and phenomena by using a set of particles each of which has simple behavioural rules. Construction of the behavioural model is based on empirical and intuitive knowledge wihtout taking a standard apparoch to numerical simulation. While a large volume could potentially be modeled with particles, the number of particles needed to fill a volume will grow with the cube of the resolution of the scene. Moreover, despite the simplified equations for the dynamics this model can be computationally expensive for a large volume of particles. On the other end of the spectrum lies the Eulerian volume tracking method, which has already enjoyed great success in the computational fluid dynamics, but not yet widely used in the computer graphics community. An Eulerian method is characterized by a stationary grid partitioning the whole computational domain, whileas using particles is a Lagrangian method characterizing by a “moving” grid, which uses the particles to be the grid nodes. In the volume tracking method, the scene is partitioned by a Cartesian grid (regular or irregular). The initial fluid configuration is converted to a scalar field in the grid. A scalar value between zero and unity, known as the “volume fraction,” indicates the proportion of liquid and air in the cell. Exact interface information is then discarded in favor of the discrete volume fraction data. Interfaces are subsequently tracked by evolving fluid volumes in time with the solution of a standard convection equation. The major negative aspect of particle methods is the cost in terms of CPU time and storage. To update a particle’s position, the intergration of the sum of the forces acting on this particle by others needs be to calculated. Using volume fractions is more economical than particles in both space and time, as only one value (the volume fraction) needs to be stored for each grid cell. Only a scalar convective equation needs to be solved to propagate the volume fractions through the 2.

(3) computational domain. Another drawback of particles method is that there is no straightforward way to extract a smooth surface from the particles, since deducing the connectivity of particles is difficult or impossible. The volume data results from the volume tracking method can be easily visualized by employing an isosurface extraction algorithm. In complex flows, a liquid surface can be stretched and torn in a dynamics fashion, the use of only an initial seeding of particles will not capture these effects well, since regions will form that lack a sufficient number of particles to represent the fluid bodies. This would require adding extra particles when the surface becomes too sparsely resolved, and removing them as the surface folds, or splashes back over itself. Another advantage of an Eulerian volume tracking method over a Lagrangian particles method is the ability to withstand severe topological change robustly, since no topological or connectivity information is retained and the arbitrary complex interface can be inferred from the volume fractions. This enables the modeling of interfaces which are exposed to large deformations. The merging or breakup of the interface is handled in a natural way. Vorticity driven stretching, tearing and splitting near the surface typifies the topological changes that must be captured. These surface features may be under-resolved for volume tracking methods with a coarse computational grid. While volume tracking is a fully Eulerian approach, the ability of volume tracking is restricted by the resolution computational grid. Particles evolution is a fully Lagrangian approach that has strength of resolving very fine effects, such as small droplets and water splash. Chief among the strengths of particles method is the ability to continue to advect bodies faithfully even when the body can not be supported on a grid. Since they tend to have complementary strengths and weakness, a combined approach gives superior results under a wider variety of situations [7, 9, 25, 27]. We propose a method incorporating the volume tracking and smoothed particle dynamics to model complex fluid motion. Fluid surface is tracked by volume fractions and evolves over space and time. Particles are only generated in the area where the fluid motion is drastic, for example, when an object impacts the fluid surface. These particles move freely until they run into the volume fractions and are absorbed. The motions of particles are governed by smoothed particle hydrodynamics (SPH) [21], which are used by physi3.

(4) cists for cosmological fluid simulation. Smoothed particles represent sample points that enable the approximation of the values and derivatives of local physical quantities inside a medium. SPH is favored than particle systems [29], which do not model inter-particle forces, and molecular dynamics [12], which does not model the coupling of pressure and velocity characterizing a real fluid. The paper is organized as follows. In Section 2, we review the previous work on modeling fluid from three aspects: the dynamics model, the geometry model, and the illumination model. Our method is outlined in Section 3, with detailed descriptions about a numerical method to solve the 3-D Navier stokes equations; the surface kinematic condition; the necessary modification of SPH to model incompressible flow; and the rendering of our geometry model. Section 4 gives the result and discussion. A conclusion and future work are presented in Section 5.. 2 2.1. Fundamentals Eulerian and Lagrangian Reference Frame. There are two equivalent ways of describing a fluid when it is in motion. In the Eulerian frame we describe some aspects of the fluid in terms of what takes place at a fixed point as the fluid flows by. We ae thus continually observing new particles as they pass the point of observation. In the Lagrangian frame we pick out a single particle of the fluid and follow it along as it partakes of the fluid motion. We can then describe the same aspect of the fluid that is of interest before but now in terms of what happens to the fluid in the immediate vicinity of a specific moving particle.. 2.2. Navier-Stokes Equations. In general, we would be interested in fluid properties that arise in the expression of the laws of conservation of mass, momentum, energy, etc. The transportation equation for the conservation of mass is derived from Reynolds’ Transport Theorem as follows. Dρ + ρ∇ · v = 0. Dt 4. (1).

(5) This relationship, also known as the continuity equation, states that the mass of a fluid in a close domain can only be changed by flow across the boundaries. The following transport equation for the conservation of momentum is also derived from Reynolds’ Transport Theorem. 1 2 Dv 1 + ∇p − ∇ v − g = 0. Dt ρ Re. (2). The only external force acting on fluid volume considered here is g, the force due to gravity. (1/Re)∇2 v is the viscous drag. This is the key to producing the realistic flow of water–the component of motion leading directly to self-propagation rotation. Rotation motion is the most important visual part of a fluid’s motion and must be hte ingredient of any realistic graphics model. Re is called the Reynolds number Re ≡. V ∗ L∗ ρV∗ L∗ = , µ ν. where V ∗ and L∗ are the characteristic velocity and length, respectively. Re may be interpreted as the magnitude of the ratio of of inertial to viscous force in a flow. While not obvious, pressure p plays an important role in the behavior of the surface of a volume of water. Pressure waves internal to the volume (acoustic waves) can travel much faster than the surface waves we are all familiar with. In (2), the velocity and pressure of fluid element are highly coupled, so they can characterize the effects that are just not possible using surface-based or particle-force techniques. Thus, much of what we see at the surface has been influenced by subtle (and sometimes not-so-subtile) changes in pressure within the volume.. 2.3. Smoothed Particle Hydrodynamics. One reason that the use of smoothed particle hydrodynamics (SPH) [21] over other numerical methods such as finite difference method (FDM) is that SPH does not depend on the boundary conditions. The FDM is based on breaking up the computational region into a grid. The reason for creating a grid is that in order to calculate the spatial derivative of the field variables, one must divide the difference of the filed variables at two points by the distance between the two points. In SPH, the fluid is sampled by a set of elements called particles (which may also be regarded as 5.

(6) interpolation points). A particle j has a fixed mass mj , a position rj , a velocity vj , and a density ρj depending on the local density of particles. As a sample point, it can also carry physical field values like pressure or temperature. Then in a way very similar to Monte-Carlo techniques, these fields and their derivatives can be approximated by a discrete sum. At the heart of SPH is an interpolation method which allows any function to be expressed in terms of its values at a set of disordered points. The integral interpolant of any function ϕ(r) is defined by ϕ(r) =. Ω. ϕ(r )Wh (r − r , h)dr ,. and Wh is an interpolating kernel which has the two properties Ω. Wh (r − r , h)dr = 1,. and lim Wh (r − r , h) = δ(r − r ),. h→0. where the limit is to be interpreted as the limit of the corresponding integral interpolants. This normalized kernel gives the spatial mass distribution profile over a smoothing length h. Then the smoothed values and derivatives of a continuous field ϕ known only at particle locations can be approximated by . ϕj Wh (r − rj ), ρj j ϕj < ∇r ϕ(r) > = mj ∇Wh (r − rj ). ρj j < ϕ(r) > =. mj. (3) (4). The essential point is that we can construct a differentiable interpolant of a function from its values at particles (interpolation points) by using a kernel which is differentiable. Derivatives of this interpolant can be obtained by ordinary differentiation; there is no need to use finite difference and no need for a computational grid. This is one of the main reasons SPH is so popular. It removes the need for a mesh to calculate spatial derivatives. With equation (3) and (4), we are ready to write down our numerical equations for the smoothed values of each point. From here on, we will remove 6.

(7) the brackets and place a subscript i on the value. This implies that the value being calculated is the smoothed value for particle i. The rewritten version of our two SPH equations are then: . ϕj Wij , ρj j ϕj < ∇r ϕ(r) > ≡ ∇i ϕi = mj ∇i Wij , ρj j < ϕ(r) > ≡ ϕi =. mj. where Wij = Wh (ri − rj , h) and ∇i implies the spatial derivative with respect to particle i’s coordinates. For example the density at particle i may be evaluated using ρi =. . mj Wij .. (5). j=i. Writing the continuity (1) in the form Dρ = −∇ · (ρv) + v · ∇ρ, Dt and using SPH particle interpolants for the RHS, the rate of change of the density of particle i becomes. dρi = mj vij · ∇ij Wij , dt j=i. where vij = vi − vj . The momentum equation for particle i becomes pi dvi pj =− mj + 2 + Πij ∇ij Wij + g, 2 dt ρ ρj i j=i. (6). where Πij is the viscous term given by  2   −αcµij +βµij , vij · rij < 0; ρij Πij =   0, otherwise and µij =. hvij · rij . rij + (0.01h)2. In these expressions the notation ρij = (ρi + ρj )/2 and rij = ri − rj has been used, and c is the speed of sound. Because of its symmetry the viscous term conserves linear and angular momentum. The first term involving α is analogous to a shear and bulk viscosity. The second one, comparable to the von Neumann-Richtmyer artificial viscosity used in grid-based method, prevents particle interpenetration at high speed. 7.

(8) 3. Previous Work. Any computer graphics techniques for visualizing a fluid has to address the issues of [1]: • capturing the dynamic behavior, creating a dynamics model ; • representation of complex fluid geometry, creating a geometry model ; and • modeling the light interaction with the fluid, creating an illumination model These three models are tightly coupled. A dynamics model may restrict the possible geometry models on representing the fluid, and the illumination model may have to be adapted to a specific geometry model.. 3.1 3.1.1. Dynamics Models Particle Methods. Goss used a particle system to model ship wakes in real-time, but this approach was not based on physics [13]. O’Brien et.al. used a spray particle system to model water drops which have broken free of the main body of fluid [27]. The dynamics of the spary subsystem are quite simple since the particles are not interacting with one another; thus the model is simply a collection of particles moving under gravity. The model can deal with an impact problem in a natural way. The impact of the object on the free surface is translated to a pressure on the affected columns. Particles are created when the upward velocity of a portion of the surface exceeds an empirically predefined threshold. This threshold is defined heuristicly to be adjusted to appearance of the splash. Low threshold value will cause more particles to be generated over a larger area. The volume of each particle is subtracted from the column from which it was created, and is removed from the system when it falls back onto the surface, and the volume of each particle is added to the column that absorbed it. One of the main problems with the splash model is its symmetric. An extended particle system based on O’Brien et.al. [27] is proposed in [25]. The condition for droplets created from the surface are now twofold: first, the absolute vertical velocity must exceeds 8.

(9) a threshold, as in [27]; second, the vertical velocity adjusted for water depth must exceed another threshold. A droplet may split into two smaller droplets, these descendent droplet may further split into even smaller droplets. The splitting mechanism depends upon surface energy. There are two values relevant to determining whether a droplet should split: the free energy, the energy contained in the droplet which can be devoted to splitting it, and the threshold energy, the energy requirement for dividing. To avoid identifiable pattern and add realism, the model incorporated some stochastic elements in order to properly capture the random look of splashing water. For example, The split conditions are tested stochasticly and the radii of two resultant droplets are determined by trial-and-error. A very economic model for splashing in thin films of water is given in [2] with the thin film constraint. It assume a stationary fluid film and that ripples are created only during an elliptical foot impact. The spray model has been segmented into the main splash and the escaping free droplets. When the foot lands, the main splash (termed the primary funnel) is simulated by a torus whose base radius is determined by the foot radius. The secondary radius of the torus grows and decays under the influence of gravity, and the crest curve angle also grows with time, to simulate the falling over of the funnel crest. Drops are generated around the circumference of the primary funnel crest when the funnel is in the rising phase. The space between droplets is stochastically distributed and is also a function of the region of the funnel with respect to the foot. Globular dynamics, or molecular dynamics, is a connected particle system with the uses a dynamics model for particle to particle interaction and interaction with other environmental constraints. Each particle has an associated radial force field which leads to a dynamically changing topology of interactions. We will refer to an element of the connected particles systems as a globule, thus avoiding the established connotations associated with words such as particle or blob. While a large volume could potentially be modeled with globules, the number of globules needed to fill a volume will grow with the cube of the resolution of the model, and despite the simplified equations for the dynamics this model can be computationally expensive for a large volume of particles. SPH is modified to simulating high deformable bodies in [5]. Stora et.al. extended SPH to include 9.

(10) a viscosity parameter depending upon temperature for each particle [32].. 3.1.2. 3-D Navier-Stokes Equations. Previous methods do not model the coupling between velocity and pressure in a fluid. It’s this coupling that is responsible for most of the behaviors we associate with fluids. The lack of velocity and pressure information also makes accurate control of the surface very difficult. There is not enough information upon which to apply the equation of motion, so control is only possible using ad-hoc method, such as raising a column of the height field to generate a wave. Without the pressure it’s also difficult to include buoyant objects in a scene because the information necessary to calculate their motions is unavailable. To provide a physical foundation for general fluid animation, one must employ the 3-D NavierStokes equations (1) and (2), which are the embodiment of Newton’s second law in fluids, and the governing equations of general fluid flow. The simulation of complex water effects using the full 3-D Navier-Stokes equations has been based upon the large amount of research done by the computational fluid dynamics (CFD) community over the past 50 years. Foster et.al. utilized the work of [15] in developing a 3-D Navier-Stokes methodology for the realistic animation of liquids [10]. Velocity and pressure are coupled to achieve more realistic three dimensional behavior compared to previous models, be calculating momentum throughout the whole volume. This does allow for subtle fluid effects to manifest at the surface, but does not fully exploit the fact that the calculation is made for the entire environment. Also, although the motion of buoyant objects within the flow could be calculated, the objects could not affect the flow in any way. The similar technique was also used to model hot gases [11]. A semi-Lagrangian“stable fluids” treatment of the convection portion of the Navier-Stokes equations was introduced to the computer graphics community by [31] in order to allow the use of significantly larger time steps without hindering the stability. Such a semi-Lagrangian advection has been used widely in recent researches [9, 7, 8, 19]. Yngve et.al. used the finite difference method to solve the compressible Navier-Stokes equations to model shock wave and convection effects generated by an explosion [33]. 10.

(11) 3.2. 3.2.1. Geometry Models. Implicit Surfaces. The problem of free surface of fluids was treated as a transport problem of a density function (volume fractions) in [18], which describes the ratio of the liquid in each voxel in the fluid simulation space. To visualize these values, the high-resolution Marching Cubes algorithm [20] was employed to extract the iso-surface. An inherent problem with many iso-surface extraction algorithm, including the Marching Cubes, is the inability to express density values smaller than the iso-value. The authors suggested a future work to display small density values as particles or semitransparent polygons based on some rules.. Foster et.al. made significant contributions to the simulation and control of three dimensional fluid through the introduction of a hybrid liquid volume model that combines an implicit surface generated from the level set method (LSM) and massless marker particles; the formulation of plausible boundary conditions for moving objects in a liquid; a time step sub-cycling scheme for the particle and implicit surface evolution equations in order to reduce the amount of visual error inherent to the large semi-Lagrangian stable fluid time-step [9].. Inspired by the hybrid liquid volume model proposed in [9], Enright et.al. introduced a ”thickened” front tracking technique, called the ”particle level set method [6],” a hybrid surface tracking method that uses massless marker particles combined with a dynamic implicit surface, to accurately represent the water surface, exhibiting significantly more realistic surface behavior [7]. This effect is achieved by focusing on modeling the surface as opposed to the liquid volume as was done in [9]. This shift in philosophy away from volume modeling and toward surface modeling is the key idea resulting in better photo-realistic behavior of the water surface. However, this method does not preserve the volume of fluid. 11.

(12) 4. Modeling Highly-Deformable Liquid. 4.1. Overview. Figure 1 illustrates the computation flow in a timestep. Initially the physical domain is partitioned by a computation grid, and this partition is only performed once before the simulation begins. The fluid dynamic solver then calculates the velocities and pressure for the next time step. The volume fractions in the current time step and new obtained velocities are used to convect the fluid volume fractions. In the volume tracking method, the accuracy of fluid interface depends upon the resolution of the computational grid. Some surface detailed features may be under-resolved if the grid is too coarse. Particles are introduced in the area where the fluid surface is potentially exposed to large deformations. Dynamics of particles is governed by smoothed particle hydrodynamics (SPH). Partition the scene by a computional grid. pn and v n. grid. Fluid dynamics solver. v n+1. fn. Evolve volume fractions. f n+1 Generate particles in the drastic motion regions existing particles. new particles. Move particles by SPH. Figure 1: Simulation flowchart.. 12.

(13) 4.2. Solving Fluid Dynamics. The restriction to incompressible flow introduces the computational difficulty that (??) contains only velocity components, and there is no explicit link with the pressure as there is for compressible flow through density. One of the earliest and most widely used methods for solving (??) and (??) is the marker and cell (MAC) method [15]. The method is characterized by the use of a staggered grid and the solution of a Poisson equation for the pressure at every time step. Although the original form of the MAC method has certain weakness, the use of a staggered grid and a Poisson equation for the pressure has been refined in many modern methods derived from MAC.. To track the free surface, the initial (known) fluid interface geometry is used to compute fluid volume fractions in each computational cell. A scalar indicator function between zero and one, known as the volume fraction, f , is used to distinguish between two different fluids, liquid and air. A value of zero indicates the absence of liquid and a value of unity corresponds to a cell full of liquid. On a computational grid, volume fraction values between these two limits indicates the presence of the interface and the value itself gives an indication of the relative proportion of liquid and air occupying the cell volume. Exact interface information is then discarded in favor of the discrete volume fraction data. Interfaces are subsequently tracked by evolving fluid volumes in time with the solution of a standard convection equation.. Special care must be taken to approximate the convection equation to keep the interface compact, since using even the most well designed finite difference methods to advect interface results in an interface thickness of at least two to three cells in width. We employ the improved donor-acceptor scheme [16] that satisfies local boundedness while keeping the transitional area between the liquid and air restricted to one cell width. The consistent discretization of the continuity and momentum equations with respect to the discretization of the volume fraction equation is used to ensure a solution algorithm which satisfies the conservation of the flow properties at all times. 13.

(14) 4.3. Modeling Splash. To compensate the motion under-resolved by the resolution of the Eulerian grid, we place Lagrangian particles near the region in the flow field having drastic motion. Section 4.3.1 describes how these regions are identified and how particles are positioned inside the cell’s “cut-volume”, defined by the volume cut by a 3-D plane approximating the liquid interface in the cell. Once generated, the particles may escape the liquid surface and their dynamics are governed by SPH. SPH is favored than ordinary particle systems since the former can approximate the NavierStokes equations and enables particles to interact with one another in the range of the kernel, thus providing a more realistic simulation. Section 4.3.2 describes some extensions that are necessary for SPH to be used in the computer grapihcs.. 4.3.1. Particle Seeding. An obvious way to place particles is to put them in the surface cells, whose volume fractions are great than zero but less than unity. These cells are where the grid potentially under-resolve the interface. But not all surface cells have violent changes, particles in these cells do not resolve the interface much better than the volume fractions. To save computation time, particles are generated only in the surface cells meeting any one of following conditions: uδt > γδx,. vδt > γδy,. wδt > γδz.. (7). where δx, δy, and δz are cell sizes, δt is the time step, u, u, w are velocity components, and γ ∈ (0, 1) is the convection threshold. Inequality (7) can be related to the CFL (Courant-Friedrichs-Lewy) condition, which physically states a particle of fluid should not travel more than one spatial step-size δx in one time step to keep the numerical method stable. A high γ values tends to produce less particles than a lower one. Each smoothed particle initially carries a fixed mass and will never split or merge with other particles. The total volume of particles in a cell should not excced the volume of liquid in the cell indicated by the volume fraction. The mass carried by a particle can be arbitrarily defined. The less the mass, the more particles are necessary to approximate a same volume of liquid in the cell. 14.

(15) There is another restriction on the maximum number of particles in a cell. If the number of particles in the cell have exceeded a predefined maximum, particles will not be seeded, since the existing particles should provide accurate approximation of the interface and in this way we prevent from potentially introducing excessive particles to the SPH simultion. To seed particles in the a surface cell, the interface in the cell must be reconstructed, a process similar to what we do in volume fraction convection. Particles will then be randomly placed below the interface. Initial velocities of particles are linearly interpolated from the velocities defined on the cell faces according to their positions, as shown in Figure 2. Conversely, particles entering the area below the interface are removed from the system. vi+1/2,j,k. cut-plane ui−1/2,j,k. ui+1/2,j,k. vi−1/2,j,k. Figure 2: Interpolate particle velocity from the cell velocity. When a piecewise constant/“stair-stepped” approximation is used in convection volume fractions, the interface is inadequate for seeding particles. Because reconstructed piecewise constant interface is parallel to one of the axes, particles positioned below these axis-aligned planes will form a regular pattern, which is not natural in a realistic animation. We approximate the interface by a Piecewise Linear Interface Calculation (PLIC) method [34]. Note that the PLIC method is not used for convecting volume fractions for efficiency, while a piecewise constant/stair-stepped approximation is accurate enough for simulation. In a PLIC method, the interface between materials is approximated in each computational cell with a linear interface defined by the equation n · x = λ, 15. (8).

(16) where n is the normal of the interface plane and λ is a scalar. Any cell having volume fractions f between zero and one will possess an interface defined by (8). We choose n to point outward the fluid, hence (8) will be positive for any point x lying within the fluid, zero for any point x lying on the line, and negative for any point x lying outside of the fluid. Determining λ is the most difficult reconstruction task because the value of λ is constrained by mass conservation. In other words, the value of λ is constrained such that the resulting plane passes through the cell with a truncation volume equal to the cell material volume fi,j,k . This determination requires inverting a V (λ) relation. However, the determination of n is not constrained by volume conservation consideration. We follow the approach in [34] by using finite difference stencils, yielding an explicit expression of n. The three components of the normal in the cell (i, j, k) is given by,. nx = fi+1/2,j,k − fi−1/2,j,k ny = fi,j+1/2,k − fi,j−1/2,k nz = fi,j,k+1/2 − fi,j,k−1/2. This method is shown to be first-order accurate, as measured with the least square criteria. We distinguish the “forward problem”, that is, to find the volume fraction Vf = fi,j,k occupied by one material given λ, from the “inverse problem”, which consists of finding λ given the volume fraction. The problem is essentially geometric in nature. The relation λ = V −1 (Vf ) is nonlinear and varies in each mixed cell, due to its dependence upon local data. A case-by-case implementation results that is not efficient, general, concise, or easily to maintained and understood. Therefore we employ an analytical expression connecting the volume fraction f and λ in this thesis. We assume that the three components of n are all positive and we need to determine the cut volume ABGHKNML of the rectangular cell which is also below the given plane, as shown in Figure 3. [14] have shown that the volume ABGHKNML is given by the expression  V =. 1  λ3 − 6nx ny nz. . Ψ(λ − ni δi) +. i∈{x,y,z}. i∈{x,y,z}. 16.  Ψ(λ − λmax + ni δi) ,. (9).

(17) z. C. D K. N F. E A. B. J. y. L M H. G. P. O I x. Figure 3: The “cut volume” is the region inside the parallelepiped ABCEFGH and below the plane IJK where the function Ψ(y) is defined as. Ψ(y) =.    y 3, for y > 0,   0,. otherwise,. and λmax = nx + ny + nz For the moment, we restrict our analysis to a unity cube whose δx = δy = δz = 1; then volume V and volume fraction Vf coincide. We also normalize the plane equation (8) by dividing it by (nx + ny + nz ); then λmax = 1. Note that expression for V is invariant with respect to a permutation of the indices, so we need to consider only one case, say nx ≤ ny ≤ nz . The graph of V has odd symmetry with respect to the point (V, λ) = (1/2, 1/2), so we restrict the analysis to the range 0 ≤ λ ≤ 1/2. Let nxy = nx + ny , and n = min(nxy , nz ). The range of λ, [0, 1/2], is partitioned into four 17.

(18) subintervals [0, V1), [V1 , V2 ), [V2 , V3 ) , and [V3 , 1/2]), where V1 =. n2x , max(6,ny nz ,). V2 = V1 +. ny −nx , 2nz. V3 = min(V31 , V32 ), where is an arbitrary small number and V31 =. n2z (3nxy −nz )+n2x (nx −3nz )+n2y (ny −3nz ) , 6nx ny nz. if n = nz ,. V32 =. nxy , 2nz. if n = nxy .. The inverse problem is given by [30] as follows:. 3. 6nx ny nz V , for 0 ≤ V < V1 , . λ = 12 nx + n2x + 8ny nz (V − V1 ) , for V1 ≤ V < V2 ,. λ =. P1 (λ) = a3 λ3 + a2 λ2 + a1 λ + a0 = 0,. for V2 ≤ V < V3 ,. and for V3 ≤ V ≤ 1/2, we have P2 (λ) = a3 λ3 + a2 λ2 + a1 λ + a0 = 0, V31 ≤ V32 λ = nz V +. nxy , 2. V32 ≤ V31. For coefficients of the two cubic polynomials we have a3 = −1, a2 = 3nxy , a1 = −3(n2x + n2y ), a0 = n3x + n3y − 6nx ny nz V, a3 = −2, a2 = 3, a1 = −3(n2x + n2y + n2z ), a0 = n3x + n3y + n3z − 6nx ny nz V. In previous relations V1 is an approximation of the value n2x /6ny nz . This approximation is needed because the limit for V1 , as nx , ny → 0, is zero; however numerically we must prevent the denominator of V1 from becoming zero.. 4.3.2. Particles Dynamics. After the initial positions and velocities of smoothed particles are specified, their trajectories and acceleration are governed by SPH. 18.

(19) Extension to SPH If (5) is used for incompressible fluids like water, where the density falls discontinuously to zero at the surface, the density, however, will be smoothed over the length 2h and surface particles will have a low density [22]. The equation of state will then introduce incorrect pressures and degrade the calculation. It is therefore preferable to depart from normal practice and approximate the rate of change of the density. All particles are then assigned the same initial density which only changes when particles are in relative motion. There are at least two ways that SPH might be extended to incompressible flow. The first is to work directly with the constraint of constant density. It’s possible to include these constraints easily in the SPH formulation by using the Gibbs–Appel equation [28] which are generalized versions of Gauss’ principle of least constraint. Unfortunately, the resulting equations are cumbersome, and it has not been possible to solve them efficiently without further approximation. The second approach [4] is based on the observation that real fluid such as water are compressible, but with a speed of sound which is very much greater than the speed of bulk flow. The equation of state has the form. p = p0. ρ ρ0. . γ −1. (10). with γ = 7, and a zero subscript denotes reference quantities. The choice of γ = 7 in (10) causes pressure to respond strongly to variations in density. Thus, perturbations to the density field remain small. However, as the density fluctuations increase, small errors in density correspond to increasingly larger errors in pressure [24]. In previous work involving incompressible fluids, the subtraction of 1 in (10) was introduced to remove spurious boundary effects at a free surface. It is well established that SPH is unstable when attractive forces act between particles [3]. Consequently, for the simulation presented in thesis, we use p = c2 ρ,. (11). as suggested by [24]. The speed of sound must be chosen carefully to ensure an efficient and accurate solution of a given problem. The value of c must be large enough that the behavior of the corresponding quasi19.

(20) incompressible fluid is sufficiently close to that of the real fluid, yet it should not be so large as to make the time step prohibitively small. One way to find appropriate c is to balance the force in Navier-Stokes momentum equation (2) to come up with several relations to which the the speed of sound is in proportion[24]. The speed of sound should be comparable with the largest of V∗ 2 V∗ gL∗ , , , δ ReL∗ δ δ. (12). where δ = ∆ρ/ρ0 . The first term in (12) corresponds to that derived by [22]. The second and third terms ensure that pressure forces are comparable with viscous and body forces, respectively. The expression of pressure in (11) resulted in positive, i.e., purely repulsive, forces expressing the natural expansion of the fluid. However, we would like to animate material with constant density at rest. Consequently, the material should exhibit some internal cohesion, resulting in attractionrepulsion forces as in the Lennard-Jones model. To keep density at ρ0 , we replace the equation of state by p = k(ρ − ρ0 ).. (13). (13) designed to maintain density close to a constant value, has a double advantage [5]. First, if particles have the same mass they will tend to be evenly distributed inside the object. This is essential since we are using them as sample points for approximating continuous function. Moreover, constant density results in a constant volume. The material will then tend to naturally come back to is initial volume after a compression. Replace p by its value in (6) leads to ρi − ρ0 ρi − ρ0 dvi =− + mj + Πij ∇ij Wij + g. kmj dt ρ2i ρ2j j=i. (14). The first term in parentheses of (14) is a density gradient descent, that tends to minimize the difference between current and desired densities. The second is a symmetry term that ensures the action-reaction principle. The parameter k determines the strength of the density recovery. Its plays the same role as a stiffness parameter in a standard molecular dynamics. A larger k will simulate a stiff material while a smaller k models a soft one. 20.

(21) After the acceleration of the particles is obtained from (14), the velocities and positions are updated by the Leapforg scheme: n− 21. vin = vi n+ 21. vi. n+ 21. ri. n− 21. = vi. δt dvin , 2 dt dvn + δt i , dt +. n+ 21. = rni + δtvi. .. The Kernel The use of different kernels in SPH is analogue to the use of different difference schemes in finite difference method. The advantage of SPH is that the kernel can be calculated in a subroutine, or a table, and it is then trivial to change a code with one kernel into a code with another. The kernel based on spline functions [23]    1 − 32 s2 + 34 s2 , if 0 ≤ s ≤ 1,   σ  1 W (r, h) = ν (2 − s)3 , if 1 ≤ s ≤ 2, 4 h      0, otherwise where s = r/h, ν is the number of dimensions and σ is a normalization constant with the values 2 10 1 , , 3 7π π. for one, two and three dimensions respectively, has advantages. Thos kernel has compact. support; the second derivative is continuous, and the dominant error term in the integral interpolant is O(h2 ).. 4.4. Rendering. The liquid in the simulation is represented by volume fractions and smoothed particles. Although the volume fractions can be easily rendered by extracting an iso-surface from the density volume, the inability of representing density values less than the iso-value will prevent the animation sequence from preserving a constant volume [18]. Particles, however, when serving as the point skeletons of an implicit surface, can not form a smooth surface unless the number is very large and evenly distributed. To produce the smooth surface of realistic liquid, another volume fractions of higher resolution (the “fine” density volume) are interpolated from the original volume fractions (the “coarse” density 21.

(22) Figure 4: Point skeletons can not easily form a smooth surface volume), and the volumes of smoothed particles are converted to volume fractions of the occupying cells in the fine density volume. In the cells where both volume fractions and smoothed particles exist, the converted volume fractions of smoothed particles always predominate under the assumption that particles should give accurate values in the regions with drastic fluid motion, since they tracks the fluid surface independent of the resolution of the computational grid. Smoothed particles are eliminated when they fall back into the liquid volume or outside the simulation domain. At last, an iso-surface extraction alogrithm (Marching Cubes here) is employed to extract a smooth surface. While the volume fractions are interpolated when rendering the liquid surface, they are not used in the future simulation process. Figure 5 illustrates the rendering procedure and Figure 6 demonstrates the results.. 5. Results. All images in this thesis are rendered by a ray-tracer at the resolution of 640×480 pixels. All computation and rendering are performed on a PC armed with an AMD 1000 MHz processor and 512 MB RAM. Figure 7 shows an orange liquid is discharge from a pipe at the speed of 1.8 into a container of dimension 6 × 5 × 5 and of resolution 30 × 24 × 24. Figure (a)-(c) shows the frame at roughly the 42th, 92th and 124th minute of the simulation, and takes 12, 10, and 8 minutes to render, respectively. 22.

(23) original volume fractions Interpolation high resolution volume fractions Convert particle volumes to volume fraction high resolution volume fractions with particle volumes Iso-surface extraction mesh Ray-tracing rendering. Figure 5: Rendering flowchart.. Figure 8 (a)-(c) show a sequence of selected frames of throwing a sphere into a tank of water at the 79th, 118th, and 138th minute of the simulation, respectively. The sphere has a initial horizontal 1.1 and vertical velocity of −0.8. The dimension of the computational grid is 8 in width and height and 12 in depth. The computation time depends upon the complexity of the scene and varies substantially throughout the simulation. It takes a long time to reach the second frame, when the sphere just hits the water surface, due to many iterations need to solve the pressure Poisson equation at each cycle. The computation is quick after the water surface is getting calm. For example, only 20 minutes are need to evolve from 12th frame to 24th frame. It costs about 10 minutes to render a frame. Figure 9 shows a sequence of selected frames of throwing a cube. The cube has a initial horizontal and vertical velocity of −2.5. The dimension of the computational grid is 7 in width and height and 6 in depth. The computation time for each frame is 30 seconds in average. The average rendering time for a frame is three minutes. 23.

(24) (a) Original volume fractions. (b) Interpolated volume fractions (c) Converting particles to volume fractions. Figure 6: Process to form a smooth fluid surface. (a) Frame 16. (b) Frame 32. (c) Frame 48. Figure 7: Liquid discharged from a pipe (30 × 24 × 24 gridl cells). Similar hybrid models in [27] and [25] permits a specific error when a solid object strike the surface at some angle to the veritical. Because they modeled the liquid volume by vertical columns, the force form the solid object is applied in a strictly vertical and the droplets escaping from the liquid surface can only be given an initial vertial velocity. Consequently, the resulting wave pattern and splash drops have no directional preference, which is incorrect. In this thesis, since the liquid volume is simulated by 3-D Navier-Stokes equations, the resulting waves naturally travel in the biased direction of the impact and the initial velocities of the droplets are interpolated near the surface. A feature characteristic of piecewise constant volume tracking methods is the unphysical and numerical creation of floatsam (“floating wreckage”) and jetsam (“jettisoned goods”) [26], as found in Figure 9 (c)-(f). These terms are appropriate for isolated, sub-mesh size material bodies that 24.

(25) (a) Frame 2. (b) Frame 12. (c) Frame 24. Figure 8: Selected frames of throwing a shpere into a tank of water (32 × 32 × 48 grid cells) separate from the main material body because of errors induced by the volume tracking algorithm. These material remnant tend to be ejected from interfaces in piecewise constant volume tracking methods when the flow has significant vorticity and shear near the interface [17]. The calculation time grows more than cubicly with the resolution of the grid. When the cell size decreases, the simulation time-step also need to decrease to satisfy the CFL condition (Section ??) and therefore more cycles are needed for a frame than in the grid of a lower resolution.. 6. Conclusion. To model liquid with dramatic deformation, a hybrid approach combining the volume tracking method and smoothed particle hydrodynamics (SPH) has been proposed. The 3-D Navier-Stokes equations are solved by the marker-and-cell (MAC) method and a density volume representing the liquid is evolved over time and space by the volume-of-fluid (VOF) method. The VOF method uses a finite-volume discretization and maintaining a value for the liquid volume in each grid cell. The advection of the volumes is calculated by using cell-face fluxes, in which the liquid volume that leaves one cell is exactly the same as the liquid volume that enters the adjacent cells. The proposed water splash modeling is a major improvement over previous hybrid volume-particle methods. Rules are proposed to identify regions with drastic motion and to position particles adequately. To add fidelity of the particle trajectory, the dynamics of the particles are governed by SPH, while previous methods treated splash as interaction-less particles. A smooth implicit surface 25.

(26) (a) Frame 19. (b) Frame 71. (c) Frame 103. (d) Frame 127. (e) Frame 159. (f) Frame 197. Figure 9: Selected frames of throwing a cube into a tank of water (28 × 24 × 28 grid cells). is extracted from the liquid density volume which interpolates the volume fractions in a higher resolution and converts particle volumes to volume fractions.. Modeling liquid by volume fractions is much more efficient and economical than traditional particle systems. Besides, the high realism of the animation is enabled by the coupling of pressure and velocity in 3-D Navier-Stokes equations. Potentially under-resolved subtle features are captured by smoothed particles that are adequately placed in the neighborhood of the highly-deformable regions. Close scrutiny reveals that SPH can provide more physically-correct particle motions because it approximates the Navier-Stokes equations. Consequently, it can continue to take into account the coupling between pressure and velocity terms, which is lacked in other particle simulation methods. Particle volumes are incorporated into the interpolated volume fractions before an unified iso-surface is extracted. As a result the liquid surface is more smooth than those of previous researches which rendered the splash by coating the particles with a field functions or as hard spheres. 26.

(27) References [1] N. Adabala and S. Manohar. Techniques for realistic visualization of fluids: A survey. Computer Graphics Forum, 21(1):65–81, 2002. [2] G. Ashraf and K Wong. Dust and water splashing models for hopping figures. Journal of Visualization and Computer Animation, 10(4):193–213, 1999. [3] D.S. Balsara.. Von-Neumann stability analysis of Smoothed Particle Hydrodynamics–. suggestions for optimal algorithms. Journal of Computational Physics, 121(2):357, 1995. [4] G.K. Batchelor. An Introduction to Fluid Dynamics. Cambridge University Press, Cambridge, 1967. [5] M. Desbrun and M.-P. Gascuel. Smoothed particles: A new paradigm for animating highly deformable bodies. In Proceedings of The 6th Eurographics Workshop on Animation and Simulation, pages 61–76, 1996. [6] D. Enright, R. Fedkiw, J. Ferziger, and I. Mitchell. A hybrid particle level set method for improved interface capturing. Journal of Computational Physics, 2002. To appear. [7] D. Enright, S. Marschner, and R. Fedkiw. Animation and rendering of complex water surfaces. In Proceedings of ACM SIGGRAPH ’02, 2002. To appear. [8] R. Fedkiw, H.W. Jensen, and J. Stam. Visual simulation of smoke. In Proceedings of ACM SIGGRAPH ’01, pages 15–22, 2001. [9] N. Foster and R. Fedkiw. Practical animation of liquids. In Proceedings of ACM SIGGRAPH ’01, pages 23–30, 2001. [10] N. Foster and D. Metaxas. Realistic animation of liquids. Graphical Model and Image Processing, 58(5):471–483, 1996. [11] N. Foster and D. Metaxas. Modelling the motion of a hot, turbulent gas. In Proceedings of ACM SIGGRAPH ’97, pages 181–188, 1997. 27.

(28) [12] G.Miller. Globular dynamics: A connected particle system for animating viscous fluids. Computer & Graphics, 13(3):305–309, 1989. [13] M.E. Goss. A real-time particle system for display of ship wakes. IEEE Computer Graphics and Applications, 10(3):30–35, 1990. [14] D. Gueyffier, A. Nadim, J. Li, R. Scardovelli, and S. Zaleski. Volume of fluid interface tracking with smoothed surface stress methods for three-dimensional flows. Journal of Computational Physics, 152:423–456, 1999. [15] F.H. Harlow and J.E. Welch. Numerical calculation of time-dependent viscous incompressible flow of fluid with free surface. Physics of Fluids, 8(12):2182–2189, 1965. [16] C.W. Hirt and B.D. Nichols. Volume of Fluid (VOF) method for the dynamics of free boundaries. Journal of Computational Physics, 39:201–225, 1981. [17] D.B. Kothe and W.J. Rider. Comments on modeling interfacial flows with Volume-of-Fluid methods. Technical Report LA-UR-94-3384, Los Alamos National Laborator, 1995. [18] A. Kunimatsu, Y. Watanabe, H. Fujii, T. Saito, K. Hiwada, and H. Ueki. Fast simulation and rendering techniques for fluid objects. Computer Graphics Forum, 20(3):57–67, 2001. [19] A.T. Layton and M. van de Panne. A numerically efficient and stable algorithm for animating water waves. The Visual Computer, 18(1):41–53, 2002. [20] W.E. Lorensen and H.E. Cline. Marching cubes: a high resolution 3D surface construction algorithm. Computer Graphics, 21(4):163–169, 1987. [21] J.J. Monaghan. Smoothed particle hydrodynamics. Annual Review of Astronomy and Astrophysics, 30:543–574, 1992. [22] J.J. Monaghan. Simulating free surface flows with SPH. Journal of Computational Physics, 110:399–406, 1994. 28.

(29) [23] J.J. Monaghan and J.C. Lattanzio. A refined method for astrophysical problems. Astron. Astrophys, 149:135–143, 1985. [24] J.P. Morris, P.J. Fox, and Y. Zhu. Modeling low reynolds number incompressible flows using SPH. Journal of Computational Physics, 136(1):214–226, 1997. [25] D. Mould and Y.H Yang. Modeling water for computer graphics. Computers & Graphics, 21(6):801–814, 1997. [26] W.F. Noh and P.R. Woodward. Slic (simple line interface method). In A.I. van de Vooren and P.J. Zandbergen, editors, Lecture Notes in Physics, volume 59, pages 330–340. Springer Verlag, Berlin, 1976. [27] J.F. O’Brien and J.K. Hodgins. Dynamic simulation of splashing fluids. In Computer Animation ’95, pages 198–205, 1995. [28] L.A. Pars. A treatise on analytical dynamics. Wiley, New York, 1965. [29] W.T. Reeves. Particle systems–a technique for modeling a class of fuzzy objects. ACM Transactions on Graphics, 2(2):91–108, 1983. [30] R. Scardovelli and S. Zaleski. Analytical relations connecting linear interfaces and volume fractions in rectangular grids. Journal of Computational Physics, 164(1):228–237, 2000. [31] J. Stam. Stable fluids. In Alyn Rockwood, editor, Proceedings of SIGGRAPH 99, pages 121–128, 1999. [32] D. Stora, P.-O. Agliati, M.-P. Gascuel, F. Neyret, and J.-D. Gascuel. Animating lava flows. In Proceedings of Graphics Interface ’99, pages 203–210, 1999. [33] G.D. Yngve, J.F. O’Brien, and J.K. Hodgins. Animating explosions. In Proceedings of ACM SIGGRAPH ’00, pages 29–36, 2000. 29.

(30) [34] D.L. Youngs. Time dependent multi-material flow with large fluid distortion. In K.W. Morton and M.J. Baines, editors, Numerical methods for fluid dunamics, pages 273–285. Academic Press, 1982.. 30.

(31)