渦の回転方向で色分け
自作の非圧縮粘性流体シミュレーションです｡
中国の湖みたいできれいですね｡
あんまり早い流速を与えると発散します｡
// forked from knd's Lake of china
package
{
import flash.display.Bitmap;
import flash.display.BitmapData;
import flash.display.Sprite;
import flash.events.Event;
import flash.events.MouseEvent;
/**
* 自作の非圧縮粘性流体シミュレーションです｡
* 中国の湖みたいできれいですね｡
* あんまり早い流速を与えると発散します｡
*/
[SWF(width="465", height="465", backgroundColor="0x0", frameRate="90")]
public class TestFluid extends Sprite
{
private var map:FluidMap;
private var dat:BitmapData;
private var bmp:Bitmap;
private var wid:uint = 31;
private var hei:uint = 31;
private var lastX:Number;
private var lastY:Number;
private var layer:Sprite;
private const SCALE:uint = 15;//セルの大きさ
private const RE:Number = 10;//レイノルズ数。流体の粘性に関わる
public function TestFluid() {
map = new FluidMap(wid, hei, RE);
dat = new BitmapData(wid, hei);
bmp = new Bitmap(dat);
layer = new Sprite();
bmp.scaleX = bmp.scaleY = SCALE;
addChild(bmp);
addChild(layer);
//map.next(200);
lastX = mouseX;
lastY = mouseY;
stage.addEventListener(MouseEvent.CLICK, reset);
mapColor();
addEventListener(Event.ENTER_FRAME, loop);
}
private const N :uint = 16;
private const V: Number = 0.01;//1フレームあたりのマウス移動量に対し与える流速の比
private function loop(e:Event):void {
var ix:Number = lastX / SCALE;
var iy:Number = lastY / SCALE;
var dx:Number = (mouseX - lastX) /SCALE / N;
var dy:Number = (mouseY - lastY) /SCALE / N;
if ( dx == 0 && dy == 0) {
map.next(50);//マウスの移動がなければGS法の計算量を多めにする
mapColor();
return;
}
ix = ix >= wid ? wid - 1: (ix < 0? 0: ix);
iy = iy >= hei ? hei - 1: (iy < 0? 0: iy);
for (var n:uint = 0; n < N; n++) {
map.setVx( (mouseX - lastX)*V, ix, iy);
map.setVy( (mouseY - lastY)*V, ix, iy);
ix += dx;
iy += dy;
}
ix += dx;
iy += dy;
lastX = mouseX;
lastY = mouseY;
map.next(20);
mapColor();
}
private function mapColor():void {
layer.graphics.clear();
layer.graphics.lineStyle(0, 0, 0.3);
var vort:Number;
dat.lock();
var vx:Number;
var vy:Number;
var b:Number = 0.1;//bias
var color:uint;
for (var ix:int = 0; ix < wid; ix++) {
for (var iy:int = 0; iy < hei; iy++) {
vx = map.getVx(ix, iy) ;
vy = map.getVy(ix, iy) ;
layer.graphics.moveTo((0.5 + ix) * SCALE, (0.5 + iy) * SCALE );
layer.graphics.lineTo((0.5 + ix) * SCALE + 50.0/b*vx, (0.5 + iy) * SCALE + 50.0/b*vy)
vx += b;
vy += b;
/*
color = (vx <= 0.0)? 0: ((vx >= 2*b) ? 0xff: uint(vx / b * 128.0));
color <<= 16;
color |= 0x8000;
color |= (vy <= 0.0)? 0: ((vy >= 2*b) ? 0xff: uint(vy / b * 128.0));
*/
/**/
vort = map.getVorticity(ix, iy);
color = vort >= 0.0?
(vort < b? (uint(vort/b * 255.0) << 16) ^ 0xffffff : 0x00ffff) :
(vort > -b? (uint(-vort/b * 255.0) ^ 0xffffff): 0xffff00);
/**/
dat.setPixel(ix, iy, color);
}
}
dat.unlock();
}
private function reset(e:MouseEvent):void {
removeEventListener(Event.ENTER_FRAME, loop);
map.reset();
map.next(20);
mapColor();
addEventListener(Event.ENTER_FRAME, loop);
}
}
}
class FluidMap
{
protected var width:uint;
protected var height:uint;
public var Re:uint;
protected var cells :Vector.<Vector.<FluidCell>>;
/**
* 流体マップを境界条件とサイズを引数にして新規作成する｡
* @param boundary
* @param width
* @param height
*/
public function FluidMap(width:uint, height:uint, r:Number = 1000) {
this.Re = r;
this.width = width;
this.height = height;
init();
}
/**
* FluidCellを生成する｡
*/
protected function init():void {
var x:uint;
var y:uint;
var c:FluidCell;
cells = new Vector.<Vector.<FluidCell>>;
for (x = 0; x < width+3; x++) {
cells[x] = new Vector.<FluidCell>;
for (y = 0; y < height + 3; y++) {
if (x < width && y < height) {
cells[x][y] = new FluidCell(Re);
} else {
cells[x][y] = new NoFluidCell(Re);//流速がゼロの境界条件セル
}
}
}
// [width+2], [height+2]は[-1]インデクスに対応しているという変なトリック
var xInc:uint;
var xDec:uint;
var yInc:uint;
var yDec:uint;
for (x = 0; x < width+3; x++) {
for (y = 0; y < height+3; y++) {
c = cells[x][y];
if (x == (width + 2)) xInc = 0;
else xInc = x + 1;
if (x == 0) xDec = width + 2;
else xDec = x - 1;
if (y == (height + 2)) yInc = 0;
else yInc = y + 1;
if (y == 0) yDec = height + 2;
else yDec = y - 1;
if(y != height + 2) c.setNorthCell(cells[x][yDec]);
if(x != width + 1) c.setEastCell(cells[xInc][y]);
if((x != width + 1)&&(y != height + 1))c.setSouthEastCell(cells[xInc][yInc]);
if(y != height + 1)c.setSouthCell(cells[x][yInc]);
if(x != width + 2)c.setWestCell(cells[xDec][y]);
if((x != width + 2)&&(y != height + 2))c.setNorthWestCell(cells[xDec][yDec]);
}
}
}
//セルの流速をゲット･セット｡
//引数のチェックは呼び出し側でやってください｡
public function setVx(vx:Number, x:uint, y:uint):void {
cells[x][y].u = vx;
}
public function setVy(vy:Number, x:uint, y:uint):void {
cells[x][y].v = vy;
}
public function getVx(x:uint, y:uint):Number {
return cells[x][y].u;
}
public function getVy(x:uint, y:uint):Number {
return cells[x][y].v;
}
public function setPressure(p:Number, x:uint, y:uint):void {
cells[x][y].p = p;
}
public function getPressure(x:uint, y:uint):Number {
return cells[x][y].p;
}
public function getVorticity(x:uint, y:uint):Number {
return cells[x][y].vort;
}
public function reset():void {
var cs:Vector.<FluidCell>;
var c:FluidCell;
for each(cs in cells) for each(c in cs) {
c.reset();
}
}
/**
* マップ上の時間を進めて状態を更新する
*/
public function next(repeat:uint=1):void {
var cs:Vector.<FluidCell>;
var c:FluidCell;
//GS法の計算回数の決定にエラーの大きさを利用してもいいかも
//var e:Number;
//var eMax:Number = 0.0;
for each(cs in cells) for each(c in cs) {
c.updateDiffV();
}
for each(cs in cells) for each(c in cs) {
c.updateDiffP();
}
while (repeat > 0) {
for each(cs in cells) for each(c in cs) {
c.stepGSM();
}
for each(cs in cells) for each(c in cs) {
c.updatePressure();
//e = c.getGSMError();
//if (e < 0) e = -e;
//if (eMax < e) eMax = e;;
}
//if (eMax < 0.01) break;
repeat--;
}
for each(cs in cells) for each(c in cs) {
c.solveNSE();
}
for each(cs in cells) for each(c in cs) {
c.updateVelocity();
}
}
}
/**
* 粘性･非圧縮流体のスタッガード格子
*/
class FluidCell
{
/*
*簡単な計算の流れについて
* 現在の流速から流速の､tについての偏微分を求める｡
* このとき一旦圧力の勾配は無視することにする｡
* 得られた流速は当然連続の式を満足しないので､
* 新しく圧力を設定することで補正する｡
*
* updateDiffV 流速まわりでの偏微分値を更新
* updateDiffP 圧力まわりでの偏微分値を更新
* stepGSM 適当な回数繰り返して圧力についてのPoisson方程式を解く｡
* getGMSError 現在の解の精度を得る｡
* updatePressure 圧力の値を確定
* solveNSE 求めた圧力を用いてNavier-Stokes方程式から流速のtについての偏微分を求める
* updateVelocity 現在の流速に求めた偏微分値を足す｡
*
* 以上のながれにて｡
* */
/*
* スタッガード格子なるものを使用するので
* 流速と圧力で基準とする位置が異なる｡
* セルの真ん中に流速を配置､
* その左上(NorthWest)に圧力を配置した｡
* */
//以下は流速の位置を基準にとる
public var u:Number ; //ｘ方向の流速
public var v:Number ; //y方向の流速
public var px:Number; //pのx偏微分
public var py:Number; //pのy偏微分
public var vort:Number; //渦度 rot v
private var uD:Number; // - u(du/dx) - v(du/dy) + Δu/Re
private var vD:Number; // - u(dv/dx) - v(dv/dy) + Δv/Re
//以下は圧力の位置を基準にとる
private var ux:Number; //uのｘ偏微分
private var uy:Number; //uのy偏微分
private var vx:Number; //ｖのx編微分
private var vy:Number; //vのy編微分
private var uA:Number; //uの平均
private var vA:Number; //vの平均
private var ut:Number; //
private var vt:Number; //
private var pPsn:Number; //圧力についてのポアソン方程式の値 (-Δp)
protected var pNex:Number ; //圧力の値を一時的に格納
public var p:Number ; //確定済みの圧力
private var RI:Number;//レイノルズ数の逆数で計算します｡
/**
* とりあえず､コンストラクタはレイノルズ数のみを指定する｡
* @param Re
*/
public function FluidCell(Re:Number) {
reset();
RI = 1.0 / Re;
}
//隣接するセルに関してすべてnullのままにしないこと
protected var nc: FluidCell;
protected var sc: FluidCell;
protected var wc: FluidCell;
protected var ec: FluidCell;
protected var nwc:FluidCell;
protected var sec:FluidCell;
public function setNorthCell(c:FluidCell):FluidCell {
nc = c;
return this;
}
public function setSouthCell(c:FluidCell):FluidCell {
sc = c;
return this;
}
public function setWestCell(c:FluidCell):FluidCell {
wc = c;
return this;
}
public function setEastCell(c:FluidCell):FluidCell {
ec = c;
return this;
}
public function setNorthWestCell(c:FluidCell):FluidCell {
nwc = c;
return this;
}
public function setSouthEastCell(c:FluidCell):FluidCell {
sec = c;
return this;
}
/**
* 流速座標における偏微分の値を更新
*/
public function updateDiffV():void {
uD = -(u * (ec.u - wc.u) + v * (sc.u - nc.u) ) * 0.5 + RI * (nc.u + sc.u + wc.u + ec.u - 4.0 * u);
vD = -(u * (ec.v - wc.v) + v * (sc.v - nc.v) ) * 0.5 + RI * (nc.v + sc.v + wc.v + ec.v - 4.0 * v); }
/**
* 圧力座標における偏微分の値を更新
*/
public function updateDiffP():void {
ux = (this.u - wc.u + nc.u - nwc.u) * 0.5;
uy = (this.u + wc.u - nc.u - nwc.u) * 0.5;
vx = (this.v - wc.v + nc.v - nwc.v) * 0.5;
vy = (this.v + wc.v - nc.v - nwc.v) * 0.5;
// (du/dx+dv/dy) + d/dt(du/dx+dv/dy) = 0 の関係を使ってNavier-Stokes方程式を書き換える
pPsn = -ux - vy - 0.5 * ( this.uD - wc.uD + nc.uD - nwc.uD + this.vD + wc.vD - nc.vD - nwc.vD );
vort = vx - uy;
}
//SOR法に使用する定数
private static const FOR :Number = 0.95;//Factor of relaxation
/**
* Gauss-Seidel method + SOR method 1ステップ
*/
public function stepGSM():void {
pNex += FOR * (0.25 * (pPsn + nc.p + sc.p + wc.p + ec.p) -p);
}
public function updatePressure():void {
p = pNex;
}
/**
* GS法による現在の値の誤差を-Δpを計算して求める｡
* @return
*/
public function getGSMError():Number {
return (4.0 * pNex - nc.pNex - sc.pNex - wc.pNex - ec.pNex) - pPsn;
}
/**
* 圧力の勾配を再計算
* Solve the Navier-Stokes Equations
*/
public function solveNSE():void {
px = (ec.p - this.p + sec.p - sc.p) * 0.5;
py = (sc.p - this.p + sec.p - ec.p) * 0.5;
ut = -px + uD;
vt = -py + vD;
}
/**
*
*/
public function updateVelocity():void {
u += ut;
v += vt;
}
public function reset():void {
u = 0.0;
v = 0.0;
ux = 0.0;
uy = 0.0;
vx = 0.0;
vy = 0.0;
uD = 0.0;
vD = 0.0;
p = 0.0;
pNex = 0.0;
pPsn = 0.0;
px = 0.0;
py = 0.0;
ut = 0.0;
vt = 0.0;
vort = 0.0;
}
}
/**
* 流速ゼロのスタッガード格子
*/
class NoFluidCell extends FluidCell
{
public function NoFluidCell(Re:Number) {
super(Re);
}
//基本的に流速を更新するような操作は一切行わない｡
override public function updateDiffV():void
{
if ((nc != null) && (wc != null) && (sc != null) && (ec != null)) {
super.updateDiffV();
}
}
override public function updateDiffP():void {
if ((nc != null) && (wc != null) && (nwc != null)) {
super.updateDiffP();
}
}
override public function stepGSM():void
{
//上下左右にセルがない場合には圧力についての境界条件を設定する
if (nc == null) this.pNex = sc.p;
else if (wc == null) this.pNex = ec.p;
else if (sc == null) this.pNex = nc.p;
else if (ec == null) this.pNex = wc.p;
//上下左右にセルがある場合にはポアソン方程式を解く
else super.stepGSM();
}
override public function getGSMError():Number
{
if ((nc != null) && (wc != null) && (sc != null) && (ec!=null)) {
return super.getGSMError();
} else {
return 0.0;
}
}
override public function solveNSE():void { }
override public function updateVelocity():void { }
}