基于粒子的流体模拟交互式应用
In contrast to Eulerian grid-based approaches, the particle-based approach makes mass conservation equations and convection terms dispensable which reduces the complexity of the simulation. In addition, the particles can directly be used to render the surface of the fluid. We propose methods to track and visualize the free surface using point splatting and marching cubes-based surface reconstruction. Our animation method is fast enough to be used in interactive systems and to allow for user interaction with models consisting of up to 5000 particles.
逼真的动画流体可以为虚拟手术模拟器或电脑游戏等交互式应用程序增加实质性的真实感。本文提出了一种基于光滑粒子流体动力学（SPH）的交互方法来模拟具有自由表面的流体。该方法是Desbrun基于sph的技术的扩展，用于动画高度可变形的物体。我们通过直接从Navier-Stokes方程推导力密度场，并通过添加一项来模拟表面张力效应，将方法推向流体模拟。 与基于欧拉网格的方法相比，基于粒子的方法省去了质量守恒方程和对流项，从而降低了模拟的复杂性。此外，颗粒可以直接用于渲染流体的表面。我们提出了使用点飞溅和基于行进立方体的表面重建来跟踪和可视化自由表面的方法。我们的动画方法足够快，可以在交互式系统中使用，并允许用户与多达5000个粒子组成的模型进行交互。
介绍
动机
Fluids (i.e. liquids and gases) play an important role in every day life. Examples for fluid phenomena are wind, weather, ocean waves, waves induced by ships or simply pouring of a glass of water. As simple and ordinary these phenomena may seem, as complex and difficult it is to simulate them. Eventhough Computational Fluid Dynamics (CFD) is a well established research area with a long history, there are still many open research problems in the field. The reason for the complexity of fluid behavior is the complex interplay of various phenomena such as convection, diffusion, turbulence and surface tension. Fluid phenomena are typically simulated off-line and then visualized in a second step e.g. in aerodynamics or optimization of turbines or pipes with the goal of being as accurate as possible.
流体（即液体和气体）在日常生活中起着重要作用。流体现象的例子有风、天气、海浪、船舶引起的波浪或仅仅是倒一杯水。这些现象看似简单和普通，但模拟它们却是复杂和困难的。尽管计算流体力学（CFD）是一个有着悠久历史的成熟研究领域，但该领域仍存在许多尚未解决的问题。流体行为复杂的原因是对流、扩散、湍流和表面张力等各种现象的复杂相互作用。流体现象通常是离线模拟，然后在第二步中可视化，例如在空气动力学或涡轮机或管道的优化中，目标是尽可能准确.
Less accurate methods that allow the simulation of fluid effects in real-time open up a variety of new applications. In the fields mentioned above real-time methods help to test whether a certain concept is promising during the design phase. Other applications for real-time simulation techniques for fluids are medical simulators, computer games or any type of virtual environment
不太精确的方法允许实时模拟流体效应，开辟了各种新的应用。在上述领域中，实时方法有助于在设计阶段测试某个概念是否有希望。流体实时模拟技术的其他应用包括医疗模拟器、电脑游戏或任何类型的虚拟环境
Related Work 关联工作
Computational Fluid Dynamics has a long history. In 1822 Claude Navier and in 1845 George Stokes formulated the famous Navier-Stokes Equations that describe the dynamics of fluids. Besides the Navier-Stokes equation which describes conservation of momentum, two additional equations namely a continuity equation describing mass conservation and a state equation describing energy conservation are needed to simulate fluids. Since those equations are known and computers are available to solve them numerically, a large number of methods have been proposed in the CFD literature to simulate fluids on computers.
Since about two decades, special purpose fluid simulation techniques have been developed in the field of computer graphics. In 1983 T. Reeves17 introduced particle systems as a technique for modeling a class of fuzzy objects. Since then both, the particle-based Lagrangian approach and the gridbased Eulerian approach have been used to simulate fluids in computer graphics. Desbrun and Cani2 and Tonnesen22 use particles to animate soft objects. Particles have also been used to animate surfaces7, to control implicit surfaces23 and to animate lava flows20. In recent years the Eulerian approach has been more popular as for the simulation of fluids in general18, water 4,3,21, soft objects14 and melting effects1. So far only a few techniques optimized for the use in interactive systems are available. Stam’s method18 that is grid based is certainly an important step towards real-time simulation of fluids. For the special case of fluids that can be represented by height fields, interactive animation techniques are available as well6. Here we propose a particle-based approach based on Smoothed Particle Hydrodynamics to animate arbitrary fluid motion.
计算流体力学有着悠久的历史。1822年克劳德·纳维尔和1845年Claude Navier（乔治·斯托克斯）提出了著名的描述流体动力学的纳维-斯托克斯方程 (N-S方程)。除了描述动量守恒的Navier-Stokes方程外，还需要两个附加方程，即描述质量守恒的连续性方程和描述能量守恒的状态方程来模拟流体。由于这些方程是已知的，并且有计算机可以对它们进行数值求解，因此在CFD文献中提出了大量在计算机上模拟流体的方法。
近二十年来，专用流体模拟技术在计算机图形学领域得到了发展。1983年，T. Reeves17引入粒子系统作为对一类模糊物体建模的技术。从那时起，基于粒子的拉格朗日方法和基于网格的欧拉方法都被用于计算机图形学中模拟流体。Desbrun和Cani2和Tonnesen22使用粒子来动画软物体。粒子也被用来使表面动起来，控制隐含的表面和使熔岩流动起来。近年来，欧拉方法在一般流体、水、软物体和熔化效应的模拟中更为流行。到目前为止，只有少数技术是为在交互式系统中使用而优化的。斯塔姆基于网格的方法无疑是朝着实时模拟流体迈出的重要一步。对于可以用高度场表示的流体的特殊情况，也可以使用交互式动画技术。本文提出了一种基于光滑粒子流体力学的基于粒子的方法来动画化任意流体运动。
Our Contribution SPH的贡献
We propose a method based on Smoothed Particles Hydrodynamics (SPH) to simulate fluids with free surfaces. Stam and Fiume first introduced SPH to the graphics community to depict fire and other gaseous phenomena19. Later, Desbrun used SPH to animate highly deformable bodies2.We extend his method focussing on the simulation of fluids. To this end, we derive the viscosity and pressure force fields directly from the Navier-Stokes equation and propose a way to model surface tension forces. For the purpose of interactivity, we designed new special purpose smoothing kernels. Surface tracking and surface rendering at interactive rates are difficult problems for which we describe possible solutions
提出了一种基于光滑粒子流体力学（SPH）的模拟自由表面流体的方法。斯塔姆和富姆首先将SPH引入图形界，用来描绘火和其他气体现象。后来，Desbrun用SPH制作了高度可变形的物体动画。我们将他的方法扩展到流体的模拟上。为此，我们直接从Navier-Stokes方程推导出黏度力场和压力力场，并提出了一种模拟表面张力的方法。为了实现交互，我们设计了新的专用平滑核。表面跟踪和表面渲染在交互速率是困难的问题，我们描述了可能的解决方案
Smoothed Particle Hydrodynamics （SPH）
Although Smoothed Particle Hydrodynamics (SPH) was developed by Lucy and Gingold and by Monaghan for the simulation of astrophysical problems, the method is general enough to be used in any kind of fluid simulation. For introductions to SPH we refer the reader to Monaghan or Münzel
虽然光滑粒子流体力学（SPH）是由Lucy和Gingold以及Monaghan为模拟天体物理问题而开发的，但该方法足够通用，可以用于任何类型的流体模拟。关于SPH的介绍，我们建议读者参考Monaghan或Münzel
SPHis an interpolation method for particle systems. With SPH, field quantities that are only defined at discrete particle locations can be evaluated anywhere in space. For this purpose, SPH distributes quantities in a local neighborhood of each particle using radial symmetrical smoothing kernels. According to SPH, a scalar quantity A is interpolated at location r by a weighted sum of contributions from all particles
spi是粒子系统的一种插值方法。使用SPH，仅在离散粒子位置定义的场量可以在空间的任何位置进行评估。为此，SPH利用径向对称平滑核在每个粒子的局部邻域内分配数量。根据SPH，标量a通过所有粒子贡献的加权和插值到位置r
where j iterates over all particles, mj is the mass of particle j, rj its position, ρj the density and Aj the field quantity at rj.
其中j遍历所有粒子，mj为粒子j的质量，rj为粒子j的位置，ρj为密度，Aj为rj处的场量。
The function W(r,h) is called the smoothing kernel with core radius h. Since we only use kernels with finite support, we use h as the radius of support in our formulation. If W is even (i.e. W(r,h)=W(−r,h)) and normalized, the interpolation is of second order accuracy. The kernel is normalized if
函数W（r,h）被称为平滑核，其核半径为h。由于我们只使用具有有限支持的核，因此我们在公式中使用h作为支持半径。如果W是偶数(即W（r,h)=W(- r,h)）并且归一化，则插值具有二阶精度。核是标准化的，如果
The particle mass and density appear in Eqn. because each particle i represents a certain volume Vi = mi / ρi. While the mass mi is constant throughout the simulation and, in our case, the same for all the particles, the density ρi varies and needs to be evaluated at every time step. Through substitution into Eqn. we get for the density at location r
粒子质量和密度用方程表示。因为每个粒子i代表一定的体积vi = mi/ρi。在整个模拟过程中，质量mi是恒定的，在我们的例子中，对所有粒子都是相同的，而密度ρi是变化的，需要在每个时间步长进行评估。通过代入到Eqn。我们得到了r点的密度
这里消元了，密度倍被约分了
In most fluid equations, derivatives of field quantities appear and need to be evaluated. With the SPH approach, such derivatives only affect the smoothing kernel. The gradient of A is simply
在大多数流体方程中，都会出现场量的导数，需要对其进行计算。在SPH方法中，这种导数只影响平滑核。a的梯度很简单
while the Laplacian of A evaluates to
而A的拉普拉斯函数等于
It is important to realize that SPH holds some inherent problems. When using SPH toderive fluid equations for particles, these equations are not guaranteed to satisfy certain physical principals such as symmetry of forces and conservation of momentum. The next section describes our SPH-based model and techniques to solve these SPH-related problems.
重要的是要认识到SPH存在一些固有的问题。当使用SPH推导粒子的流体方程时，这些方程不能保证满足某些物理原理，如力的对称性和动量守恒。下一节描述我们基于sph的模型和解决这些sph相关问题的技术。
Modelling Fluids with Particles 用颗粒模拟流体
In the Eulerian (grid based) formulation, isothermal fluids are described by a velocity field v, a density field ρ and a pressure field p. The evolution of these quantities over time is given by two equations. The first equation assures conservation of mass
在欧拉（基于网格的）公式中，等温流体由速度场v、密度场ρ和压力场p来描述。这些量随时间的演变由两个方程给出。第一个方程保证了质量守恒
while the Navier-Stokes equation formulates conservation of momentum
而Navier-Stokes方程则表述了动量守恒(注意 密度场ρ和压力场p别看错了)
where g is an external force density field and µ the viscosity of the fluid. Many forms of the Navier-Stokes equation appear in the literature. Eqn. represents a simplified version for incompressible fluids.
其中g为外力密度场，µ为流体粘度。文献中出现了许多形式的纳维-斯托克斯方程。表示不可压缩流体的简化版本。
The use of particles instead of a stationary grid simplifies these two equations substantially. First, because the number of particles is constant and each particle has a constant mass, mass conservation is guaranteed and Eqn. can be omitted completely. Second, the expression ∂v/∂t +v·∇v on the left hand side of Eqn. can be replaced by the substantial derivative Dv/Dt. Since the particles move with the fluid, the substantial derivative of the velocity field is simply the time derivative of the velocity of the particles meaning that the convective term v·∇v is not needed for particle systems.
用粒子代替固定网格大大简化了这两个方程。首先，因为粒子的数量是恒定的，每个粒子的质量是恒定的，所以质量守恒是有保证的，Eqn。可完全省略。第二，表达式∂v/∂t +v·∇v在Eqn的左边。可以用实质导数Dv/Dt代替。由于粒子随流体运动，速度场的实质导数就是粒子速度的时间导数，这意味着粒子系统不需要对流项v·∇v。
There are three force density fields left on the right hand side of Eqn. modeling pressure (−∇p), external forces (ρg) and viscosity (µ∇^2v). The sum of these force density fields f = −∇p+ρg+µ∇2v determines the change of momentum ρDv/Dt of the particles on the left hand side. For the acceleration of particle i we, thus, get
在Eqn的右侧有三个力密度场。模拟压力（−∇p）、外力（ρg）和粘度（µ∇^2v）。这些力密度场f =−∇p+ρg+µ∇2v的和决定了左边粒子的动量ρDv/Dt的变化。对于粒子i的加速度，我们得到
where vi is the velocity of particle i and fi and ρi are the force density field and the density field evaluated at the location of particle i, repectively. We will now describe how we model the force density terms using SPH
其中vi为质点I的速度，fi和ρi分别为力密度场和质点I位置处的密度场。现在我们将描述如何使用SPH对力密度项进行建模
Pressure 压力
Application of the SPH rule described in Eqn. to the pressure term −∇p yields
Eqn中描述的SPH规则的应用。到压力项−∇p产率
Unfortunately, this force is not symmetric as can be seen when only two particles interact. Since the gradient of the kernel is zero at its center, particle i only uses the pressure of particle j to compute its pressure force and vice versa. Because the pressures at the locations of the two particles are not equal in general, the pressure forces will not be symmetric. Different ways of symmetrization of Eqn. have been proposed in the literature. We suggest a very simple solution which we found to be best suited for our purposes of speed and stability
不幸的是，当只有两个粒子相互作用时，这种力不是对称的。由于核的梯度在其中心为零，粒子i只使用粒子j的压力来计算它的压力，反之亦然。因为两个粒子所在位置的压强一般不相等，所以压强不对称。方程的不同对称方式。已经在文献中提出。我们建议一个非常简单的解决方案，我们发现最适合我们的目的，速度和稳定性
The socomputed pressure force is symmetric because it uses the arithmetic mean of the pressures of interacting particles.
所计算的压力是对称的，因为它使用了相互作用粒子压力的算术平均值。
Since particles only carry the three quantities mass, position and velocity, the pressure at particle locations has to be evaluated first. This is done in two steps. Eqn. yields the density at the location of the particle. Then, the pressure can be computed via the ideal gas state equation
由于粒子只携带质量、位置和速度这三个量，因此必须首先计算粒子位置处的压力。这分两步完成。得到粒子所在位置的密度。然后，可以通过理想气体状态方程来计算压力
where k is a gas constant that depends on the temperature. In our simulations we use a modified version of Eqn. suggested by Desbrun^2
k是一个气体常数它取决于温度。在我们的模拟中，我们使用了一个改进版本的Eqn。由Desbrun^2提出
where ρ0 is the rest density. Since pressure forces depend on the gradient of the pressure field, the offset mathematically has not effect on pressure forces. However, the offset does influence the gradient of a field smoothed by SPH and makes the simulation numerically more stable.
其中ρ0是静止密度。由于压力与压力场的梯度有关，因此数学上的偏置对压力没有影响。然而，偏移量确实会影响SPH平滑场的梯度，使数值模拟更加稳定。
Viscosity 粘度
Application of the SPH rule to the viscosity term µ∇^2v again yields asymmetric forces
将SPH规则应用于黏度项µ∇^2v再次产生不对称力
because the velocity field varies from particle to particle. Since viscosity forces are only dependent on velocity differences and not on absolute velocities, there is a natural way to symmetrize the viscosity forces by using velocity differences:
因为速度场因粒子而异。由于黏性力只依赖于速度差，而不是绝对速度，有一个自然的方法来对称黏性力通过使用速度差：
A possible interpretation of Eqn. is to look at the neighbors of particle i from i’s own moving frame of reference. Then particle i is accelerated in the direction of the relative speed of its environment.
对公式的可能解释。就是从I自身的运动参照系中观察粒子I的邻居。然后粒子i在其环境的相对速度方向上加速。
Surface Tension 表面张力
We model surface tension forces (not present in Eqn.) explicitly based on ideas of Morris. Molecules in a fluid are subject to attractive forces from neighboring molecules. Inside the fluid these intermolecular forces are equal in all directions and balance each other. In contrast, the forces acting on molecules at the free surface are unbalanced. The net forces (i.e. surface tension forces) act in the direction of the surface normal towards the fluid. They also tend to minimize the curvature of the surface. The larger the curvature, the higher the force. Surface tension also depends on a tension coefficient σ which depends on the two fluids that form the surface
我们根据莫里斯Morris的思想明确地模拟了表面张力（不存在于方程中）。流体中的分子受到邻近分子的吸引力。在流体中，这些分子间的力在各个方向上都相等，并且相互平衡。相反，作用在自由表面的分子上的力是不平衡的。合力（即表面张力）作用于表面法向流体的方向。它们也倾向于最小化表面的曲率。曲率越大，力就越大。表面张力还取决于张力系数σ，而张力系数σ取决于形成表面的两种流体
The surface of the fluid can be found by using an additional field quantity which is 1 at particle locations and everywhere else. This field is called color field in the literature. For the smoothed color field we get:
流体的表面可以通过在粒子位置和其他地方使用额外的场量为1来找到。这个领域在文献中被称为色域。对于平滑的颜色场，我们得到：
The gradient field of the smoothed color field
平滑色场的渐变场（gradient field，这个很重要，是斥力的表现）
yields the surface normal field pointing into the fluid and the divergence of n measures the curvature of the surface
产生指向流体的表面法向场，n的散度测量表面的曲率
To distribute the surface traction among particles near the surface and to get a force density we multiply by a normalized scalar field δs = |n| which is non-zero only near the surface. For the force density acting near the surface we get
为了在表面附近的粒子之间分配表面引力，并得到力密度，我们乘以一个归一化标量场δs = |n|，它只在表面附近是非零的。对于作用在表面附近的力密度，我们得到
Evaluating n/|n| at locations where |n| is small causes numerical problems. We only evaluate the force if |n| exceeds a certain threshold.
在|n|较小的位置计算n/|n|会导致数值问题。我们只在b| n|超过一定阈值时才计算力。
External Forces 额外力
Our simulator supports external forces such as gravity, collision forces and forces caused by user interaction. These forces are applied directly to the particles without the use of SPH. When particles collide with solid objects such as the glass in our examples, we simply push them out of the object and reflect the velocity component that is perpendicular to the object’s surface
我们的模拟器支持外力，如重力、碰撞力和由用户交互引起的力。这些力直接作用于粒子而不使用SPH。当粒子与固体物体（如我们例子中的玻璃）碰撞时，我们只是将它们推出物体，并反射垂直于物体表面的速度分量
Smoothing Kernels 平滑核
Stability, accuracy and speed of the SPH method highly depend on the choice of the smoothing kernels. The kernels we use have second order interpolation errors because they are all even and normalized. In addition, kernels that are zero with vanishing derivatives at the boundary are conducive to stability. Apart from those constraints, one is free to design kernels for special purposes. We designed the following kernel
SPH方法的稳定性、精度和速度在很大程度上取决于光滑核的选择。我们使用的核函数有二阶插值误差，因为它们都是偶数且归一化的。此外，在边界处导数为零的核有利于稳定性。除了这些限制之外，人们可以自由地为特殊目的设计内核。我们设计了以下内核
其对应的Python代码为(可用于绘制图像):
import numpy as np
import matplotlib.pyplot as plt
def W_poly6(r, h):
"""Poly6 核函数"""
coef = 315 / (64 * np.pi * h**9)
return coef * (h**2 - r**2)**3 * (np.abs(r) <= h)
def grad_W_poly6(r, h):
"""Poly6 核函数的梯度"""
coef = -6 * 315 / (64 * np.pi * h**9)
return coef * r * (h**2 - r**2)**2 * (np.abs(r) <= h)
def laplacian_W_poly6(r, h):
"""Poly6 核函数的拉普拉斯"""
coef = 315 / (64 * np.pi * h**9)
return coef * (12 * r**2 - 6 * h**2) * (h**2 - r**2) * (np.abs(r) <= h)
def W_spiky(r, h):
"""Spiky 核函数"""
coef = 15 / (np.pi * h**6)
return coef * (h - np.abs(r))**3 * (np.abs(r) <= h)
def grad_W_spiky(r, h):
"""Spiky 核函数的梯度（结果取绝对值）"""
coef = -3 * 15 / (np.pi * h**6)
return coef * (h - np.abs(r))**2 * (np.abs(r) <= h)
def laplacian_W_viscosity(r, h):
"""Viscosity 核函数的拉普拉斯"""
coef = 15 / (2 * np.pi * h**3)
return coef * (h - np.abs(r)) * (np.abs(r) <= h)
# 平滑长度 h 和距离 r，x 轴在 [-2, 2] 范围内
h = 1.0
r = np.linspace(-2.0, 2.0, 1000) # 数据更密集，确保负半轴完整
# 计算核函数、梯度和拉普拉斯
poly6 = W_poly6(r, h)
grad_poly6 = np.abs(grad_W_poly6(r, h)) # 取绝对值
laplacian_poly6 = laplacian_W_poly6(r, h)
spiky = W_spiky(r, h)
grad_spiky = np.abs(grad_W_spiky(r, h)) # 取绝对值
viscosity_laplacian = laplacian_W_viscosity(r, h)
# 可视化
fig, axes = plt.subplots(1, 3, figsize=(12, 4))
# Poly6 核函数
axes[0].plot(r, poly6, label='Poly6 Kernel', linewidth=2)
axes[0].plot(r, grad_poly6, label='Gradient (|grad|)', linestyle='--')
axes[0].plot(r, laplacian_poly6, label='Laplacian', linestyle=':')
axes[0].set_title("W_poly6")
axes[0].set_xlim([-2, 2])
axes[0].legend()
# Spiky 核函数
axes[1].plot(r, spiky, label='Spiky Kernel', linewidth=2)
axes[1].plot(r, grad_spiky, label='Gradient (|grad|)', linestyle='--')
axes[1].set_title("W_spiky")
axes[1].set_xlim([-2, 2])
axes[1].legend()
# Viscosity 核函数拉普拉斯
axes[2].plot(r, viscosity_laplacian, label='Laplacian', linewidth=2)
axes[2].set_title("W_viscosity")
axes[2].set_xlim([-2, 2])
axes[2].legend()
plt.tight_layout()
plt.show()
and use it in all but two cases. An important feature of this simple kernel is that r only appears squared which means
除了两种情况外，都可以使用它。这个简单核函数的一个重要特征是r只出现平方，这意味着
The three smoothing kernels W_poly6,W_spiky and W_viscosity (from left to right) we use in our simulations. The thick lines show the kernels, the thin lines their gradients in the direction towards the center and the dashed lines the Laplacian. Note that the diagrams are differently scaled. The curves show 3-d kernels along one axis through the center for smoothing length h = 1.
我们在模拟中使用的三个平滑核W_poly6，W_spiky和W_viscosity（从左到右）。粗线表示核，细线表示向中心方向的梯度，虚线表示拉普拉斯函数。注意，图表的比例是不同的。当平滑长度为h = 1时，曲线显示沿一个轴穿过中心的三维核。
that is can be evaluated without computing square roots in distance computations. However, if this kernel is used for the computation of the pressure forces, particles tend to build clusters under high pressure. As particles get very close to each other, the repulsion force vanishes because the gradient of the kernel approaches zero at the center. Desbrun2 solves this problem by using a spiky kernel with a non vanishing gradient near the center. For pressure computations we use Debrun’s spiky kernel
它可以在距离计算中不计算平方根而求值。然而，如果用这个核来计算压力，粒子倾向于在高压下形成团簇。当粒子彼此非常接近时，斥力消失，因为核的梯度在中心趋近于零。Desbrun2通过使用在中心附近具有不消失梯度的尖状核解决了这个问题。对于压强的计算，我们使用Debrun的尖核
另外两种内核的公式如下, Wspiky(r,h)：Spiky pow3核函数
另外还有一个Spiky pow2的核函数: *(注意半径和直径的区别)
that generates the necessary repulsion forces. At the boundary where it vanishes it also has zero first and second derivatives.
这就产生了必要的斥力。在它消失的边界上一阶导数和二阶导数都为零。
Viscosity is a phenomenon that is caused by friction and, thus, decreases the fluid’s kinetic energy by converting it into heat. Therefore, viscosity should only have a smoothing effect on the velocity field. However, if a standard kernel is used for viscosity, the resulting viscosity forces do not always have this property. For two particles that get close to each other, the Laplacian of the smoothed velocity field (on which viscosity forces depend) can get negative resulting in forces that increase their relative velocity. The artifact appears in coarsely sampled velocity fields. In real-time applications where the number of particles is relatively low, this effect can cause stability problems. For the computation of viscosity forces we, thus, designed a third kernel:
粘度是一种由摩擦引起的现象，因此，它通过将流体的动能转化为热量来降低流体的动能。因此，粘度应该只对速度场起平滑作用。然而，如果使用标准核来表示粘度，则得到的粘度力并不总是具有这种性质。对于两个相互靠近的粒子，平滑速度场（黏性力所依赖的）的拉普拉斯函数可以为负，从而产生增加它们相对速度的力。伪影出现在粗采样的速度场中。在粒子数量相对较低的实时应用中，这种效应可能会导致稳定性问题。因此，对于粘性力的计算，我们设计了第三核：
W_viscosity(r,h)：Viscosity 核函数
whose Laplacian is positive everywhere with the following additional properties:
其拉普拉斯式处处为正，并具有以下附加性质：
The use of this kernel for viscosity computations increased the stability of the simulation significantly allowing to omit any kind of additional damping.
使用这个核粘度计算增加了稳定性的模拟显著允许省略任何类型的附加阻尼。
Simulation 模拟
For the integration of the Eqn. we use the Leap-Frog scheme. As a second order scheme for which the forces need to be evaluation only once, it best fits our purposes and in our examples allows time steps up to 10 milliseconds. For the examples we used constant time steps. We expect even better performance if adaptive time steps are used based on the Courant-Friedrichs-Lewy condition
对于Eqn的积分。我们使用跳蛙方案。作为只需要评估一次力的二阶方案，它最适合我们的目的，并且在我们的示例中允许时间步长达到10毫秒。对于例子，我们使用常数时间步长。如果使用基于Courant-Friedrichs-Lewy条件的自适应时间步长，我们期望有更好的性能
跳蛙模拟算法
常用于数值计算的时间步进方法，特别是在求解偏微分方程（如波动方程和流体动力学方程）时被广泛应用。
基本原理：
Leap-Frog Scheme 通过交错的时间步进，分别计算不同时间点上的物理量。在时间步长上，物理量的计算会“跳跃”到下一个时间点，从而得到新的值，这种方式类似于“蛙跳”的行为，因此得名 Leap-Frog。
在时间离散化时，Leap-Frog Scheme 采用中心差分法，具有以下特点：
时间上是二阶精度的。
在空间上也可以配合中心差分，达到二阶精度。
假设在时间 ttt 和 t+Δtt + \Delta tt+Δt 的两个位置上已知物理量，可以“跳跃”计算出 t+2Δtt + 2\Delta tt+2Δt 的值。
数学公式
Leap-Frog Scheme 主要应用于 时间导数项 的数值计算。例如，对于时间演化方程
离散化为：
Surface Tracking and Visualization 表面跟踪和可视化
从这部分开始，计算的干货没那么多了，主要是一些思维上的扩展想法
The color field cS and its gradient field n = ∇cS defined in section 3.3 can be used to identify surface particles and to compute surface normals. We identify a particle i as a surface particle if
3.3节中定义的颜色场cS及其梯度场n =∇cS可用于识别表面粒子和计算表面法线。我们把粒子i定义为表面粒子if
where l is a threshold parameter. The direction of the surface normal at the location of particle i is given by
其中l为阈值参数。粒子i位置的表面法线方向由
Point Splatting 泼洒
We now have a set of points with normals but without connectivity information. This is exactly the type of information needed for point splatting techniques24. However, these methods are designed to work with point clouds obtained from scanners that typically contain at least 10,000 to 100,000 points. We only use a few thousand particles a fraction of which are identified as being on the surface. Still surface splatting yields plausible results as shown in the results section.
We are currently working on ways to upsample the surface of the fluid. Hereby, the color field information of surface particles is interpolated to find locations for additional particles on the surface only used for rendering
我们现在有一组有法线但没有连通性信息的点。这正是点溅射技术所需要的信息类型。然而，这些方法设计用于从扫描仪获得的点云，通常包含至少10,000到100,000个点。我们只使用几千个粒子，其中一小部分被识别为在表面上。仍然表面溅射产生似是而非的结果，如结果部分所示。 我们目前正在研究对液体表面进行上采样的方法。在此，对表面粒子的颜色场信息进行插值，以查找仅用于渲染的表面上的附加粒子的位置
Marching Cubes
Another way to visualize the free surface is by rendering an iso surface of the color field cS. We use the marching cubes algorithm8 to triangulate the iso surface. In a grid fixed in space the cells that contain the surface are first identified.
另一种可视化自由曲面的方法是渲染颜色场cS的iso曲面。我们使用行进立方体算法对等值面进行三角剖分。在固定在空间中的网格中，首先识别包含表面的单元格。
We start searches from all the cells that contain surface particles and from there recursively traverse the grid along the surface. With the use of a hash table we make sure that the cells are not visited more than once. For each cell identified to contain the surface, the triangles are generated via a fast table lookup.
我们从包含表面粒子的所有细胞开始搜索，并从那里沿着表面递归地遍历网格。通过使用哈希表，我们可以确保单元格不会被访问超过一次。对于每个被识别为包含表面的单元格，三角形是通过快速表查找生成的。
Implementation 实现
Since the smoothing kernels used in SPH have finite support h, a common way to reduce the computational complexity is to use a grid of cells of size h. Then potentially interacting partners of a particle i only need to be searched in i’s own cell and all the neighboring cells. This technique reduces the time complexity of the force computation step from O(n^2) to O(nm), m being the average number of particles per grid cell.
由于SPH中使用的平滑核具有有限的支持h，因此降低计算复杂度的一种常用方法是使用大小为h的单元格。然后只需在i自己的单元格和所有邻近的单元格中搜索粒子i的潜在相互作用伙伴。该技术将力计算步骤的时间复杂度从O（n^2）降低到O(nm)， m为每个网格单元的平均粒子数。
With a simple additional trick we were able to speed up the simulation by an additional factor of 10. Instead of storing references to particles in the grid, we store copies of the particle objects in the grid cells (doubling memory consumption). The reason for the speed up is the proximity in memory of the information needed for interpolation which dramatically increases the cash hit rate. Further speedup might be possible through even better clustering using Hilbert space filling curves11. The data structure for fast neighbor searches is also used for surface tracking and rendering
通过一个简单的附加技巧，我们能够将模拟速度提高10倍。我们不是在网格中存储对粒子的引用，而是在网格单元中存储粒子对象的副本（内存消耗加倍）。加速的原因是内存中插值所需的信息的接近性，这极大地提高了现金命中率。通过使用希尔伯特空间填充曲线进行更好的聚类，可能会进一步加快速度。快速邻居搜索的数据结构也用于表面跟踪和渲染
这里看到一个方法，就是利用GPU去排序邻接Hash点，这种方法可以更好的利用GPU，通过GPU分离主线程的方式减少这一个部分的时间消耗
Results
The water in the glass shown in Fig. 3 is sampled with 2200 particles. An external rotational force field causes the fluid to swirl. The first image (a) shows the individual particles. For the second image (b), point splatting was used to render the free surface only. In both modes, the animation runs at 20 frames per second on a 1.8 GHz Pentium IV PC with a GForce 4 graphics card. The most convincing results are produced when the iso surface of the color field is visualized using the marching cubes algorithm as in image (c). However, in this mode the frame rate drops to 5 frames per second. Still this frame rate is significantly higher than the one of most off-line fluid simulation techniques and with the next generation of graphics hardware, real-time performance will be possible.
图3所示的玻璃杯中的水被取样了2200个颗粒。外部旋转力场使流体旋转。第一张图像(a)显示了单个粒子。对于第二张图像(b)，点飞溅仅用于渲染自由表面。在这两种模式下，动画在1.8 GHz的Pentium IV PC上以每秒20帧的速度运行，并带有GForce 4显卡。最令人信服的结果是当使用行军立方体算法可视化色场的iso表面时，如图(c)所示。然而，在这种模式下，帧率下降到每秒5帧。尽管如此，这个帧率仍然明显高于大多数离线流体模拟技术之一，并且随着下一代图形硬件的发展，实时性能将成为可能。
The image sequence shown in Fig. 4 demonstrates interaction with the fluid. Through mouse motion, the user generates an external force field that cause the water to splash. The free surface is rendered using point splatting while isolated particles are drawn as single droplets. The simulation with 1300 particles runs at 25 frames per second.
图4所示的图像序列显示了与流体的相互作用。通过鼠标移动，用户产生一个使水飞溅的外部力场。自由表面使用点飞溅渲染，而孤立的粒子被绘制为单个液滴。1300个粒子的模拟以每秒25帧的速度运行。
For the animation shown in Fig. 5 we used 3000 particles and rendered the surface with the marching cubes technique at 5 frames per second.
对于图5所示的动画，我们使用了3000个粒子，并以每秒5帧的速度使用行进立方体技术渲染表面。
这个思路扩展下去，现代计算机驱动，270k粒子，每秒60次计算，可以跑到 60帧+。
个人的认为SPH算法主要：
提供了一种思路 ，通过计算粒子间斥力的方式渲染斥力效果，这个斥力通过Smoothing kernel平滑核进行驱动，并且比较多种平滑核效果的优劣，最终落地模拟的效果。
然后通过物理量的估算效果，将粒子和密度、压力等物理量进行插值估算
最后通过粒子相互作用力，压力梯度，粘滞力等方法，计算运动学力。
同时使用一个时间积分，比如 Leap-Frog或者Verlet积分，或者牛顿更新粒子的位置和速度，推进这个流体仿真。
0 条评论