前回は粒子の可視化について説明し、SPH法のプログラムを完成させました。
今回は、この講座の最終回です。回転するブレードを作って水がかき上げられていく様子をシミュレーションしてみます。
目次
ブレードの回転運動
ブレードは水槽の中央付近にあり、ブレードの中心の周りに回転するものとします。水槽の壁と同様に、ブレードも粒子で作成します。ブレード回転の角速度を$\omega$とすると、一時間ステップ$\Delta t$あたりの回転角度は、
$$\Delta \theta = \omega \Delta t \tag{12-1}$$
と表されます。
ブレードを構成する粒子の位置を$(x, y)$、ブレード中心を$(c_x, c_y)$とすると一時間ステップ後の位置$(x', y')$は、
$$\begin{pmatrix} x' \\ y' \end{pmatrix} = \begin{pmatrix} \cos \Delta \theta & -\sin \Delta \theta \\ \sin \Delta \theta & \cos \Delta \theta \end{pmatrix} \begin{pmatrix} x - c_x\\ y - c_y \end{pmatrix} + \begin{pmatrix} c_x\\ c_y \end{pmatrix} \tag{12-2}$$
となります。
また、その時の回転運動の速度は${\bf v} = {\bf \omega} \times {\bf r}$より、
$$\begin{pmatrix} v_x \\ v_y \end{pmatrix} = \begin{pmatrix} -\omega (y' - c_y) \\ \omega (x' - c_x) \end{pmatrix} \tag{12-3}$$
と書くことができます。
プログラム
回転するブレードのプログラムは、ダム崩壊問題のプログラムにブレードに関する部分を付け加えて作っていきましょう。
setParameter関数
まず、ブレードの計算条件をsetParameter関数に加えます。
// blade
regionBlade = new Rectangle(0.1, 0.25, 0.4, thicknessWall);
const rpmBlade = 60;
omegaBlade = rpmBlade / 60 * 2 * Math.PI;
angleBlade = timeDelta * omegaBlade;
const cx = regionBlade.left + 0.5 * regionBlade.width;
const cy = regionBlade.bottom + 0.5 * regionBlade.height;
centerBlade = new Vector(cx, cy);regionBladeはブレードの領域、rpmBladeは回転数[rpm]、omegaBladeは角速度$\omega$[rad/s]、angleBladeは一ステップあたりの回転角度$\Delta \theta$[rad]、centerBladeはブレードの中心座標$(c_x, c_y)$を表します。
ここで、グローバル変数は、regionBlade、omegaBlade、angleBlade、centerBladeです。
initialParticle関数
次に、initialParticle関数で、初期位置でのブレード粒子pBladeを作成します。
// blade particle
pBlade = [];
create(pBlade, regionBlade);pBladeはグローバル変数にしています。
setParticleToCell関数
setParticleToCell関数で、ブレード粒子を空間分割法のバケットに追加します。
// blade particle
cell.add(pBlade);densityPressure関数
densityPressure関数では密度・圧力の計算にブレードを追加します。
calcDP(pBlade);motionUpdate関数
そして、ブレードの回転運動はmotionUpdate関数に記述します。
for (let i = 0, n = pBlade.length; i < n; i++) {
let rv = pBlade[i].position.sub(centerBlade);
pBlade[i].position.x = rv.x * Math.cos(angleBlade) - rv.y * Math.sin(angleBlade) + centerBlade.x;
pBlade[i].position.y = rv.x * Math.sin(angleBlade) + rv.y * Math.cos(angleBlade) + centerBlade.y;
rv = pBlade[i].position.sub(centerBlade);
pBlade[i].velocity.x = -omegaBlade * rv.y;
pBlade[i].velocity.y = omegaBlade * rv.x;
}各ブレード粒子に対して、(12-2)、(12-3)式より、位置positionと速度velocityを計算します。
drawParticle関数
最後に、drawParticle関数に、ブレード粒子の可視化を追加します。
// blade particle
drawEllipse(pBlade, scale, d);回転ブレードプログラムのソースコード
こちらも同様に、JavaScriptのプログラムとp5.js、index.htmlを用意して実行してください。
こちらもソースコードはGitHubに置いています。
→プログラムはこちら（GitHub）
※bladeのフォルダーを参照してください。
まとめ
本講座では、粒子法の一種であるSPH法で流体シミュレーションのプログラムを作ってきました。今回は、JavaScriptと可視化ライブラリのp5.jsを使って、ブラウザ上で動くプログラムを作りました。クラスを使ったオブジェクト指向にもふれたり、盛りだくさんの内容だったと思います。
SPH法は数式も多く出てきてややこしかったかも知れません。これでも、かなり簡易化したモデリングをしています。実際にはもっと物理的にややこしかったりしますが、様々な書籍や論文がでていますので、さらに勉強したい方はいろいろ調べてみてください。