斷斷續續花了一個月,終於把這本書的一二兩章啃了下來,理解流體模擬的理論彷佛不難,不管是《Fluid Simulation for Computer Graphics》仍是《計算流體力學基礎及其應用》都能很好幫助程序員去理解這些原理,可在缺少實踐狀況下,這種對原理的理解其實跟死記硬背沒什麼區別。《Fluid Engine Development》提供了一個實現完成的流體模擬引擎以及它的編程實現原理,充分幫助程序員經過編程實現流體動畫引擎,以此完成流體模擬學習的第一步。這不,早在今年一月就嚷嚷研究學習流體模擬卻苦苦掙扎沒法入門的我,在抄着做者的代碼看着做者的書的狀況,終於實現了一個流體模擬引擎。我終於能夠自信地說本身已經入門了流體模擬o(╥﹏╥)o。html
我學習流體模擬的本質緣由,是由於迷戀上了複雜的流體動畫,渴望經過編程來創造它們。流體模擬主要是指結合流體模擬的物理現象、方程和計算機圖形學的方法來模擬海面、海浪、煙霧等場景。git
流體模擬的基本方法可分爲三類: 基於紋理變換的流體模擬、基於二維高度場網格的流體模擬以及基於真實物理方程的流體模擬。程序員
基於紋理變換的流體模擬只須要對水面紋理進行法向擾動後、繪製水面的倒影(反射)以及繪製水底的狀況(透射)便可繪製出通常的水面效果。 但這種方法因爲其根本上沒有水面網格,因此水面起伏的繪製效果不明顯。github
基於二維高度場的網格流體模擬方法把水面表示成爲一個連續的平面網格,而後生成一系列對應於這張網絡的連續的高度紋理-稱爲高度圖。接着每一個網格頂點對應於一個高度圖的像素,做爲水面高度,從而表示出整個水面。算法
以上2者方法相對簡單,但真實度不高,本質上只是hack,基於真實物理方程的流體模擬才能製做目前所能呈如今計算機上最真實的流體動畫。儘管這種方法,受限於太高的計算量,一直處於離線渲染領域,很難在實時遊戲中使用,但隨着實時流體模擬技術的進步以及硬件條件的進步,在遊戲中使用這種方法並不遙遠,實際上,這甚至已經成爲了現實。看看這個視頻,實時流體模擬已經能作得很好了。編程
流體模擬的基本方法主要有三種,一種使用粒子(拉格朗日方法),一種使用網格(歐拉方法),第三種則是2者的混合。本文主要講基於粒子的流體模擬。數組
流體動畫首先是動畫,動畫的本質是在給定時間序列播放一系列圖像,由此咱們能夠編寫幀的結構體和動畫的抽象類。咱們建立新的動畫時,只要編寫類繼承Animation,再重寫onUpdate函數,在指定幀更新圖像便可。網絡
struct Frame final { public: int index = 0; double timeIntervalInSeconds = 1.0 / 60.0; double GetTimeInSeconds() const; void Advance(); void Advance(int delta); }; class Animation{ public: void Update(const Frame& frame); protected: //********************************************** //the function is called from Animation:Update(); //this function and implement its logic for updating the animation state. //********************************************** virtual void onUpdate(const Frame& frame) = 0; }; void Animation::Update(const Frame & frame){ onUpdate(frame); } double Frame::GetTimeInSeconds() const{ return index * timeIntervalInSeconds; } void Frame::Advance(){ ++index; } void Frame::Advance(int delta){ index += delta; }
流體動畫屬於基於物理的動畫,它不一樣於正常的動畫,正常的動畫並不依賴外界的輸入以及其餘幀的數據和狀態,它們大可能是根據時間狀態的改變來播放不一樣的圖像,或是經過以時間爲輸入的函數(例如sin(t)圖像)來切換狀態。基於物理的動畫正好相反,當前幀的數據和狀態徹底是由上一幀決定的。由此編寫了PhysicsAnimation類app
class PhysicsAnimation : public Animation{ public: double GetCurrentTime() const { return _currentTime; } private: void initialize(); void timeStep(double timeIntervalInSeconds); Frame _currentFrame; double _currentTime = 0.0; protected: virtual void onUpdate(const Frame& frame) override final; //********************************************** //the function is called from Animation:timeStep(double); //This function is called for each time-step; //********************************************** virtual void onTimeStep(double timeIntervalInSeconds) = 0; //********************************************** //the function is called from Animation:initialize(); //Called at frame 0 to initialize the physics state.; //********************************************** virtual void onInitialize() = 0; }; void PhysicsAnimation::initialize(){ onInitialize(); } void PhysicsAnimation::onUpdate(const Frame & frame){ if (frame.index > _currentFrame.index) { if (_currentFrame.index < 0) { initialize(); } int numberOfFrames = frame.index - _currentFrame.index; for (int i = 0; i < numberOfFrames; ++i) { timeStep(frame.timeIntervalInSeconds); } _currentFrame = frame; } }
咱們能夠看到PhysicsAnimation重寫了onUpdate函數,採用漸近更新每一幀dom
咱們使用粒子模擬流體,那就不可避免的要使用粒子系統
如下是我我的流體引擎粒子系統動畫系統的部分頭文件
class ParticleSystemData3 { public: typedef std::vector<double> ScalarData; typedef std::vector<Vector3D> VectorData; ParticleSystemData3(); virtual ~ParticleSystemData3(); explicit ParticleSystemData3(size_t numberOfParticles); void Resize(size_t newNumberOfParticles); size_t GetNumberOfParticles() const; size_t AddVectorData(const Vector3D& initialVal = Vector3D::zero); size_t AddScalarData(double initialVal = 0.0); const std::vector<Vector3D>& GetPositions() const; std::vector<Vector3D>& GetPositions(); const std::vector<Vector3D>& GetVelocities() const; std::vector<Vector3D>& GetVelocities(); const std::vector<Vector3D>& GetForces() const; std::vector<Vector3D>& GetForces(); void AddParticle( const Vector3D& newPosition, const Vector3D& newVelocity, const Vector3D& newForce = Vector3D::zero); void AddParticles( const std::vector<Vector3D>& newPositions, const std::vector<Vector3D>& newVelocities, const std::vector<Vector3D>& newForces); const std::vector<double>& ScalarDataAt(size_t idx) const; std::vector<double>& ScalarDataAt(size_t idx); const std::vector<Vector3D>& VectorDataAt(size_t idx) const; std::vector<Vector3D>& VectorDataAt(size_t idx); double GetParticleRadius() const; void SetParticleRadius(double newRadius); double GetParticleMass() const; virtual void SetParticleMass(double newMass); void BuildNeighborSearcher(double maxSearchRadius); void BuildNeighborLists(double maxSearchRadius); const std::shared_ptr<PointNeighborSearcher3> & GetNeighborSearcher() const; const std::vector<std::vector<size_t>>& GetNeighborLists() const; protected: size_t _positionIdx; size_t _velocityIdx; size_t _forceIdx; size_t _numberOfParticles = 0; double _radius = 1e-3; double _mass = 1e-3; std::vector<ScalarData> _scalarDataList; std::vector<VectorData> _vectorDataList; std::shared_ptr<PointNeighborSearcher3> _neighborSearcher; std::vector<std::vector<size_t>>_neighborLists; }; class ParticleSystemSolver3 : public PhysicsAnimation{ private: void timeStepStart(double timeStepInSeconds); void timeStepEnd(double timeStepInSeconds); void timeIntegration(double timeIntervalInSeconds); void resolveCollision(); void updateCollider(double timeStepInSeconds); void updateEmitter(double timeStepInSeconds); ParticleSystemData3::VectorData _newPositions; ParticleSystemData3::VectorData _newVelocities; protected: void onTimeStep(double timeIntervalInSeconds) override; virtual void onInitialize() override; //********************************************** //the function is called in ParticleSystemSolver3:timeStepStart(double); // Called when a time-step is about to begin; //********************************************** virtual void onTimeStepStart(double timeStepInSeconds); //********************************************** //the function is called in ParticleSystemSolver3:timeStepEnd(double); // Called when a time-step is about to end; //********************************************** virtual void onTimeStepEnd(double timeStepInSeconds); //********************************************** //the function is called in ParticleSystemSolver3:onTimeStep(double); //accumulate forces //********************************************** virtual void accumulateForces(double timeIntervalInSeconds); std::shared_ptr<ParticleSystemData3> _particleSystemData; std::shared_ptr<VectorField3> _wind; Vector3D _gravity = Vector3D(0.0, -9.8, 0.0); double _dragCoefficient = 1e-4; std::shared_ptr<Collider3> _collider; double _restitutionCoefficient = 0.0; std::shared_ptr<ParticleEmitter3> _emitter; };
主要看粒子系統重寫的onTimeStep函數,它用以更新粒子系統,它依次調用了5個函數。
accumulateForces用以累加做用於粒子上的力,timeIntegration用以更新粒子速度和位置.resolveCollision處理碰撞,timeStepStart,timeStepEnd分別用做每幀邏輯更新先後的事件調用以及內部數據更新
void ParticleSystemSolver3::onTimeStep(double timeIntervalInSeconds) { timeStepStart(timeIntervalInSeconds); accumulateForces(timeIntervalInSeconds); timeIntegration(timeIntervalInSeconds); resolveCollision(); timeStepEnd(timeIntervalInSeconds); }
具體實現以下
void CalfFluidEngine::ParticleSystemSolver3::timeStepStart(double timeStepInSeconds) { auto& forces = _particleSystemData->GetForces(); tbb::parallel_for( tbb::blocked_range<size_t>(0, forces.size()), [&](const tbb::blocked_range<size_t> & b) { for (size_t i = b.begin(); i != b.end(); ++i) forces[i] = Vector3D::zero; }); updateCollider(timeStepInSeconds); updateEmitter(timeStepInSeconds); size_t n = _particleSystemData->GetNumberOfParticles(); _newPositions.resize(n); _newVelocities.resize(n); onTimeStepStart(timeStepInSeconds); } void CalfFluidEngine::ParticleSystemSolver3::accumulateForces(double timeIntervalInSeconds) { size_t n = _particleSystemData->GetNumberOfParticles(); auto& forces = _particleSystemData->GetForces(); auto& velocities = _particleSystemData->GetVelocities(); auto& positions = _particleSystemData->GetPositions(); const double mass = _particleSystemData->GetParticleMass(); tbb::parallel_for( tbb::blocked_range<size_t>(0, n), [&](const tbb::blocked_range<size_t> & b) { for (size_t i = b.begin(); i != b.end(); ++i){ // Gravity Vector3D force = mass * _gravity; // Wind forces Vector3D relativeVel = velocities[i] - _wind->Sample(positions[i]); force += -_dragCoefficient * relativeVel; forces[i] += force; } }); } void CalfFluidEngine::ParticleSystemSolver3::timeIntegration(double timeIntervalInSeconds) { size_t n = _particleSystemData->GetNumberOfParticles(); auto& forces = _particleSystemData->GetForces(); auto& velocities = _particleSystemData->GetVelocities(); auto& positions = _particleSystemData->GetPositions(); const double mass = _particleSystemData->GetParticleMass(); tbb::parallel_for( tbb::blocked_range<size_t>(0, n), [&](const tbb::blocked_range<size_t> & b) { for (size_t i = b.begin(); i != b.end(); ++i) { Vector3D& newVelocity = _newVelocities[i]; newVelocity = velocities[i]; newVelocity.AddScaledVector(forces[i] / mass, timeIntervalInSeconds); Vector3D& newPosition = _newPositions[i]; newPosition = positions[i]; newPosition.AddScaledVector(newVelocity, timeIntervalInSeconds); } }); } void CalfFluidEngine::ParticleSystemSolver3::resolveCollision() { if (_collider != nullptr) { size_t numberOfParticles = _particleSystemData->GetNumberOfParticles(); const double radius = _particleSystemData->GetParticleRadius(); tbb::parallel_for( tbb::blocked_range<size_t>(size_t(0), numberOfParticles), [&](const tbb::blocked_range<size_t> & b) { for (size_t i = b.begin(); i != b.end(); ++i) { _collider->ResolveCollision( radius, _restitutionCoefficient, &_newPositions[i], &_newVelocities[i]); } }); } } void CalfFluidEngine::ParticleSystemSolver3::timeStepEnd(double timeStepInSeconds) { size_t n = _particleSystemData->GetNumberOfParticles(); auto& positions = _particleSystemData->GetPositions(); auto& velocities = _particleSystemData->GetVelocities(); tbb::parallel_for( tbb::blocked_range<size_t>(0, n), [&](const tbb::blocked_range<size_t> & b) { for (size_t i = b.begin(); i != b.end(); ++i) { positions[i] = _newPositions[i]; velocities[i] = _newVelocities[i]; } }); onTimeStepEnd(timeStepInSeconds); }
上一節咱們提到了碰撞,這個須要特別列出一節來說。
碰撞在物理模擬乃至物理引擎領域都是一個很是複雜的問題,《Fluid Engine Development》對於這一塊並無進行深刻,碰撞檢測只是簡單檢測碰撞點是否在幾何表面的內部,以及碰撞點,粒子的中心間的距離是否小於粒子的半徑。
至於碰撞的處理,主要是改變粒子的速度,和把粒子修正到正確的位置。
粒子的正確位置爲碰撞點+碰撞點法線 * 粒子半徑,也就是緊貼碰撞面表面的位置
至於粒子的速度,這裏把粒子相對碰撞面的速度拆解切向速度和法向速度
通過碰撞,法向速度變動爲 \(- R V_{n}\) (R是恢復係數,Vn是法向速度)
切向速度變動爲 \(max\left( 1 - u \frac{\left|\Delta V_{n}\right|}{V_{t}} \right) V_{t}\)
切向速度的變動你們或許會有疑惑,這裏列出推導過程
假定\(\Delta V_{n} = V_{n}^{new} - V_{n} = (- R - 1)V_{n}\)
已知\(\Delta V_{t} = a_{t}\Delta t = F_{f}/m \Delta t = u F_{n}/m \Delta t\)
又上類推得\(\Delta V_{n} = a_{n}\Delta t = F_{n}/m \Delta t\)
因此\(\Delta V_{t} = u \Delta V_{n}\)
由於新切向速度必然小於原速度,因此\(\Delta V_{t} = min(u\Delta V_{n},V_{t})\)
可得 \(V_{t}^{new} =max\left( 1 - u \frac{\left|\Delta V_{n}\right|}{V_{t}} \right) V_{t}\)
具體代碼實現以下
struct ColliderQueryResult final { double distance; Vector3D point; Vector3D normal; Vector3D velocity; }; void CalfFluidEngine::Collider3::ResolveCollision( double radius, double restitutionCoefficient, Vector3D * position, Vector3D * velocity) { ColliderQueryResult res; getClosestPoint(_surface, *position, &res); if (isPenetrating(res, *position, radius)) { Vector3D targetNormal = res.normal; Vector3D targetPoint = res.point + radius * targetNormal; Vector3D colliderVelAtTargetPoint = res.velocity; Vector3D relativeVel = *velocity - colliderVelAtTargetPoint; double normalDotRelativeVel = Vector3D::Dot(targetNormal, relativeVel); Vector3D relativeVelN = normalDotRelativeVel * targetNormal; Vector3D relativeVelT = relativeVel - relativeVelN; if (normalDotRelativeVel < 0.0) { Vector3D deltaRelativeVelN = (-restitutionCoefficient - 1.0) * relativeVelN; relativeVelN *= -restitutionCoefficient; // Apply friction to the tangential component of the velocity // From Bridson et al., Robust Treatment of Collisions, Contact and // Friction for Cloth Animation, 2002 // http://graphics.stanford.edu/papers/cloth-sig02/cloth.pdf if (relativeVelT.SquareMagnitude() > 0.0) { double frictionScale = std::max( 1.0 - _frictionCoeffient * deltaRelativeVelN.Magnitude() / relativeVelT.Magnitude(), 0.0); relativeVelT *= frictionScale; } *velocity = relativeVelN + relativeVelT + colliderVelAtTargetPoint; } *position = targetPoint; } } bool CalfFluidEngine::Collider3::isPenetrating( const ColliderQueryResult & colliderPoint, const Vector3D & position, double radius) { return _surface->IsInside(position) || colliderPoint.distance < radius; }
粒子系統有一個名爲updateEmitter的函數須要注意一下,它在timeStepStart函數和onInitialize()函數中被調用,用於隨機生成粒子系統的粒子數據。不過在實例demo中,粒子只生成了一次,畢竟每幀從新生成粒子是很是耗時的。粒子發射器的核心代碼以下。
void CalfFluidEngine::VolumeParticleEmitter3::onUpdate(double currentTimeInSeconds, double timeIntervalInSeconds) { auto particles = GetTarget(); if (particles == nullptr) { return; } if (_numberOfEmittedParticles > 0 && _isOneShot) { return; } std::vector<Vector3D> newPositions; std::vector<Vector3D> newVelocities; std::vector<Vector3D> Forces; emit(particles, &newPositions, &newVelocities); particles->AddParticles(newPositions, newVelocities, Forces); } void CalfFluidEngine::VolumeParticleEmitter3::emit(const std::shared_ptr<ParticleSystemData3>& particles, std::vector<Vector3D>* newPositions, std::vector<Vector3D>* newVelocities) { if (_surface == nullptr) return; _surface->Update(); BoundingBox3D region = _bounds; if (_surface->IsBounded()) { BoundingBox3D surfaceBox = _surface->GetBoundingBox(); region.lowerCorner = max(region.lowerCorner, surfaceBox.lowerCorner); region.upperCorner = min(region.upperCorner, surfaceBox.upperCorner); } const double j = GetJitter(); const double maxJitterDist = 0.5 * j * _spacing; if (_allowOverlapping || _isOneShot) { _pointGenerator->ForEachPoint(region, _spacing, [&](const Vector3D& point) { Vector3D randomDir = uniformSampleSphere(random(_rng), random(_rng)); Vector3D offset = maxJitterDist * randomDir; Vector3D candidate = point + offset; if (_surface->SignedDistance(candidate) <= 0.0) { if (_numberOfEmittedParticles < _maxNumberOfParticles) { newPositions->push_back(candidate); ++_numberOfEmittedParticles; } else { return false; } } return true; }); } else { PointHashGridSearcher3 neighborSearcher( Vector3<size_t>( kDefaultHashGridResolution, kDefaultHashGridResolution, kDefaultHashGridResolution), 2.0 * _spacing); if (!_allowOverlapping) { neighborSearcher.Build(particles->GetPositions()); } _pointGenerator->ForEachPoint(region, _spacing, [&](const Vector3D& point) { Vector3D randomDir = uniformSampleSphere(random(_rng), random(_rng)); Vector3D offset = maxJitterDist * randomDir; Vector3D candidate = point + offset; if (_surface->SignedDistance(candidate) <= 0.0 && (!_allowOverlapping &&!neighborSearcher.HasNearbyPoint(candidate, _spacing))) { if (_numberOfEmittedParticles < _maxNumberOfParticles) { newPositions->push_back(candidate); neighborSearcher.Add(candidate); ++_numberOfEmittedParticles; } else { return false; } } return true; }); } newVelocities->resize(newPositions->size()); tbb::parallel_for( tbb::blocked_range<size_t>(0, newVelocities->size()), [&](const tbb::blocked_range<size_t> & b) { for (size_t i = b.begin(); i != b.end(); ++i) (*newVelocities)[i] = velocityAt((*newPositions)[i]); }); } Vector3D CalfFluidEngine::VolumeParticleEmitter3::velocityAt(const Vector3D & point) const { Vector3D r = point - _surface->transform.GetTranslation(); return _linearVel + Vector3D::Cross(_angularVel,r) + _initialVel; } void CalfFluidEngine::BccLatticePointGenerator::ForEachPoint(const BoundingBox3D & boundingBox, double spacing, const std::function<bool(const Vector3D&)>& callback) const { double halfSpacing = spacing / 2.0; double boxWidth = boundingBox.GetWidth(); double boxHeight = boundingBox.GetHeight(); double boxDepth = boundingBox.GetDepth(); Vector3D position; bool hasOffset = false; bool shouldQuit = false; for (int k = 0; k * halfSpacing <= boxDepth && !shouldQuit; ++k) { position.z = k * halfSpacing + boundingBox.lowerCorner.z; double offset = (hasOffset) ? halfSpacing : 0.0; for (int j = 0; j * spacing + offset <= boxHeight && !shouldQuit; ++j) { position.y = j * spacing + offset + boundingBox.lowerCorner.y; for (int i = 0; i * spacing + offset <= boxWidth; ++i) { position.x = i * spacing + offset + boundingBox.lowerCorner.x; if (!callback(position)) { shouldQuit = true; break; } } } hasOffset = !hasOffset; } }
咱們看到上節中粒子發射器的代碼實現裏使用了PointHashGridSearcher3這個類,這是我爲何說粒子發生器須要注意的緣由,這是個很是精妙的算法,它將3D點映射到整數,它被這樣設計的意義是爲了加速搜索。
它將長方體的包圍盒區域劃分紅多個長方形塊,將3D點的座標除以長方形塊的邊長,將它轉換成三維整形數(Vector3
代碼以下
size_t CalfFluidEngine::PointHashGridSearcher3::GetHashKeyFromBucketIndex(const Vector3<size_t>& bucketIndex) const { Vector3<size_t> wrappedIndex = bucketIndex; wrappedIndex.x = bucketIndex.x % _resolution.x; wrappedIndex.y = bucketIndex.y % _resolution.y; wrappedIndex.z = bucketIndex.z % _resolution.z; return (wrappedIndex.z * _resolution.y + wrappedIndex.y) * _resolution.x + wrappedIndex.x; } Vector3<size_t> CalfFluidEngine::PointHashGridSearcher3::GetBucketIndex(const Vector3D & position) const { Vector3<size_t> bucketIndex; bucketIndex.x = static_cast<size_t>( std::floor(position.x / _gridSpacing)); bucketIndex.y = static_cast<size_t>( std::floor(position.y / _gridSpacing)); bucketIndex.z = static_cast<size_t>( std::floor(position.z / _gridSpacing)); return bucketIndex; }
PointHashGridSearcher3內部有一個二維數組,或者能夠稱它爲桶數組,桶數組的鍵值就是由3D座標轉換而來的key,表明一個包圍盒區域內的長方形塊。桶數組存儲所對應長方塊內的全部粒子在粒子數組中的索引。
當須要查找鄰域粒子時,就能夠經過粒子座標所對應的key值查找獲取粒子所在長方塊附近其餘長方形塊的key,從而獲取粒子所在長方形塊和粒子附近長方形塊內其餘粒子。
獲取附近key值代碼實現以下
void CalfFluidEngine::PointHashGridSearcher3::getNearbyKeys(const Vector3D & position, size_t * nearbyKeys) const { Vector3<size_t> originIndex = GetBucketIndex(position), nearbyBucketIndices[8]; for (int i = 0; i < 8; i++) { nearbyBucketIndices[i] = originIndex; } if ((originIndex.x + 0.5f) * _gridSpacing <= position.x) { nearbyBucketIndices[4].x += 1; nearbyBucketIndices[5].x += 1; nearbyBucketIndices[6].x += 1; nearbyBucketIndices[7].x += 1; } else { nearbyBucketIndices[4].x -= 1; nearbyBucketIndices[5].x -= 1; nearbyBucketIndices[6].x -= 1; nearbyBucketIndices[7].x -= 1; } if ((originIndex.y + 0.5f) * _gridSpacing <= position.y) { nearbyBucketIndices[2].y += 1; nearbyBucketIndices[3].y += 1; nearbyBucketIndices[6].y += 1; nearbyBucketIndices[7].y += 1; } else { nearbyBucketIndices[2].y -= 1; nearbyBucketIndices[3].y -= 1; nearbyBucketIndices[6].y -= 1; nearbyBucketIndices[7].y -= 1; } if ((originIndex.z + 0.5f) * _gridSpacing <= position.z) { nearbyBucketIndices[1].z += 1; nearbyBucketIndices[3].z += 1; nearbyBucketIndices[5].z += 1; nearbyBucketIndices[7].z += 1; } else { nearbyBucketIndices[1].z -= 1; nearbyBucketIndices[3].z -= 1; nearbyBucketIndices[5].z -= 1; nearbyBucketIndices[7].z -= 1; } for (int i = 0; i < 8; i++) { nearbyKeys[i] = GetHashKeyFromBucketIndex(nearbyBucketIndices[i]); } }
獲取粒子附近其餘粒子代碼以下
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) { ...... } } }
下一節將進入流體力學的部分,在那以前先重溫或學習一下關於矢算場的數學知識,若是不瞭解的話,會對Navier-Stocks方程的一些數學符號感到疑惑。
《Fluid Engine Development》有介紹這部分的知識,更詳細的你們能夠閱讀另外一位博主推薦的《失算場論札記》,不過這本書我淘寶也買不到,在學校圖書館借了看,後來花9塊錢買了矢量分析與場論-工程數學-第四版
\(\frac{\partial f}{\Delta } = \frac{f(x + \Delta,y,z) - f(x,y,z)}{ \Delta}\)
\(\nabla f(x,y,z) = (\frac{\partial f}{\partial x} , \frac{\partial f}{\partial y} , \frac{\partial f}{\partial z} )\)
\(\nabla\) 是梯度算子,用以測量標量的變化率和變換方向,
梯度是標量場增加最快的方向,梯度的長度是這個最大的變化率
梯度本意是一個向量(矢量),當某一函數在某點處沿着該方向的方向導數取得該點處的最大值,即函數在該點處沿方向變化最快,變化率最大(爲該梯度的模)
對程序開發者而言,它是由一階偏導組成的矢量
當咱們須要知道函數在特定方向上的變化率,經過將梯度與該方向的向量點乘,咱們就能夠計算出函數在這一點的方向導數 $\frac{\partial f}{\partial n} =\nabla f \cdot n $
\(\nabla \cdot f(x) = (\frac{\partial }{\partial x} , \frac{\partial }{\partial y} , \frac{\partial }{\partial z} ) \cdot f(x)\)
$\nabla \cdot $ 是散度算子,散度就是通量密度,通量即經過一個面的某物理量
在流體力學中,速度場的散度是體積膨脹率,散度(divergence)算子用於描述向量場中在某點的彙集或發散程度
對程序開發者,這就是梯度算子與向量場的點乘
流體模擬中,流體每每有着難以壓縮的性質,一般認爲速度場的散度爲零
\(\nabla\times f(x) = (\frac{\partial }{\partial x} , \frac{\partial }{\partial y} , \frac{\partial }{\partial z} ) \times f(x)\)
\(\nabla\times\) 是旋度算子,旋度就是環量密度,以水旋渦爲例,沿着圓周關於速度的線積分就是環量,環量除以圓面積,而後讓面積無限小,比值就是旋度,旋度用以描述某個點的旋轉程度
對程序開發者,這就是梯度算子與向量場的叉積
$\nabla^{2} f(x) = \nabla \cdot \nabla f(x) =\frac{\partial^{2} f}{\partial x^{2} } +\frac{\partial^{2} f}{\partial y^{2} } +\frac{\partial^{2} f}{\partial z^{2} } $
$\nabla^{2} $是拉普拉斯算子,又稱拉氏算子,它的本質是梯度的散度
拉普拉斯算子一般用於描述一個函數中某點的值與它鄰域的平均值之差,或者說用於描述某點 在它的鄰域中是「多麼不同凡響」
拉普拉斯算子能夠檢測圖像的邊緣和峯值,也能夠用以模糊矢量場的尖銳部分(拉普拉斯算子 加上原始矢量場)
如今開始進入正題,前面講的都是編程相關的細節,如今開始流體力學的部分。
流體運動的本質是源於多種力的組合。
在保證密度不變的狀況下,流體受到重力,壓力,粘度力三種力的做用。
重力不須要我太多介紹,主要是壓力和粘度力
風是高壓區吹到低壓區的,一樣的規則適用於流體。
壓力實際上是壓力梯度力,當你潛水時,隨着深度的增長,收到的壓力也隨之增長。這種壓力差會沿着深度在與重力相反的方向產生向上的壓力。正由於這種力,水能夠保持其狀態而不是收縮
假設一個立方體大小的流體,單位大小是L,密度爲 \(\rho\), 假設壓力沿x軸方向,左右之間的壓力差是\(\Delta p\),壓力大小就等於\(\Delta p L^{2}\),因此 $F_{p} = -\Delta p L^{2}= ma_{p} = L^{3}\rho a_{p} $
可推導得 \(a_{p} = \frac{-\Delta p}{\rho L}\),假定立方體很是小,可得 $a_{p} = -\frac{\partial p}{\partial x} \frac{1}{\rho} $
接着咱們將它推導到三維上 ,可得 $a_{p} = -\frac{\nabla p}{\rho} $
因此$ a = g -\frac{\nabla p}{\rho} + \cdots$
流體模擬中粘度相似於阻尼,它試圖模糊2點間的速度差別,它致使液體有了必定的粘稠感。正好,拉普拉斯算子能作到模糊矢量場。
\(V_{new} = V + \Delta t u \nabla^{2}V\) u 是粘度係數
由上式能夠獲得粘度影響的加速度$ a_{V} =u \nabla^{2}V$
因此$ a = g -\frac{\nabla p}{\rho} + u \nabla^{2}V$ 這就是 著名的 納維斯托克斯(Navier--Stokes, NS)方程