《Fluid Engine Development》學習筆記3-光滑粒子流體動力學

寂滅萬乘發表於2019-06-19

用粒子表示流體最熱門的方法就是就是光滑粒子流體動力學（Smoothed Particle Hydrodynamics (SPH).）

這種方法模糊了流體的邊界，用有限數量的粒子代表流體，該方法的基本思想是將視作連續的流體(或固體)用相互作用的質點組來描述，各個物質點上承載各種物理量，包括質量、速度等，通過求解質點組的動力學方程和跟蹤每個質點的運動軌道，求得整個系統的力學行為

經典核函式

SPH演算法涉及到“光滑核”的概念，可以這樣理解這個概念，粒子的屬性都會“擴散”到周圍，並且隨著距離的增加影響逐漸變小，這種隨著距離而衰減的函式被稱為“光滑核”函式，最大影響半徑為“光滑核半徑”。

書中提到的經典核函式有 $W_{std}(r) = \frac{315}{64\pi h^{3}}(1 -\frac{r^{2}}{h_{2}})^{3} (0 \leq r \leq h) $，其他情況為0

SPH插值

SPH插值的基本思想是通過查詢附近的粒子來測量任意給定位置的任何物理量。它是一個加權平均，權重是質量乘以核函式除以相鄰粒子的密度。

質量除以密度就是體積，因此這個插值，將更多的權重放在離原點更近的值上

相關程式碼實現如下

Vector3D CalfFluidEngine::SphSystemData3::Interpolate(const Vector3D & origin, const std::vector<Vector3D>& values) const
{
    Vector3D sum = Vector3D::zero;
    auto& d = GetDensities();
    SphStandardKernel3 kernel(_kernelRadius);
    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);

    const double queryRadiusSquared = radius * radius;

    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();
            if (rSquared <= queryRadiusSquared) {
                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;
    SphStandardKernel3 kernel(_kernelRadius);
    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];
    SphSpikyKernel3 kernel(_kernelRadius);
    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] *
                kernel.Gradient(dist, dir);
        }
    }

    return sum;
}

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

然而這種梯度的實現是不對稱的，相鄰的粒子可能會因為擁有不同的價值和密度而擁有不同的梯度，這也意味著2個粒子將被施加不同的力。根據牛頓第三運動定律，每一個作用力都有一個相等且相反的作用力

為解決這個問題，需要修改梯度實現。

書所使用的公式是 $\nabla \phi（x）= \rho _{j}m \sum_{j}(\frac{\phi_{i}}{\rho _{i} ^{2}} + \frac{\phi_{j}}{\rho _{j} ^{2}}) \nabla W(|x - x_{j}|)$

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];
    SphSpikyKernel3 kernel(_kernelRadius);
    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])) *
                kernel.Gradient(dist, dir);
        }
    }

    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];
    SphSpikyKernel3 kernel(_kernelRadius);
    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);
}

遺憾的是這般計算拉普拉斯運算元在即便所有場值都是相同的非零值時，也不會輸出零場

拉普拉斯正確的計算方法如下 $\nabla^{2} \phi（x）=m \sum_{j}(\frac{\phi_{j} - \phi_{i}}{\rho _{j} } ) \nabla^{2} W(|x - x_{j}|)$

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];
    SphSpikyKernel3 kernel(_kernelRadius);
    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;
}

Spiky核函式

梯度運算元是用來計算壓力梯度的，粒子太接近，壓力就會把粒子推開，然而經典核函式即使粒子越來越接近，也會出現壓力越來越小的情況，甚至還會出現負值

如下圖是原書中的圖,a是經典核函式，實線是原核函式，虛線是一階偏導，點線是二階導

《Fluid Engine Development》學習筆記3-光滑粒子流體動力學

為解決這個問題，Spiky核函式誕生了，如上圖b

公式為$W_{spiky}(r) = \frac{15}{\pi h^{3}}(1 -\frac{r^{3}}{h_{3}})^{3} (0 \leq r \leq h) $其他情況為0

我們插值獲取權重時使用經典核函式，計算拉普拉斯運算元和梯度時使用Spiky核函式

主體程式碼結構

這裡給出SPH系統的標頭檔案

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函式，分別用以新增粘度壓力計算，更新密度，抑制噪聲

以下是accumulateForces函式的程式碼結構

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

可以看到了相比粒子動畫，多了粘度和壓力的計算

以下是onTimeStepStart函式，用以更新粒子集合的密度

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

    particles->BuildNeighborSearcher(particles->GetKernelRadius());
    particles->BuildNeighborLists(particles->GetKernelRadius());
    particles->UpdateDensities();
}

以下是onTimeStepEnd函式

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

計算壓強

狀態方程（Equation-of-State ，EOS）描述了狀態變數間的關係，我們通過狀態方程 $p = \frac{\kappa}{\gamma}( \frac{\rho}{\rho_{0}}- 1)^{\gamma}$ 將密度對映為壓強

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;
}

觀察上公式，我們發現density 小於 targetDensity會出現負壓強的情況，而液體表面附近的確會出現密度過小的情況

為防止負壓強的引入，我們需要夾緊壓強，具體如下

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);
        }
    });
}

這裡注意到eosScale引數的計算，並不是我們想象中那樣隨便取個值，需要通過公式 $\kappa =\rho_{0} \frac{c_{s}}{\gamma} $ cs是流體中的聲速，實踐中可以用較低的值跟蹤效能。

計算壓力

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

回憶我們之前提到的梯度運算元計算方法，我們可以得到$f_{p}= m^{2} \sum_{j}(\frac{p_{i}}{\rho _{i} ^{2}} + \frac{p_{j}}{\rho _{j} ^{2}}) \nabla W(|x - x_{j}|)$

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;
    const SphSpikyKernel3 kernel(particles->GetKernelRadius());

    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]))
                        * kernel.Gradient(dist, dir);
                }
            }
        }
    });
}

計算粘度

粘度力公式為$f_{v} = - m \mu \nabla^{2}u$ 代入之前拉普拉斯運算元的計算方法，可得公式$\nabla^{2} \phi（x）=m^{2} \mu\sum_{j}(\frac{u_{j} - u_{i}}{\rho _{j} } ) \nabla^{2} W(|x - x_{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;
    const SphSpikyKernel3 kernel(particles->GetKernelRadius());

    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);
            }
        }
    });
}

降低噪聲

降低噪聲的方法很簡單，以引數_pseudoViscosityCoefficient線性插值速度場和加權平均速度即可

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();
    const SphSpikyKernel3 kernel(particles->GetKernelRadius());

    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);
        }
    });
}

聲速引數和時間步長

之前我們計算壓強時使用了聲速cs，為什麼會有聲速呢，因為在一個時間步長內，壓力傳播不能大於粒子核半徑h，而水中傳播的最快速度就是聲速，所以時間步長的理想步長是h/cs

最後，根據幾位科學家的研究成果，時間步長需要做如下的限制

$\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}$ 是力向量的最大大小

然後時間步長因為這種限制可能會非常小，導致巨大的計算成本，而且實際上我們也無法評估最大速度和最大力是多少

為從根本解決這個問題，Solenthaler 和Pajarola提出一種預測-校正模型，消除了對聲速的依賴。這個新的模型將在下一篇筆記中闡述。

演示模擬結果

《Fluid Engine Development》學習筆記3-光滑粒子流體動力學

《Fluid Engine Development》學習筆記2-基礎
2019-06-18
UIdev筆記
結構動力學教材-學習筆記
2024-06-19
筆記
vue學習筆記3-事件處理
2024-06-30
Vue筆記事件
【Go】Go語言學習筆記-3-包
2019-02-16
Go筆記
網路流學習筆記
2024-08-13
筆記
【學習筆記】網路流
2022-07-11
筆記
基金訓練營學習筆記3-股票基金
2021-01-04
筆記
數字人粒子系統學習筆記之一
2018-07-19
筆記
Java學習筆記09（IO流）
2020-12-10
Java筆記
[計算化學]分子動力學筆記
2024-03-22
筆記
JVM學習筆記——自動記憶體管理
2019-02-28
JVM筆記記憶體
力學筆記
2024-05-09
筆記
流體力學深度學習建模技術研究進展
2022-09-12
深度學習
李沐大佬-動手學深度學習筆記-注意力機制
2024-11-19
深度學習筆記
學習筆記1熱力學與統計物理
2018-11-13
筆記
Activiti工作流基礎學習筆記
2019-03-03
筆記
Node.js stream(流) 學習筆記(一)
2018-06-05
Node.js筆記
CAN匯流排協議學習筆記
2018-09-01
協議筆記
numpy的學習筆記\pandas學習筆記
2018-03-18
筆記
React學習手記3-事件處理
2019-03-02
React事件
流體力學8-3
2024-10-07
學習 PixiJS — 粒子效果
2019-01-22
JS
網路流最大流、最小割學習筆記
2021-01-21
筆記
強化學習-學習筆記13 | 多智慧體強化學習
2022-07-10
強化學習筆記智慧體
【軟體測試】學習筆記
2024-06-06
筆記
GO 學習筆記->結構體
2020-01-10
Go筆記結構體
HexMap學習筆記(八)——水體
2019-05-20
筆記
AC 自動機學習筆記
2024-07-28
筆記
動態規劃學習筆記
2024-08-03
動態規劃筆記
AC自動機學習筆記
2020-10-25
筆記
QT學習筆記4（動畫）
2020-11-30
QT筆記動畫
【學習筆記】數學
2024-07-15
筆記
學習筆記
2024-04-14
筆記
【java學習】JDK(Java Development Kit)
2018-09-28
JavaJDKdev
機器學習學習筆記
2021-06-01
機器學習筆記
《JAVA學習指南》學習筆記
2019-03-24
Java筆記
機器學習框架ML.NET學習筆記【9】自動學習
2019-06-10
機器學習框架筆記
004:ZYNQ_AXI匯流排學習筆記（1）
2021-06-05
筆記

《Fluid Engine Development》 學習筆記3-光滑粒子流體動力學

經典核函式

SPH插值

梯度運算元

拉普拉斯運算元

Spiky核函式

主體程式碼結構

計算壓強

計算壓力

計算粘度

降低噪聲

聲速引數和時間步長

演示模擬結果

相關文章

《Fluid Engine Development》學習筆記3-光滑粒子流體動力學