模擬流體粒子運動


源碼:http://files.cnblogs.com/flash3d/fluid.rar

拖動框框,可以讓流體出現慣性現象(是這說法么?+_+)。

今兒在一個小日本網站上看到人家用as模擬流體運動,炫目無比,於是模仿做了個類似的效果。

下面主要講其中涉及到的數學物理知識。

關於流體運動,上網查了下,撈到一些公式。由於用用在程序用模擬粒子運動,而不是進行精確的科學計算,所以,本人特意將公式化為最簡,將一些能默認的系數默認,能忽悠的函數忽悠,得出效率和效果比較平衡的算法(其中兩個光滑核函數參考小日本的)。

 

所謂模擬流體運動,並不是說每個粒子代表一個分子,而是比如你將一定量的羽毛放在風中,讓其飛揚,羽毛的運動就能近似模擬風的運動。在下在構思的時候,曾試圖通過范德華力這些分子層面的東西來模擬,走了一些彎路。

作為一個放入風中的羽毛,它會受到什么力呢?

羽毛會受到重力,來自其他羽毛帶來的壓力(並非直接接觸,而是通過驅動風),風的阻力。

將這些力合成,就是羽毛所受合力。

知道所受合力,根據牛頓定律,就能求出速度和位移。

問題I在於如何求這三個力?

重力好求的,壓力和阻力怎么求呢?

網上找了一些公式。

一個羽毛受到所有影響范圍內的羽毛的作用,我們認為被作用的那個羽毛,它的各個屬性被影響后的值,一般遵循以下規律(科學家研究的):

假設被作用的那個羽毛叫F,所有影響那個羽毛分別為F0 F1 F2...FN。

那么羽毛FJ對羽毛F的A屬性的影響是這么計算的:F.AJ=FJ.A*FJ.M/FJ.RO*W(DIST,RANGE)

其中M屬性為質量,RO屬性為密度,DIST為兩個羽毛之間的距離,RANGE為影響范圍的半徑,W函數為光滑核函數。

光滑核函數是一類滿足特定性質的函數,其通式為W(H,R),其函數值隨H的絕對值增大而衰減,直到H的絕對值大於等於R的時候,函數值就為0啦。光滑核函數還有一些性質,與我們的程序無關,比如他是一個偶函數,他的積分為1。要注意下,這個通式,他的自變量只有H一個,R是一個常數。

光滑核函數的意義在於,它可以描述一個粒子對周圍影響的強烈程度和周圍事物與該粒子的距離的關系。

於是,F.A=F.A0+F.A1+F.A2...+F.AJ+...+F.AN

首先計算粒子所受壓力F.FP。

一個粒子所受壓力是周圍粒子固有的壓力的累加。那么粒子固有壓力是什么呢?就是壓強啦。

高中學過氣體壓強公式,P=K(RO-RO'),RO-RO'是內外密度差,K是某熱力學常數,P就是壓強。

所以,要求所受壓力,還必須求出每個粒子所在位置的密度。動用剛才那個屬性疊加公式。

F.ROJ=FJ.RO*FJ.M/FJ.RO*W(DIST,RANGE)

密度能約掉,得到F.ROJ=FJ.M*W(DIST,RANGE)

其中,每個粒子的質量默認為1,於是F.ROJ=W(DIST,RANGE)

這里,我再講一下光滑核函數,其是一類函數,滿足其給定性質的均被稱為光滑核函數。所以每個屬性的影響疊加所用的光滑核函數也是不一樣的,具體取哪個函數,是要看要疊加的屬性的性質,或者看你需要的效果了。

這里我參考了小日本程序里面的光滑核函數:W(DIST,RANGE)=(1-DIST/RANGE)^2 自變量范圍0到RANGE

這樣 我們就能講兩個粒子的距離求出來,從而算出粒子間相互影響的量,並累加到粒子的密度中去。

代碼實現

for (i = 0; i < numParticles; i++)

{

pi = particles[i];//取外循環要操作的粒子

pi.numNeighbors = 0;//將粒子鄰居至空

pi.density = 0;//將粒子密度屬性初始化

/*內循環,用握手規則減少循環次數*/

for (j = 0; j < i; j++)

{

pj = particles[j];//取內循環粒子

dist = Math.pow(pi.x - pj.x, 2) + Math.pow(pi.y - pj.y, 2);//計算兩個粒子的距離平方

/*如果兩個粒子距離在范圍內進入判斷*/

if (dist < RANGE2)

{

dist = Math.sqrt(dist);//算出距離

weight = Math.pow(1 - dist / RANGE, 2);//算出光滑核函數的值

//將值累加到兩個粒子的密度里

pi.density += weight;

pj.density += weight;

//互加為鄰居

pi.neighbors[pi.numNeighbors++] = pj;

pj.neighbors[pj.numNeighbors++] = pi;

}

}

有了密度,我們就能算壓強啦,P=K(RO-RO')

RO'是流體在靜態的時候全局密度,在代碼中對應為DENSITY

代碼實現

/*算出粒子的壓力*/

for (i = 0; i < numParticles; i++)

{

pi = particles[i];

if (pi.density < DENSITY) pi.density = DENSITY;

pi.pressure = pi.density - DENSITY;//粒子壓力等於K(粒子密度-靜態密度),K默認為1

}

有了壓強,終於能計算粒子對其他粒子的壓力了。

這里聰明的孩子馬上意識到,壓強是標量,壓力是矢量,標量累加怎么就成了矢量了呢?

公式上的做法是將壓強做一介微分和,也就是哈密頓算子(這個可以去百度),不過,我們程序上的做法是先算標量,再處理矢量。因為力的方向,其實我們早就知道了,從發力粒子指向受力粒子。

所以得到的公式就是F.PJ=FJ.P*FJ.M/FJ.RO*W(DIST,RANGE)

不過不幸的是,這個公式是“不平衡”的,也就是說,位於不同壓強區的兩個粒子之間的作用力不等,所以計算中一般使用雙方粒子壓強的算術平均值代替單個粒子的壓力

F.PJ=(FJ.P+F.P)/2*FJ.M/FJ.RO*W(DIST,RANGE)

其中光滑核函數W(DIST,RANGE)=1 - DIST / RANGE

質量為1

F.PJ=(FJ.P+F.P)/2/FJ.RO*(1 - DIST / RANGE)

求出作用力之后,將其分成x和y兩個分量分別累加

其中還要乘上壓力系數也就是表面積。

F.PJ=F.PJ*PRESSURE

代碼實現

dx = pi.x - pj.x;

dy = pi.y - pj.y;

dist = Math.sqrt(Math.pow(dx, 2) + Math.pow(dy, 2));//計算兩個粒子的距離

weight = 1 - dist / RANGE;//光滑核函數

pressure = weight * (pi.pressure + pj.pressure) / (2 * pj.density) * PRESSURE;//計算粒子對粒子的壓力標量

/*將壓力分成向量的兩個值累計到受力中*/

dist = 1 / dist;

dx *= dist;

dy *= dist;

pi.fx += dx * pressure;

pi.fy += dy * pressure

然后是計算阻力,阻力是粘度系數乘上運動速度差。所以,我們對運動速度差進行屬性疊加。

F.VJ=FJ.V*FJ.M/FJ.RO*W(DIST,RANGE)

F.VJ=(F.v-FJ.v)*FJ.M/FJ.RO*W(DIST,RANGE)

其中光滑核函數W(DIST,RANGE)=1 - DIST / RANGE

於是化簡得到F.VJ=(F.v-FJ.v)/FJ.RO*(1 - DIST / RANGE)

最后得到阻力F.VJ=F.VJ*VISCOSITY

在代碼中,由於速度差是矢量,所以,乘上速度的時候,分別在x和y兩個方向處理。

代碼實現

viscosity = weight / pj.density * VISCOSITY;//計算除速度外的阻力系數

//計算粒子相對速度的向量

dx = pi.vx - pj.vx;

dy = pi.vy - pj.vy;

//乘入速度因子得出粒子對粒子的壓力累加入受力

pi.fx -= dx * viscosity;

pi.fy -= dy * viscosity

最后剩下其他外力作用了

其中有重力,瓶子被移動引擎的加速度,還有碰到牆壁反彈的力

重力直接累加到受力內

瓶子的加速度反向累加到受力內

牆壁的彈力,這里將牆壁想象成彈簧,彈性系數K,那么如果粒子剛接觸到牆壁被視為開始接觸彈簧,那么粒子超出牆壁的范圍,被視為彈簧的壓縮量L,那么產生的力F=-K*L,這個力是被即時累加到速度里面的。

重力和瓶子加速度累加的代碼

for (i = 0; i < numParticles; i++)

{

pi = particles[i];

//考慮外部瓶子受力,將其反向累計

pi.fx -= bottle_fx;

pi.fy -= bottle_fy;

pi.fy += Main.GRAVITY;

pi.move();//將受力作用成速度和位移

}

撞牆后反饋力的代碼(彈性系數為1)

//如果撞牆,速度為0,如果穿過牆,則受到和超過部分相同的力的反饋

if (x < 4) vx += 4 - x;

if (y < 4) vy += 4 - y;

if (x > 461) vx += 461 - x;

if (y > 461) vy += 461 - y;

至於其他部分,如力怎么累加到速度,速度怎么累加到位移之類的問題,我就先不講了。祝大家周末愉快!!

感謝四樓提醒,本文部分公式參考自http://www.thecodeway.com/blog/?p=1557


免責聲明!

本站轉載的文章為個人學習借鑒使用,本站對版權不負任何法律責任。如果侵犯了您的隱私權益,請聯系本站郵箱yoyou2525@163.com刪除。



 
粵ICP備18138465號   © 2018-2025 CODEPRJ.COM