# 《Fluid Engine Development》 学习笔记3-光滑粒子流体动力学

## 经典核函数

SPH算法涉及到“光滑核”的概念，可以这样理解这个概念，粒子的属性都会“扩散”到周围，并且随着距离的增加影响逐渐变小，这种随着距离而衰减的函数被称为“光滑核”函数，最大影响半径为“光滑核半径”。

## SPH插值

SPH插值的基本思想是通过查找附近的粒子来测量任意给定位置的任何物理量。它是一个加权平均，权重是质量乘以核函数除以相邻粒子的密度。

Vector3D CalfFluidEngine::SphSystemData3::Interpolate(const Vector3D & origin, const std::vector<Vector3D>& values) const
{
Vector3D sum = Vector3D::zero;
auto& d = GetDensities();
const double m = GetParticleMass();

GetNeighborSearcher()->ForEachNearbyPoint(
origin, _kernelRadius, [&](size_t i, const Vector3D& neighborPosition)
{
double dist = Vector3D::Distance(origin,neighborPosition);
double weight = m / d[i] * kernel(dist);
sum += weight * values[i];
}
);

return sum;
}

double CalfFluidEngine::SphStandardKernel3::operator()(double distance) const
{
if (distance * distance >= h2) {
return 0.0;
}
else {
double x = 1.0 - distance * distance / h2;
return 315.0 / (64.0 * kPiD * h3) * x * x * x;
}
}

void CalfFluidEngine::PointHashGridSearcher3::ForEachNearbyPoint(const Vector3D & origin, double radius, const std::function<void(size_t, const Vector3D&)>& callback) const
{
if (_buckets.empty()) {
return;
}

size_t nearbyKeys[8];
getNearbyKeys(origin, nearbyKeys);

for (int i = 0; i < 8; i++) {
const auto& bucket = _buckets[nearbyKeys[i]];
size_t numberOfPointsInBucket = bucket.size();

for (size_t j = 0; j < numberOfPointsInBucket; ++j) {
size_t pointIndex = bucket[j];
double rSquared = (_points[pointIndex] - origin).SquareMagnitude();
callback(pointIndex, _points[pointIndex]);
}
}
}
}


void CalfFluidEngine::SphSystemData3::UpdateDensities()
{
auto& p = GetPositions();
auto& d = GetDensities();
const double m = GetParticleMass();

tbb::parallel_for(
tbb::blocked_range<size_t>(0, GetNumberOfParticles()),
[&](const tbb::blocked_range<size_t> & b) {
for (size_t i = b.begin(); i != b.end(); ++i)
{
double sum = SumOfKernelNearby(p[i]);
d[i] = m * sum;
}
});
}

double CalfFluidEngine::SphSystemData3::SumOfKernelNearby(const Vector3D & origin) const
{
double sum = 0.0;
GetNeighborSearcher()->ForEachNearbyPoint(
origin, _kernelRadius, [&](size_t, const Vector3D& neighborPosition) {
double dist = Vector3D::Distance(origin, neighborPosition);
sum += kernel(dist);
});



## 梯度算子

Vector3D CalfFluidEngine::SphSystemData3::GradientAt(size_t i, const std::vector<double>& values) const
{
Vector3D sum;
auto& p = GetPositions();
auto& d = GetDensities();
const auto& neighbors = GetNeighborLists()[i];
Vector3D origin = p[i];
const double m = GetParticleMass();

for (size_t j : neighbors) {
Vector3D neighborPosition = p[j];
double dist = Vector3D::Distance(origin, neighborPosition);
if (dist > kEpsilonD) {
Vector3D dir = (neighborPosition - origin) / dist;
sum += m * values[i] / d[j] *
}
}

return sum;
}

Vector3D ...::Gradient(double distance, const Vector3D & directionToParticle) const
{
return -firstDerivative(distance) * directionToParticle;
}


Vector3D CalfFluidEngine::SphSystemData3::GradientAt(size_t i, const std::vector<double>& values) const
{
Vector3D sum;
auto& p = GetPositions();
auto& d = GetDensities();
const auto& neighbors = GetNeighborLists()[i];
Vector3D origin = p[i];
const double m = GetParticleMass();

for (size_t j : neighbors) {
Vector3D neighborPosition = p[j];
double dist = Vector3D::Distance(origin, neighborPosition);
if (dist > kEpsilonD) {
Vector3D dir = (neighborPosition - origin) / dist;
sum += d[i] * m *
(values[i] / (d[i] * d[i]) + values[j] / (d[j] * d[j])) *
}
}

return sum;
}


## 拉普拉斯算子

double CalfFluidEngine::SphSystemData3::LaplacianAt(size_t i, const std::vector<double>& values) const
{
double sum = 0.0;
auto& p = GetPositions();
auto& d = GetDensities();
const auto& neighbors = GetNeighborLists()[i];
Vector3D origin = p[i];
const double m = GetParticleMass();

for (size_t j : neighbors) {
Vector3D neighborPosition = p[j];
double dist = Vector3D::Distance(origin, neighborPosition);
sum += m * values[j]  / d[j] * kernel.Laplacian(dist);
}

return sum;
}

double ...::Laplacian(double distance) const
{
return secondDerivative(distance);
}


double CalfFluidEngine::SphSystemData3::LaplacianAt(size_t i, const std::vector<double>& values) const
{
double sum = 0.0;
auto& p = GetPositions();
auto& d = GetDensities();
const auto& neighbors = GetNeighborLists()[i];
Vector3D origin = p[i];
const double m = GetParticleMass();

for (size_t j : neighbors) {
Vector3D neighborPosition = p[j];
double dist = Vector3D::Distance(origin, neighborPosition);
sum += m * (values[j] - values[i]) / d[j] * kernel.Laplacian(dist);
}

return sum;
}


## 主体代码结构

class SphSystemSolver3 : public ParticleSystemSolver3
{
public:
SphSystemSolver3();
virtual ~SphSystemSolver3();
void SetViscosityCoefficient(
double newViscosityCoefficient) {
_viscosityCoefficient = std::max(newViscosityCoefficient, 0.0);
}
void SetPseudoViscosityCoefficient(
double newPseudoViscosityCoefficient) {
_pseudoViscosityCoefficient
= std::max(newPseudoViscosityCoefficient, 0.0);
}
void SetTimeStepLimitScale(double newScale) {
_timeStepLimitScale = std::max(newScale, 0.0);
}
std::shared_ptr<SphSystemData3> GetSphData() const;
protected:
virtual void accumulateForces(double timeIntervalInSeconds) override;
virtual void onTimeStepStart(double timeStepInSeconds) override;
virtual void onTimeStepEnd(double timeStepInSeconds) override;
virtual unsigned int getNumberOfSubTimeSteps(
double timeIntervalInSeconds) const override;
private:
void accumulateViscosityForce();
void accumulatePressureForce(double timeStepInSeconds);
void computePressure();
void accumulatePressureForce(
const std::vector<Vector3D>& positions,
const std::vector<double>& densities,
const std::vector<double>& pressures,
std::vector<Vector3D>& pressureForces);
void computePseudoViscosity(double timeStepInSeconds);

//! Exponent component of equation - of - state(or Tait's equation).
double _eosExponent = 7.0;

//! Speed of sound in medium to determin the stiffness of the system.
//! Ideally, it should be the actual speed of sound in the fluid, but in
//! practice, use lower value to trace-off performance and compressibility.
double _speedOfSound = 100.0;

//! Negative pressure scaling factor.
//! Zero means clamping. One means do nothing.
double _negativePressureScale = 0.0;

double _viscosityCoefficient = 0.01;

//Scales the max allowed time-step.
double _timeStepLimitScale = 1.0;

//! Pseudo-viscosity coefficient velocity filtering.
//! This is a minimum "safety-net" for SPH solver which is quite
//! sensitive to the parameters.
double _pseudoViscosityCoefficient = 10.0;
};


SPH系统相比正常的粒子动画系统，重写了accumulateForces函数和onTimeStepStart函数以及onTimeStepEnd函数，分别用以添加粘度压力计算，更新密度，抑制噪声

void CalfFluidEngine::SphSystemSolver3::accumulateForces(double timeIntervalInSeconds)
{
ParticleSystemSolver3::accumulateForces(timeIntervalInSeconds);
accumulateViscosityForce();
accumulatePressureForce(timeIntervalInSeconds);
}


void CalfFluidEngine::SphSystemSolver3::onTimeStepStart(double timeStepInSeconds)
{
auto particles = GetSphData();

particles->UpdateDensities();
}


void CalfFluidEngine::SphSystemSolver3::onTimeStepEnd(double timeStepInSeconds)
{
computePseudoViscosity(timeStepInSeconds);
}


## 计算压强

inline double computePressureFromEos(
double density,
double targetDensity,
double eosScale,
double eosExponent,
double negativePressureScale) {
// Equation of state
// (http://www.ifi.uzh.ch/vmml/publications/pcisph/pcisph.pdf)
double p = eosScale / eosExponent
* (std::pow((density / targetDensity), eosExponent) - 1.0);

return p;
}


inline double computePressureFromEos(
double density,
double targetDensity,
double eosScale,
double eosExponent,
double negativePressureScale) {
// Equation of state
// (http://www.ifi.uzh.ch/vmml/publications/pcisph/pcisph.pdf)
double p = eosScale / eosExponent
* (std::pow((density / targetDensity), eosExponent) - 1.0);

// Negative pressure scaling
if (p < 0) {
p *= negativePressureScale;
}

return p;
}


void CalfFluidEngine::SphSystemSolver3::computePressure()
{
auto particles = GetSphData();
size_t numberOfParticles = particles->GetNumberOfParticles();
auto& d = particles->GetDensities();
auto& p = particles->GetPressures();

// See Equation 9 from
// http://cg.informatik.uni-freiburg.de/publications/2007_SCA_SPH.pdf
const double targetDensity = particles->GetDensity();
const double eosScale
= targetDensity * (_speedOfSound * _speedOfSound) / _eosExponent;

tbb::parallel_for(
tbb::blocked_range<size_t>(0, numberOfParticles),
[&](const tbb::blocked_range<size_t> & b) {
for (size_t i = b.begin(); i != b.end(); ++i)
{
p[i] = computePressureFromEos(
d[i],
targetDensity,
eosScale,
_eosExponent,
_negativePressureScale);
}
});
}


## 计算压力

$f_{p} = - m \frac{\nabla p}{\rho}$

void CalfFluidEngine::SphSystemSolver3::accumulatePressureForce(const std::vector<Vector3D>& positions, const std::vector<double>& densities, const std::vector<double>& pressures, std::vector<Vector3D>& pressureForces)
{
auto particles = GetSphData();
size_t numberOfParticles = particles->GetNumberOfParticles();

double mass = particles->GetParticleMass();
const double massSquared = mass * mass;

tbb::parallel_for(
tbb::blocked_range<size_t>(0, numberOfParticles),
[&](const tbb::blocked_range<size_t> & b) {
for (size_t i = b.begin(); i != b.end(); ++i)
{
const auto& neighbors = particles->GetNeighborLists()[i];
for (size_t j : neighbors) {
double dist = Vector3D::Distance(positions[i], positions[j]);

if (dist > kEpsilonD) {
Vector3D dir = (positions[j] - positions[i]) / dist;
pressureForces[i] -= massSquared
* (pressures[i] / (densities[i] * densities[i])
+ pressures[j] / (densities[j] * densities[j]))
}
}
}
});
}


## 计算粘度

void CalfFluidEngine::SphSystemSolver3::accumulateViscosityForce()
{
auto particles = GetSphData();
size_t numberOfParticles = particles->GetNumberOfParticles();
auto& x = particles->GetPositions();
auto& v = particles->GetVelocities();
auto& d = particles->GetDensities();
auto& f = particles->GetForces();

double mass = particles->GetParticleMass();
const double massSquared = mass * mass;

tbb::parallel_for(
tbb::blocked_range<size_t>(0, numberOfParticles),
[&](const tbb::blocked_range<size_t> & b) {
for (size_t i = b.begin(); i != b.end(); ++i)
{
const auto& neighbors = particles->GetNeighborLists()[i];
for (size_t j : neighbors) {
double dist = Vector3D::Distance(x[i],x[j]);

f[i] += _viscosityCoefficient * massSquared
* (v[j] - v[i]) / d[j]
* kernel.Laplacian(dist);
}
}
});
}


## 降低噪声

void CalfFluidEngine::SphSystemSolver3::computePseudoViscosity(double timeStepInSeconds)
{
auto particles = GetSphData();
size_t numberOfParticles = particles->GetNumberOfParticles();
auto& x = particles->GetPositions();
auto& v = particles->GetVelocities();
auto& d = particles->GetDensities();

const double mass = particles->GetParticleMass();

std::vector<Vector3D> smoothedVelocities(numberOfParticles);

tbb::parallel_for(
tbb::blocked_range<size_t>(0, numberOfParticles),
[&](const tbb::blocked_range<size_t> & b) {
for (size_t i = b.begin(); i != b.end(); ++i)
{
double weightSum = 0.0;
Vector3D smoothedVelocity;

const auto& neighbors = particles->GetNeighborLists()[i];
for (size_t j : neighbors) {
double dist = Vector3D::Distance(x[i],x[j]);
double wj = mass / d[j] * kernel(dist);
weightSum += wj;
smoothedVelocity += wj * v[j];
}

double wi = mass / d[i];
weightSum += wi;
smoothedVelocity += wi * v[i];

if (weightSum > 0.0) {
smoothedVelocity /= weightSum;
}

smoothedVelocities[i] = smoothedVelocity;
}
});

double factor = timeStepInSeconds * _pseudoViscosityCoefficient;
factor = Clamp(factor, 0.0, 1.0);

tbb::parallel_for(
tbb::blocked_range<size_t>(0, numberOfParticles),
[&](const tbb::blocked_range<size_t> & b) {
for (size_t i = b.begin(); i != b.end(); ++i)
{
v[i] = Lerp(
v[i], smoothedVelocities[i], factor);
}
});
}


## 声速参数和时间步长

$\Delta t _{v} =\frac{ \lambda _{v} h}{c_{s}} ,\Delta t_{f} = \lambda_{f}\sqrt{\frac{hm}{F_{Max}}}, \Delta \leq(\Delta t_{v},\Delta t_{f})$

$\lambda_{v},\lambda_{f}$是2个预设好的标量，大概0.25~0.4之间，$F_{max}$ 是力向量的最大大小

## 演示模拟结果

posted @ 2019-06-19 00:34  寂灭万乘  阅读(497)  评论(0编辑  收藏