p5.jsで流体シミュレーション風

この記事はProcessing Advent Calendar 2018 20日目の記事です。

今回は2次元流体シミュレーションっぽいものを作ってみました。

作ったもの

マウスでドラッグすると流れを作ることができ、キーボードを押すとリセットできます。

解説

初期化

まず、シミュレーションの空間を格子に区切り、各格子点にそれぞれ流体の速度と色を持たせます。

// variables for simulations
let colors;
let velocities;
...
function initialize() {
	colors = new Array(GRID_NUM + 1);
	velocities = new Array(GRID_NUM + 1);
	const colorFunc = getInitialColorFunction();
	for(let xi = 0; xi < GRID_NUM + 1; xi++) {
		colors[xi] = new Array(GRID_NUM + 1);
		velocities[xi] = new Array(GRID_NUM + 1);
		for (let yi = 0; yi < GRID_NUM + 1; yi++) {
			colors[xi][yi] = colorFunc(xi, yi);
			velocities[xi][yi] = createVector(0.0, 0.0);
		}
	}
	...
}

イメージにすると、速度と色それそれ以下のようになります。

更新

シミュレーションの1ステップごとに格子点の速度と色を1ステップ前の値をもとに更新していきます。

function update() {
	const now = millis() * 0.001;
	const deltaTime = (now - lastSeconds) * TIME_SCALE;
	lastSeconds = now;
	updateVelocities(deltaTime);
	updateColors(deltaTime);
}

速度の更新

速度の更新では、外力の追加、移流、減衰を順番に行っていきます。

function updateVelocities(deltaTime) {
	let nextVels = addForceToVelocities(velocities, deltaTime);
	nextVels = advectVelocities(nextVels, deltaTime);
	nextVels = attenuateVelocities(nextVels, deltaTime);	
	velocities = nextVels;
}

まず、外力を加えます。今回はインタラクションを行うためにマウスの位置と動きに応じて速度を加えるようにしています。

function addForceToVelocities(vels, deltaTime) {
	if (!mouseIsPressed) {
		return vels;
	}
	const mousePos = createVector(mouseX, mouseY).div(width).mult(GRID_NUM).mult(GRID_SIZE);
	const prevMousePos = createVector(pmouseX, pmouseY).div(width).mult(GRID_NUM).mult(GRID_SIZE);
	const force = p5.Vector.sub(mousePos, prevMousePos).mult(FORCE_SCALE).limit(MAX_FORCE).mult(deltaTime);

	const nextVels = new Array(GRID_NUM + 1);
	for (let xi = 0; xi < GRID_NUM + 1; xi++) {
		nextVels[xi] = new Array(GRID_NUM + 1);
		for (let yi = 0; yi < GRID_NUM + 1; yi++) {
			const pos = createVector(xi, yi).mult(GRID_SIZE);
			const d = p5.Vector.mag(p5.Vector.sub(pos, mousePos));
			const forceScale = max((FORCE_RANGE - d) / FORCE_RANGE, 0.0);
			nextVels[xi][yi] = p5.Vector.add(vels[xi][yi], force.copy().mult(forceScale));
		}
	}
	return nextVels;
}

次に移流という流体を流す処理を行います。時刻$t$にある格子点$p(t)$上に存在する流体はその格子点上の速度を$v(t)$とすると、$\Delta t$秒が経過後に$p(t) + \Delta tv(t)$に移動します。
ある格子点にある流体が$\Delta t$秒後にどこに移動するかはわかりましたが、必要なのは$\Delta t$秒後にある格子点に流れてくる流体の位置と速度です。
これを直接求めるのは難しいので、$p(t - \Delta t)\simeq p(t) - \Delta tv(t)$として現在の速度の反対方向から流れてきたとして、近似的に計算します。
$p(t - \Delta t)$がちょうど格子点上になることはほぼないので、4近傍から補間することで次の速度を求めています。

function advectVelocities(vels, deltaTime) {
	const nextVels = new Array(GRID_NUM + 1);
	for (let xi = 0; xi < GRID_NUM + 1; xi++) {
		nextVels[xi] = new Array(GRID_NUM + 1);
		for (let yi = 0; yi < GRID_NUM + 1; yi++) {
			const pos = createVector(xi * GRID_SIZE, yi * GRID_SIZE);
			const vel = vels[xi][yi].copy();
			vel.mult(-1).mult(deltaTime);
			pos.add(vel);
			pos.div(GRID_SIZE);
			nextVels[xi][yi] = getVelocity(vels, pos.x, pos.y);
		}
	}
	return nextVels;
}

function getVelocity(vels, x, y) {
	if (x < 0.0 || y < 0.0 || x >= GRID_NUM || y >= GRID_NUM) {
		return createVector(0.0, 0.0);
	}
	const ix = floor(x);
	const iy = floor(y);
	
	const fx = x - ix;
	const fy = y - iy;
	
	return p5.Vector.lerp(
		p5.Vector.lerp(vels[ix][iy], vels[ix + 1][iy], fx),
		p5.Vector.lerp(vels[ix][iy + 1], vels[ix + 1][iy + 1], fx),
		fy
	);
}

最後に速度を適当に減衰させることで、何も操作しないと流れが落ち着くようにしています。

function attenuateVelocities(vels, deltaTime) {
	const nextVels = new Array(GRID_NUM + 1);
	for (let xi = 0; xi < GRID_NUM + 1; xi++) {
		nextVels[xi] = new Array(GRID_NUM + 1);
		for (let yi = 0; yi < GRID_NUM + 1; yi++) {
			nextVels[xi][yi] = vels[xi][yi].copy().mult(1.0 - (1.0 - ATTENUATION_RATE) * deltaTime);
		}
	}
	return nextVels;
}

色の更新

色も速度と同様に移流させます。

function updateColors(deltaTime) {
	const nextCols = new Array(GRID_NUM + 1);
	for (let xi = 0; xi < GRID_NUM + 1; xi++) {
		nextCols[xi] = new Array(GRID_NUM + 1);
		for (let yi = 0; yi < GRID_NUM + 1; yi++) {
			const pos = createVector(xi * GRID_SIZE, yi * GRID_SIZE);
			const vel = velocities[xi][yi].copy();
			vel.mult(-1).mult(deltaTime);
			pos.add(vel);
			pos.div(GRID_SIZE);
			nextCols[xi][yi] = getColor(pos.x, pos.y);		
		}
	}
	colors = nextCols;
}

function getColor(x, y) {
	if (x < 0.0 || y < 0.0 || x >= GRID_NUM  || y >= GRID_NUM) {
		return color(255);
	}
	const ix = floor(x);
	const iy = floor(y);
	
	const fx = x - ix;
	const fy = y - iy;
	
	return lerpColor(
		lerpColor(colors[ix][iy], colors[ix + 1][iy], fx),
		lerpColor(colors[ix][iy + 1], colors[ix + 1][iy + 1], fx),
		fy
	);
}

描画

各格子を四角形で描画していきます。色は格子点が持っているので、格子の色は四角形の４つの頂点の値をgetColor(xi + 0.5, yi + 0.5)で補間して求めます。

function render() {
	background(255);
	noStroke();
	for (let xi = 0; xi < GRID_NUM; xi++) {
		const xj = xi + 1;
		for (let yi = 0; yi < GRID_NUM; yi++) {
			const yj = xi + 1;
			if (!RENDER_VELOCITIES) {
				const c = getColor(xi + 0.5, yi + 0.5);
				fill(c);
			} else {
				...
			}
		}
	}
	rect(xi * gridStep, yi * gridStep, gridStep, gridStep);
}

余談

今回、「流体シミュレーション風」とわざわざ言っているのは、真面目に流体シミュレーションをするならば考慮しなければならないいくつかの点を無視しているためです。

流体シミュレーションは以下のナビエ・ストークス方程式と連続の式を満たす必要があります。

\frac{\partial \bf{v}}{\partial t} + (\bf{v} \cdot \nabla) \bf{v} = -\frac{1}{\rho} \nabla p + \nu \nabla^2 \bf{v} + \bf{f}

\nabla \cdot \bf{v} = 0

1つ目の式の左辺第2項が移流項、右辺第1項が圧力項、第2項が粘性項、第3項が外力項となっています。今回実装したものは移流項と外力項しか考慮しておらず、圧力項、粘性項、そして連続の式を考慮していません。すべての式と項をちゃんと計算すると、よりリアルな流体になるはずです。

流体シミュレーションの手法には流体を格子に区切って計算するオイラー的な方法と、流体を粒子とみなして計算するラグランジュ的な方法に分けることができます。

オイラー的な方法でよくCG分野で使われるものにStable Fluidsがあります。今回解説した移流の計算はStable Fluidsを参考にしました。Stable Fluidsの論文にはC言語で実装が紹介されているので、それを参考に愚直に実装するとそれっぽいものが作れます。

私が以前にp5.jsで実装してみたものも置いておきます。
https://www.openprocessing.org/sketch/516944

圧力項や連続の式を真面目に解こうとすると、ポアソン方程式を解いたりと複雑な処理がたくさんあるので、この記事では解説しませんでした（というか、できませんでした...）。

さらに余談

流体を粒子とみなすラグランジュ的な手法の一つにSPH法があります。私が以前にp5.jsで実装してみたものも置いておきます。
https://www.openprocessing.org/sketch/537549