球状星団1個の星を対象として、それを流体とみなして、流体力学で用いられるナビエ・ストークス方程式を適用させます。球状星団は全体的に丸い形をしていて、中心では重力が強く働いていることを考えるとそれに対抗する圧力が必要です。この圧力は熱勾配による圧力なのですが、熱勾配による圧力と重力の力が等しいものであると考えます。よって今回は重力ポテンシャルを今回用いる運動方程式の圧力勾配に用います。また速度に比例する粘性抵抗が星にかかるとします。そして外力は自己重力多体系を簡単に扱うために、一次元調和振動子（いわゆるばねポテンシャル）を用います。
こうして得られたナビエ・ストークス方程式は非線形2階微分方程式になります。線形的に解けないため、数値計算が必要になります。今回は1次元上で星がどのように運動するかをシミュレーションしましたが、その数値計算にルンゲクッタ法を用いました。
以下にそのグラフを載せます。(gnuplotにより作成)
図1 横軸は時間、縦軸は距離を表します。2本の青い線は初期値0.5であり、大きく外側に動いている線はあまり大きく動いていない線に比べて圧力勾配の項目を2倍にしています。2本の緑色の線は初期値0.2であり、圧力勾配を青い線と同じように設定しています。この2本ずつの青い線及び、緑色の線は軽い物体を想定しています。それに対して2本の赤い線は初期値を0.5、0.2にそれぞれ設定し、圧力勾配の項を消しています。圧力勾配の項目は流体密度に反比例して値が小さくなるため、重い星は圧力勾配がかかりづらいとして圧力勾配の項目を消しています。上の図からわかることは各星において減衰振動をしていることです。また圧力勾配がかかっていない流体密度が高くて重い星は、最終的に球状星団の中心に沈み込み、圧力勾配がかかっている流体密度が低くて軽い星は、その圧力勾配の値によって最終的な位置が決まることを今回のシミュレーションでは示しています。
図2 縦軸は速度、横軸は時間を表します。線の色は図1に対応させています。
圧力勾配のかかっている軽い星は速度が大きく変化します。それに対して圧力勾配のかかっていない重い星は速度の変化分が少ないことがわかります。
