エクセルでのシミュレーション

ここでは、Excelを使って数値計算を実施してみましょう。 数値計算といっても比較的簡単なものを演習しますから、難しいと決めつけないで挑戦してみて下さい。

ここで、演習するのは、慣性振動のシミュレーションです。 皆さんご存知の通り地球は自転しています。 フーコーの振り子は、地上にいる私たちにとっては、見かけ上、触れる方向が地球の自転にともなって変化し、地球が 自転していることを教えてくれます。

この、回転している系ではたらく、みかけの力をコリオリの力と呼びます。コリオリの力については、Youtubeの動画でも学ぶ事が出来ます。

コリオリの力は、水平面に立てた垂直の軸の回転が、軸の真上から見て反時計回りである北半球では、運動の方向に対して直角右方向にはたらきます。 南半球では、軸の回転が時計回りになるので、コリオリの力は、運動の方向に直角左方向にはたらきます。 このコリオリの力によって、風が回転して吹く場合、表層の流れも、しばしば、地球の自転と同じ周期で（水平面に立てた垂直の軸の回転周期）くるくると回る様な流れになります。下図(図の参照 )は海面に浮かべたブイの軌跡を時間を追って示したものです。

Figure: Van Meurs 1998からの引用　表層ドリフターの軌跡

　　

　　

この地球の自転に伴う、慣性振動を数値的に解く事を演習してみましょう。 高校から学んできた、運動の方程式は、$Ma=F$だったと思います。 ここで、$M$は質量、$a$は加速度、$F$は外力です。外力$F$によって 質量$M$の物体が、加速度$a$で運動している事を示しています。

この物体を、流体粒子、海水の一部分としましょう。 ここで扱う海水の密度は一定であるため、両辺を密度で割り、 質量$M$は、現れません。この流体を追跡した運動方程式は、

$$\frac{dU}{dt}=fV$$ (1)

$$\frac{dV}{dt}=-fU$$

となります。ここで、$U$は東西流、$V$は南北流、$f$はコリオリパラメータで、 $f=2\Omega \sin(\phi)$で定義され、 $\Omega$は地球の自転角速度、$\phi$は緯度(ラジアン)になります。ここでは、簡単のために、$f=10^{-4}$とします。 この式の導出は、２、３年生に受講できる授業で習うはずですので割愛します。 右辺に現れた項が、コリオリの力を表します。この式を頭の中で計算してみますと、$V$が正で北への流れがあれば、 $U$（東西流）は、右辺が東西流の時間に対する傾きを示す事から、時間とともに増加します。 また、東西流$U$が負なら南北流$V$は時間とともに増加する事が分かります。つまり、流速は、時計回りに回転 しそうだということがわかると思います。

さて、ここでは、北方向に時間的に変化しない平均地衡流$\bar{v}$が存在すると仮定しましょう。即ち、 もとの$x$方向の運動方程式には、圧力傾度力が加わり

$$\frac{dU}{dt}=fV-\frac{1}{\rho_0}\frac{\partial p}{\partial x}$$ (2)

$$\frac{dV}{dt}=-fU$$

となります。流速場は、以下の様に分けられる事が出来ます。即ち、

$$U=u$$ (3)

$$V=\overline{v}+v$$

で、東西流は、慣性振動する$u$のみで、南北流は、定常的な流れ$\bar{v}$と慣性振動する成分 $v$の和で表されるとしましょう。すると式()は、以下の様に表されます。

$$\frac{du}{dt}=f(\bar{v}+v)-\frac{1}{\rho_0}\frac{\partial p}{\partial x}$$ (4)

$$\frac{dv}{dt}=-fu$$

ここで、以下の地衡流のバランスを用いれば、

$$f\bar{v}=\frac{1}{\rho_o}\frac{\partial p}{\partial x}$$ (5)

運動方程式は、

$$\frac{du}{dt}=fv$$ (6)

$$\frac{dv}{dt}=-fu$$

となり元の式の形に戻ります。 いま、$v$を上式から削除すると、式は、

$$\frac{d^2u}{dt^2}=-f^2u$$ (7)

となり、これはバネの振動等を表す、単振動の式です。式は、線形で同次、定数係数をもつので、解を $e^{\lambda t}$とおいて式()に導入してみましょう。時間のある人は、解を求めてみましょう。 ここで$\lambda$は定数です。

この運動方程式を離散化してみます。離散化とは、コンピューターで行う、四則演算で計算できる様、 方程式を書き直す事です。具体的には、次の様なものが上式を離散化したものになります。

$$\frac{u^{n+1}-u^{n}}{\Delta t}=\frac{f}{2}(v^{n}+v^{n+1})$$ (8)

$$\frac{v^{n+1}-v^{n}}{\Delta t}=-\frac{f}{2}(u^{n}+u^{n+1})$$

ここで、 $n=1,2,3,\cdots$は時間ステップを表し、$n+1$は1ステップ未来の時刻と言う事になります。 即ち、時刻は $t=(n-1) \times \Delta t$と書く事が出来ます。 ここでは、右辺の外力：コリオリの力を求めるときに、未来のステップと現在のステップでの 流速から算出される平均的なコリオリの力を用いるとしておきましょう。しかしながら、 右辺第二項は、未知の未来の流速を使って計算しなければならない事になります。 現実的には、これに近い方法で、次のステップの流速を求める事を考えた方が良さそうです。 例えば、1ステップ未来の流速の推測値を$u^*$として以下の離散方法でこれを求めることにします。

$$\frac{u^*-u^{n}}{\Delta t}=fv^{n}$$ (9)

$$\frac{v^*-v^{n}}{\Delta t}=-f u^{n}$$

即ち、1ステップ未来の推定値は、簡単に現在の流速に伴うコリオリの力から求める事が出来ます。

$$u^*=u^n+\Delta t fv^{n}$$ (10)

$$v^*=v^n-\Delta t fu^{n}$$

次に、これを式()の$u^{n+1}$や$v^{n+1}$に代入し、右辺を1ステップ未来の推定値と現在の 流速から求める事にしましょう。

$$u^{n+1}=u^{n}{+{\Delta t} \frac{f}{2}(v^{n}+v^*)$$ (11)

$$v^{n+1}=v^{n}-{\Delta t} \frac{f}{2}(u^{n}+u^*)$$

このような若干面倒な方法で次のステップを計算するのは、この方法が比較的 誤差が小さく(Second Order)、計算も安定し易いためです。

さて、この問題を、$u$の初期値を0.2 $m/s$とし、$v$の初期値を0として、計算しましょう。 時間ステップは、100秒で行います。前述の通り、コリオリパラメータ、$f=10^{-4}$とします。

まず、エクセルの新しいファイルを開いて数の様に各列に$u$と$u^*$、$v$と$v^*$、等とともに 計算に必要な定数を記載します。

Figure: 数値計算

　　

　　

ではまず、$u$の1ステップ未来の推定値$u*$を求めましょう。 そのためには、式()にならって、

=A2+$F$2*$E$2*C2

とします(図の参照)。マークは絶対参照を意味します。

Figure: 数値計算

　　

　　

$v^*$も同様に行いますが、符号が先ほどと異なる事に注意しましょう。

つづいて、$u$の次のステップの計算を式()に習って、先ほど計算した推定値$v^*$を用いて 以下の様に実施しましょう。

=A2+$F$2*$E$2*0.5*(C2+D2)

これと同様に$v$の次のステップの計算も行いましょう。先ほどと同様に符号が異なる点も注意して下さい。 これらが出来たら、新しく計算したステップの$u$と$v$を使って、更に1ステップ未来の推定値を計算しましょう。 絶対参照しているので、下にセルをコピーすれば可能です。計算がうまくいっていれば、3行目以降はずっと先の 未来まで計算してみましょう。864ステップ計算すれば、86400秒の計算ですので1日になります。

計算された$u$と$v$の列を選択してグラフを作成しましょう。散布図でよいです。

さて、これは流体粒子の流速を示しています。 次に、粒子の軌跡を計算しましょう。最初の方で述べた様に、背景には、平均的な 地衡流に伴う流れ$\bar{v}$が北方向にあるとしました。それを加えた南北流$V=\bar{v}+v$と 東西流$U=u$を、今度は時間的に積分します。$\bar{v}=0.05$ $m/s$としましょう。即ち、

$$x=x_o+\int U dt$$ (12)

$$y=y_o+\int V dt$$

$x_o$と$y_o$は両方とも0とします。 先ほどと同様、流速を2ステップ使ってその平均値を 時間で積分し、その時間に進んだ、$x$と$y$方向の距離を計算します。即ち

$$x=x^n+\frac{\Delta t}{2} (u^{n+1} + u^n)$$ (13)

$$y=y^n+\frac{\Delta t}{2} (v^{n+1} + v^n)$$