電影工業的流體模擬（特別篇）---渦方法（二），N體問題初探以及最簡單的流體模擬程序

01-30

本文作者@張心欣，個人主頁Xinxin Zhang UBC

本文著作權歸原作者張心欣所有，未經原作者許可的轉載均屬侵權行為，作者保留對侵權行為追究責任的權利。

本文已取得原作者授權發布。

在渦方法（一）中，我們初步介紹了渦方法的數學背景，讓我們來簡單複習一下最重要的幾個式子：

今天，我們要看一下求解這個方程組的一種方法，在進行實際的演算法求解前，我們先跳過一些數學推導並了解一強大的數學工具，格林函數（Greens function），有關數學推導在鏈接中，希望感興趣的讀者自行參閱。簡而言之，通過格林函數這個數學工具，我們有如下事實：

如果在一個無邊界的空間中有如下泊松方程，

則有， 2維情況下

3維情況下

於是在一個無邊界的空間里，格林函數恰恰給我們的流體速度提供了一個半解析的解, 如果我們將上面的兩個積分式直接離散化，並在求和符號的前邊對 $x$ 求 $nablatimes$ 運算元，就能分別得到在2維和3維情況下速度的半解析解，又稱為Biot-Savart法則：

2D: $u_{i}=sum_{j,jneq i}^{N}frac{(0,0,omega_{j})times(x_{i}-x_{j})}{2pi r_{ij}^2}$

3D： $u_{i}=sum_{j,jneq i}^{N}frac{boldsymbol{omega_j}times(x_{i}-x_{j})}{4pi r_{ij}^3}$

其中 $r_{i,j}$ 是第i個渦粒子和第j個渦粒子間的距離。在此我們注意到，實則每一個粒子的運動速度，是需要對空間中所有的其它粒子的渦量求和得到，這個問題和求解宇宙中的萬有引力有相同的形式，同樣都是N-body problem（N體問題）。

在實際工程中，我們通常不使用如上最簡單的渦粒子，而是使用一種更為穩定的表達，渦泡兒（vortex blob），其詳細的推導可在以上鏈接中找打，從實用性的角度出發來講，我們意識到當兩個渦粒子的位置無限接近時， $r_{i,j}$ 趨於0，使得上面的式子出現了奇異點，難以計算。vortex blob的推導則正好是對這兩個式子加了一些奇異點的修正（mollify），以二維為例，有

$u_{i}=sum_{j,jneq i}^{N}frac{(0,0,omega)times(x_{i}-x_{j})}{2pi r_{ij}^2}left( 1-e^{r_{ij}^2/epsilon^2} right)$

至此為止，我們其實已經可以構造一個極簡單的程序來模擬一種用普通的NS方程求解難以模擬的現象 -- vortex ring leapfrogging

為了對這種現象進行一些定性的分析，我們寫一個2D的模擬代碼來看看，從vortex的結構開始：

struct vortex2Dn{ntdouble x; double y; double vort;ntvortex2D(){}ntvortex2D(double _x,double _y, double _vort)nt{nttx=_x; y=_y; vort=_vort;nt}ntvortex2D(vortex2D &vortex)nt{nttx = vortex.x; y=vortex.y; vort=vortex.vort;nt}n};n

我們希望模擬與視頻中的實驗類似的情況，不妨在一開始設定兩對vortex pair

std::vector<vortex2D> vortex_particles;ntvortex_particles.resize(0);ntvortex_particles.push_back(vortex2D(0.0, 1.0, 1.0));ntvortex_particles.push_back(vortex2D(0.0, -1.0, -1.0));ntvortex_particles.push_back(vortex2D(0.0, 0.3, 1.0));ntvortex_particles.push_back(vortex2D(0.0, -0.3, -1.0));n

如此的設定對應的初始畫面（T=0.1時）是：

接下來，我們根據之前討論的二維速度計算公式來計算某一個坐標位置在所有渦粒子影響下的速度：

//u component of velocity from a single vortexndouble compute_u_from_single_vortex(double x, double y, vortex2D &vortex)n{ntdouble r_ij2 = (x-vortex.x)*(x-vortex.x) + (y-vortex.y)*(y-vortex.y);ntreturn vortex.vort*(vortex.y-y)/(r_ij2*3.1415926) * 0.5 *(1.0 - exp(-r_ij2/(EPS*EPS)));n}nn//v component of velocity from a single vortexndouble compute_v_from_single_vortex(double x, double y, vortex2D &vortex)n{ntdouble r_ij2 = (x-vortex.x)*(x-vortex.x) + (y-vortex.y)*(y-vortex.y);ntreturn vortex.vort*(x-vortex.x)/(r_ij2*3.1415926) * 0.5 *(1.0 - exp(-r_ij2/(EPS*EPS)));n}nn//u component of velocity from all vortex influencendouble compute_u_full_influence(double x, double y, std::vector<vortex2D> &vortex)n{ntdouble u=0;ntfor (int i=0;i<vortex.size();i++)nt{nttu+=compute_u_from_single_vortex(x,y,vortex[i]);nt}ntreturn u;n}nn//v component of velocity from all vortex influencendouble compute_v_full_influence(double x, double y, std::vector<vortex2D> &vortex)n{ntdouble v=0;ntfor (int i=0;i<vortex.size();i++)nt{nttv+=compute_v_from_single_vortex(x,y,vortex[i]);nt}ntreturn v;n}n

至此為止，我們只需要再有幾個for循環，就可以運行這個模擬了！！！！！在我們開始模擬之前，我們先初始化一些tracer粒子可以用來可視化我們的流體運動：

int num_tracer = 200000;ntdouble *pos_x = new double[num_tracer];ntdouble *pos_y = new double[num_tracer];ntntint num=0;ntwhile(num<num_tracer)nt{nttdouble x=frand(-0.5,0.5);nttdouble y=frand(-1.5,1.5);nttpos_x[num] = x;nttpos_y[num] = y;nttnum++;nt}n

接下來，我們可以開始我們的模擬了，基本上，我們只需要設一個for循環做兩件事：

double dt = 0.1;nfor (int T=0;T<300;T++)n{n //integrate tracersn //integrate vortex particlesn}n

在我們的程序中，對於給定的渦粒子分布，我們將採用RK3格式來進行積分tracer粒子：

#pragma omp parallel fornfor (int i=0;i<num_tracer;i++)n{n rk3_integrate_pos(pos_x[i],pos_y[i],vortex_particles, dt);n}n

至於積分渦粒子的運動，我們僅採用簡單的前向歐拉，詳情請見最後附上的代碼。

把今天的代碼編譯運行後並採用適當的程序可視化運行結果，我們能夠得到：

在這一節中，請注意到我們所運行的演算法之簡潔以及它所產生的結果之複雜，之具有挑戰性，我們的數值演算法在數值層面上完全的保證環量守恆，這是使用傳統的NS求解根本不可能達到的！！如果採用傳統的NS方程求解，像這樣的leapfrogging現象，根本不可能被在數值解中重現出來。（作者幾乎沒有見過通過求解NS方程完成兩個個以上leapfrogging循環的結果。）

以下是本章的附帶常式全篇 main.cpp

#include <stdio.h>n#include <random>n#include <math.h>nn//define the mollify radiusn#define EPS 0.01nn//return a random number in [a,b]ndouble frand(double a, double b)n{ntdouble r = (double)rand()/(double)RAND_MAX;ntreturn (b-a)*r+a;n}nn//a simple 2D vortexc structnstruct vortex2Dn{ntdouble x; double y; double vort;ntvortex2D(){}ntvortex2D(double _x,double _y, double _vort)nt{nttx=_x; y=_y; vort=_vort;nt}ntvortex2D(vortex2D &vortex)nt{nttx = vortex.x; y=vortex.y; vort=vortex.vort;nt}n};nn//u component of velocity from a single vortexndouble compute_u_from_single_vortex(double x, double y, vortex2D &vortex)n{ntdouble r_ij2 = (x-vortex.x)*(x-vortex.x) + (y-vortex.y)*(y-vortex.y);ntreturn vortex.vort*(vortex.y-y)/(r_ij2*3.1415926) * 0.5 *(1.0 - exp(-r_ij2/(EPS*EPS)));n}nn//v component of velocity from a single vortexndouble compute_v_from_single_vortex(double x, double y, vortex2D &vortex)n{ntdouble r_ij2 = (x-vortex.x)*(x-vortex.x) + (y-vortex.y)*(y-vortex.y);ntreturn vortex.vort*(x-vortex.x)/(r_ij2*3.1415926) * 0.5 *(1.0 - exp(-r_ij2/(EPS*EPS)));n}nn//u component of velocity from all vortex influencendouble compute_u_full_influence(double x, double y, std::vector<vortex2D> &vortex)n{ntdouble u=0;ntfor (int i=0;i<vortex.size();i++)nt{nttu+=compute_u_from_single_vortex(x,y,vortex[i]);nt}ntreturn u;n}nn//v component of velocity from all vortex influencendouble compute_v_full_influence(double x, double y, std::vector<vortex2D> &vortex)n{ntdouble v=0;ntfor (int i=0;i<vortex.size();i++)nt{nttv+=compute_v_from_single_vortex(x,y,vortex[i]);nt}ntreturn v;n}nn//a simple rk3 integratornvoid rk3_integrate_pos(double &x, double &y, std::vector<vortex2D> &vortex, double dt)n{ntdouble x0 = x;double y0 = y;ntdouble u0 = compute_u_full_influence(x0,y0,vortex);ntdouble v0 = compute_v_full_influence(x0,y0,vortex);ntdouble x1 = x0 + 0.5*dt*u0; double y1 = y0 + 0.5*dt*v0;nntdouble u1 = compute_u_full_influence(x1,y1,vortex);ntdouble v1 = compute_v_full_influence(x1,y1,vortex);ntdouble x2 = x0 + 0.75*dt*u1; double y2 = y0 + 0.75*dt*v1;nntdouble u2 = compute_u_full_influence(x2,y2,vortex);ntdouble v2 = compute_v_full_influence(x2,y2,vortex);ntx += 0.222222222222*dt*u0 + 0.333333333333*dt*u1 + 0.444444444444*dt*u2;nty += 0.222222222222*dt*v0 + 0.333333333333*dt*v1 + 0.444444444444*dt*v2;n}n//write tracer particle data to a file so that can be visulized by matlabnvoid write_file(double *pos_x, double *pos_y, int num, int frame)n{nt//I have a openGL executable to quickly draw 3D particles, so all 2D data is transfered to 3D for rendering.ntchar filename[256];ntsprintf(filename,"Particle_data%04d.bin",frame);ntfloat *data;ntdata = new float[num*4];ntfor (int i=0;i<num;i++)nt{nttdata[i*4+0] = pos_x[i];nttdata[i*4+1] = pos_y[i];nttdata[i*4+2] = 0;nttdata[i*4+3] = 1;nt}ntFILE *f;ntf = fopen(filename,"wb");ntfwrite(data,sizeof(float),4*num,f);nntdelete[]data;ntfclose(f);nn}nnn// our main function computes two vortex ring leapforggingnint main(int argc, char * argv[])n{nt//initialize computation : 2 vortex ring pair and a lot tracer particles used to visualizennt//2 vortex ring pairs;ntstd::vector<vortex2D> vortex_particles;ntvortex_particles.resize(0);ntvortex_particles.push_back(vortex2D(0.0, 1.0, 1.0));ntvortex_particles.push_back(vortex2D(0.0, -1.0, -1.0));ntvortex_particles.push_back(vortex2D(0.0, 0.3, 1.0));ntvortex_particles.push_back(vortex2D(0.0, -0.3, -1.0));nnnt//200000 tracer particles ntint num_tracer = 200000;ntdouble *pos_x = new double[num_tracer];ntdouble *pos_y = new double[num_tracer];ntntint num=0;ntwhile(num<num_tracer)nt{nttdouble x=frand(-0.5,0.5);nttdouble y=frand(-1.5,1.5);nttpos_x[num] = x;nttpos_y[num] = y;nttnum++;nt}nnt//our simulationntdouble dt = 0.1;ntfor (int T=0;T<300;T++)nt{nttwrite_file(pos_x,pos_y,num_tracer,T);nntt//few substeps, not necessarynttfor(int substep=0;substep<4;substep++)ntt{ntt//integrate tracers, we are going to use rk3 integratornttt#pragma omp parallel forntttfor (int i=0;i<num_tracer;i++)nttt{nttttrk3_integrate_pos(pos_x[i],pos_y[i],vortex_particles, dt);nttt}nttnttt//integrate vortex particlesntttstd::vector<vortex2D> vortex_particles_temp;ntttvortex_particles_temp.resize(vortex_particles.size());nttt//before we integrate, we should copy vortex particles tonttt//freeze them in space.ntttfor (int i=0;i<vortex_particles.size();i++)nttt{nttttvortex_particles_temp[i] = vortex_particles[i];nttt}ntttfor(int i=0;i<vortex_particles_temp.size();i++)nttt{ntttnttttdouble u=0;double v = 0;nttttfor (int j=0;j<vortex_particles.size();j++)ntttt{ntttttif(j!=i){nnttttttu += compute_u_from_single_vortex(vortex_particles_temp[i].x,vortex_particles_temp[i].y,vortex_particles[j]);nttttttv += compute_v_from_single_vortex(vortex_particles_temp[i].x,vortex_particles_temp[i].y,vortex_particles[j]);nttttt}ntttt}nttttvortex_particles_temp[i].x += dt*u;nttttvortex_particles_temp[i].y += dt*v;nttt}nttt//copy back datantttfor (int i=0;i<vortex_particles.size();i++)nttt{nttttvortex_particles[i] = vortex_particles_temp[i];nttt}nttntt}nttnttprintf("step %d donen",T);nt}ntntdelete[]pos_x; delete[]pos_y;ntreturn 0;n}n

請毫不猶豫地關注我們：

我們的網站：GraphiCon

知乎專欄：GraphiCon圖形控 - 知乎專欄

公眾號：GraphiCon

如果你有什麼想法，建議，或者想加入我們，你可以：

給我們發郵件：hi@graphicon.io

加入我們的QQ群：SIQGRAPH（342086343）

加入我們的slack群：Join us on Slack!

GraphiCon長期接受投稿，如果你想投稿給我們可以通過上面的方式聯繫我們！

我們是圖形學愛好者，歡迎加入我們！