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

時間  2019-12-05
標籤 fluid engine development 學習 筆記 光滑 粒子 流體動力學 简体版
原文   原文鏈接

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

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

經典核函數

SPH算法涉及到「光滑核」的概念，能夠這樣理解這個概念，粒子的屬性都會「擴散」到周圍，而且隨着距離的增長影響逐漸變小，這種隨着距離而衰減的函數被稱爲「光滑核」函數，最大影響半徑爲「光滑核半徑」。ide

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

SPH插值

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

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

相關代碼實現以下ui

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

咱們能夠看到插值函數依賴於密度，由於粒子的位置在每一個時間步長都會改變，而密度也隨之在每一個時間步長都會改。spa

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

梯度算子

相似於以前的插值，梯度能用相似的方法得到code

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個粒子將被施加不一樣的力。根據牛頓第三運動定律，每個做用力都有一個相等且相反的做用力component

爲解決這個問題，須要修改梯度實現。

書所使用的公式是 \(\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是經典核函數，實線是原核函數，虛線是一階偏導，點線是二階導

爲解決這個問題，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提出一種預測-校訂模型，消除了對聲速的依賴。這個新的模型將在下一篇筆記中闡述。

演示模擬結果

相關文章
相關標籤/搜索
每日一句
    每一个你不满意的现在，都有一个你没有努力的曾经。
本站公眾號
   歡迎關注本站公眾號,獲取更多信息
聯繫我們 最近搜索 最新文章 沪ICP备13005482号-10 MyBatis教程 SQL 教程 MySQL教程 Java 教程 Thymeleaf 教程 Hibernate教程 Spring教程 Redis教程