2D粗粒流体模拟

想做好看的东西

我试着推了一个方程，用来将理想不可压缩流体转换成粗粒流体。但由于我刚学不久，也不知道这个方法之前有没有人想到过。。。。总之先发出来吧，希望各位路过的大佬能多多指教╭(′▽`)╯

先来看下我的计算结果（512*512实时渲染，文末有完整Unity代码下载）

不可压流体NS方程：

模拟方法GPU Gems上有（https://developer.nvidia.com/gpugems/gpugems/part-vi-beyond-triangles/chapter-38-fast-fluid-dynamics-simulation-gpu）这里不细讲。

粗粒流体转换的方法与我上篇写粗粒NDF的核心思想相同，只不过由于NDF只需要计算当前时刻视线方向的结果，把方程简化了。所以还请不要在意“粗粒流体”这个名词，我只是擅自这样称呼它。

来讲我做了什么。

我认为现实中的流体内部压力变化并没有那么理想，所以想了个办法对压力场的梯度做了一个转换。为了导出我的方程，这里提出2个假设：

假设1：对于理想不可压缩流体，单位质量流体（单元格）内可看作包含了无穷多个流体微团。每个无穷小流体微团对应的梯度分量大小相等，方向以此单位质量流体的理想压力场梯度为峰值服从标准正态分布。

假设2：对于粗粒流体，单位质量流体内可看作包含了有限个（N_real个）流体微团，且与不可压流体内流体微团的性质相同。

首先，我们需要为每单位质量流体分配N_real个服从标准正态分布的随机数：

（i：0～1的随机数；R_gaus：服从标准正态分布的随机数）

然后把R_gaus映射成梯度分量的偏转角度θ：

每单位质量流体有N_real个流体微团，也就需要N_real个θ用来偏转梯度分量后相加得到此单位质量流体的粗粒梯度方向：

方向有了，接下来求大小。

先获得等间隔服从标准正态分布的偏转角度φ：

（i = 1，2，3。。。N_real）

然后偏转N_real次相加再除以梯度大小，得到梯度分量与梯度的大小之比scale：

方向和大小都求出来了，结合在一起就是粗粒流体的梯度了：

也可以这么理解：对每一单位质量流体掷N_real次骰子，在时间和空间的整体上构建出一个稳定的结果。

当N_real取无穷大时，退化为理想不可压流体。

同理，推广到三维：梯度分量以理想梯度方向为峰值在球面上服从标准正态分布。把方程中的二维偏转改成三维偏转。

代码实现：

用噪声图采样来分配服从标准正态分布的随机数给每一单位质量流体，每单位质量流体N_real个随机数。噪声图生成方法：

#include <iostream>
#include <fstream>
#include <cstdlib>
#include <cmath>
#include <limits>
#include <string>
using namespace std;

#define PI 3.1415926
double rand_back(double i,double j)
{
    double u1=double(rand()%1000)/1000,u2=double(rand()%1000)/1000,r;
    
    static unsigned int seed=0;
    
    r=i+sqrt(j)*sqrt(-2.0*(log(u1)/log(exp(1.0))))*cos(2*PI*u2);
    
    return r;
}


void WriteNoise(const string fileName)
{
    ofstream ouTT;   //流名称
    
    const string folder = "/书/Ray Tracing in weekend/";
    const string name = folder + fileName;
    
    
    ouTT.open(name);  //流写入位置

    
    int nx = 512;
    int ny = 512;
    ouTT << "P3\n" << nx << " " << ny << "\n255\n";
    
    
    for (int j = ny-1; j >= 0; j--) {
        for (int i = 0; i < nx; i++) {
            
            
            int r = (rand_back(0, 1)+3) / 6 * 255;
            while (r < 0 || r > 255)    r = (rand_back(0, 1)+3) / 6 * 255;
            
            int g = (rand_back(0, 1)+3) / 6 * 255;
            while (g < 0 || g > 255)    g = (rand_back(0, 1)+3) / 6 * 255;
            
            int b = (rand_back(0, 1)+3) / 6 * 255;
            while (b < 0 || b > 255)    b = (rand_back(0, 1)+3) / 6 * 255;
            
            
            ouTT << r << " " << g << " " << b << "\n";    //流输出内容，用法和COUT一样
        }
    }
    
    ouTT.close();
}



int main(int argc, const char * argv[]) {

    WriteNoise("noise1.ppm");
    WriteNoise("noise2.ppm");
    WriteNoise("noise3.ppm");
    WriteNoise("noise4.ppm");
    
}

不可压流体直接用的Scrawk大佬的代码（https://github.com/Scrawk/GPU-GEMS-2D-Fluid-Simulation）在SubtractGradient.shader中把压力场梯度代入我的方程转换：

Shader "FluidSim/SubtractGradient" 
{   
    //=========================
    Properties 
    {
        _Noise1 ("Noise1 (RGB)", 2D) = "white" {}
        _Noise2 ("Noise2 (RGB)", 2D) = "white" {}
        _Noise3 ("Noise3 (RGB)", 2D) = "white" {}
        _Noise4 ("Noise4 (RGB)", 2D) = "white" {}
    }
    //=========================


	SubShader 
	{

    	Pass 
    	{
			ZTest Always

			CGPROGRAM
			#include "UnityCG.cginc"
			#pragma target 3.0
			#pragma vertex vert
			#pragma fragment frag

            //=========================
            #define PI 3.1415926

            
            sampler2D _Noise1;
            sampler2D _Noise2;
            sampler2D _Noise3;
            sampler2D _Noise4;
            //=========================

			
			uniform sampler2D _Velocity;
			uniform sampler2D _Pressure;
			uniform sampler2D _Obstacles;
			uniform float _GradientScale;
			uniform float2 _InverseSize;
			
			struct v2f 
			{
    			float4  pos : SV_POSITION;
    			float2  uv : TEXCOORD0;
			};

			v2f vert(appdata_base v)
			{
    			v2f OUT;
    			OUT.pos = UnityObjectToClipPos(v.vertex);
    			OUT.uv = v.texcoord.xy;
    			return OUT;
			}
			


            //=======================================
            float2 randomRotate(float ran, float2 g)
            {
                float theta = ran * 2 * PI - PI; //旋转角度
                float2 gg = float2(g.x * cos(theta) - g.y * sin(theta), g.x * sin(theta) + g.y * cos(theta));  //旋转

                return gg;  //旋转后速度方向
            }
            //=======================================



			float4 frag(v2f IN) : COLOR
			{
			
			    // Find neighboring pressure:
			    float pN = tex2D(_Pressure, IN.uv + float2(0, _InverseSize.y)).x;
			    float pS = tex2D(_Pressure, IN.uv + float2(0, -_InverseSize.y)).x;
			    float pE = tex2D(_Pressure, IN.uv + float2(_InverseSize.x, 0)).x;
			    float pW = tex2D(_Pressure, IN.uv + float2(-_InverseSize.x, 0)).x;
			    float pC = tex2D(_Pressure, IN.uv).x;
			
			    // Find neighboring obstacles:
			    float bN = tex2D(_Obstacles, IN.uv + float2(0, _InverseSize.y)).x;
			    float bS = tex2D(_Obstacles, IN.uv + float2(0, -_InverseSize.y)).x;
			    float bE = tex2D(_Obstacles, IN.uv + float2(_InverseSize.x, 0)).x;
			    float bW = tex2D(_Obstacles, IN.uv + float2(-_InverseSize.x, 0)).x;
			
			    // Use center pressure for solid cells:
			    if(bN > 0.0) pN = pC;
			    if(bS > 0.0) pS = pC;
			    if(bE > 0.0) pE = pC;
			    if(bW > 0.0) pW = pC;
			
			    // Enforce the free-slip boundary condition:
			    float2 oldV = tex2D(_Velocity, IN.uv).xy;
			    float2 grad = float2(pE - pW, pN - pS) * _GradientScale;
			   


                //===========================================================================
                //改压力场梯度

                float2 g = grad;

                //4 * 3 = 12 个服从标准正态分布的随机数
                float3 r1 = tex2D(_Noise1, IN.uv).rgb;
                float3 r2 = tex2D(_Noise2, IN.uv).rgb;  
                float3 r3 = tex2D(_Noise3, IN.uv).rgb;
                float3 r4 = tex2D(_Noise4, IN.uv).rgb;

                
                // 12 个旋转后的向量
                float2 u1r = randomRotate(r1.r, g);
                float2 u1g = randomRotate(r1.g, g);
                float2 u1b = randomRotate(r1.b, g);
                
                float2 u2r = randomRotate(r2.r, g);
                float2 u2g = randomRotate(r2.g, g);
                float2 u2b = randomRotate(r2.b, g);
                
                float2 u3r = randomRotate(r3.r, g);
                float2 u3g = randomRotate(r3.g, g);
                float2 u3b = randomRotate(r3.b, g);
                
                float2 u4r = randomRotate(r4.r, g);
                float2 u4g = randomRotate(r4.g, g);
                float2 u4b = randomRotate(r4.b, g);


                //N_real取12时：
                //----------------------
                //float2 g_real = u1r + u1g + u1b + u2r + u2g + u2b + u3r + u3g + u3b + u4r + u4g + u4b;
                //g_real = g_real * 0.123;
                //----------------------

                //N_real取3时：
                //----------------------
                float2 g_real = u1r + u1g + u1b;
                g_real = g_real * 0.49;//文章第二个视频
                //g_real = g_real * 0.5;//文章第一个视频
                //----------------------


                grad = g_real;


                //===========================================================================


                
                float2 newV = oldV - grad;
			    
			    return float4(newV,0,1);  
			}
			
			ENDCG

    	}
	}
}

Unity完整代码下载：

https://github.com/GORK44/RoughFluid/blob/master/RoughFluid.unitypackagegithub.com

2D粗粒流体模拟

2D粗粒流体模拟

来讲我做了什么。

代码实现：

Recommend

数值计算中的分析方法

不要企图用机器学习的视角来碰瓷好莱坞的特效仿真工业体系：对一篇标题党式报道的吐槽

我们需要什么样的机器学习仿真算法？

大规模线性方程组求解简史（上）

什么是好的技术文章

关于 Todo 工具

MacOS 启动失败和恢复记录

读《学会学习》

20171118 聊下最近

2017 总结

About Joyk