下心あってラグランジュ描像のシミュレーションがしたい。流体力学のいちばん最初の授業でやるオイラー描像とラグランジュ描像のあれ。結局そこでしか出てこなくて その後の話は 99% オイラー描像になるというあれ。そういえば海洋ってあんまり混ざっていないという話をしたときにそのような論文を読んだのだった。 調べてみるとこの Haertel さん が職人芸的にこつこつとやっていてどうやら Haertel and Randall (2002) あたりにアルゴリズムが 詳しく述べられている。読んでみると「海水入り変形可能袋」(!)が流れていくという仕掛けだった。これでいいんだろか。 もそっと調べてみると Smooth Particle Hydrodynamics (Monaghan, 2012) の海洋版という位置づけと見た。
いやいや、ちゃんと定式化からどっぷりLagrangian でやってほしいと思って図書館に行って教科書を借りてきた。 流体力学のいちばん最初の授業でやったが、最初時刻 \(s\) の時の座標 \(\boldsymbol{a}\) で流体に名札を付ける。その流体が時刻 \(t\) のときの場所が \[ x_i(\boldsymbol{a}, s|t) \] である。でシミュレーションは運動方程式を解かねばいかんからそれどうなってるの、というと \[ \rho(\boldsymbol{a},s|t)\frac{\partial u_i}{\partial t} (\boldsymbol{a},s|t) = -\frac{\partial A_k}{\partial x_i}(\boldsymbol{x}, t|s)\frac{\partial p}{\partial a_k}(\boldsymbol{a}, s|t) \] となってるそうだ。式 (3.7) である。ここで \[ A_k(\boldsymbol{x}, t|s) \] というのは時刻 \(t\) で座標\(\boldsymbol{x}\) の流体がその昔時刻 \(s\) の時の場所を与える関数。 ここで右辺の微分が \(x_k\) というのが泣けるところ。だってこれ従属変数。ところが、よく考えると \(A_k\) は \(x_i\) の逆関数だから \[ \frac{\partial(x_1,x_2,x_3)}{\partial(a_1,a_2,a_3)}\frac{\partial(A_1,A_2,A_3)}{\partial(x_1,x_2,x_3)} = I\]（ヤコビアンと単位行列） ですなわち \[ \rho\frac{\partial x_i}{\partial a_k}\frac{\partial u_i}{\partial t} = -\frac{\partial p}{\partial a_k} \]（引数は省略） となる。 これを離散化すればよい理屈だが、左辺の\(\partial x_i / \partial a_k\) が大変危険な気配を醸し出す。これもともと隣同士だったグリッドがぐいぐい離れてしまう場合（リヤプノフ不安定） 爆発しそう。そのためには \(\boldsymbol{a}\) で内挿してうまいこと微分を取り直したりせねばいかんのではないか。これが "A purely Lagrangian method for simulating the shallow water equations..." の(8) 式のあたりの事情と読める。 あわてて導入部を読むと "the equations reduce to a smooth particle hydrodynamics" と書いてある。このあたりを追っていくと結局最初の段落に戻るような。車輪見つけた。
Lagrangian は運動方程式がいきなり積分できたりしてなんだかとても素敵なものに見えるのだけれど隣の芝であったか。くやしいので積分をメモしとく。 Cauchy-Weber integral というそうな。 エンタルピー\(dh = dp/\rho\) を定義してやると、 \[ \frac{\partial x_i}{\partial a_j} u_i = -\frac{\partial\phi}{\partial a_j} + v_j \] ここに \[ \phi = \int_s^t \left( h - \frac{1}{2}u_k u_k\right)\mbox{d}r \] で \(v_j\) は時刻\(s\) の流速 \(v_j(\boldsymbol{a},s)\)。(3.46) 式。