ナビエ・ストークス方程式：円柱周りの流れ

レイノルズ数を設定すると、流体がどちらに回転しているかを色で示しながら、円柱周りの流れに対するナビエ・ストークス方程式をリアルタイムで解き進めます。図の下にあるグラフは、3直径分下流における横方向の速度の軌跡です。境界付近では画像から判断しにくい「後流が一定の振幅に落ち着くのか、減衰するのか、あるいは永遠に揺れ続けるのか」を、このグラフが教えてくれます。

計算例 — ナビエ・ストークス方程式：円柱周りの流れ

ソルバが用いる粘性（格子単位）	0.01000
緩和時間 τ	0.5300
後流が落ち着く振幅	0.636 U
ストローハル数（振幅の周期から）	0.1690
抗力係数（時間平均）	1.489
うち圧力によるもの	1.157
抗力に占める表面摩擦の割合	22%
ストローハル数（Roshkoの風洞実験）	0.1671

ページを開いたときに計算機が表示する値です。上の入力を変えると再計算されます。

インタラクティブシミュレーションを読み込んでいます...

レッスン

理論 — ナビエ・ストークス方程式：円柱周りの流れ

ナビエ・ストークス方程式とは、流体に対するニュートンの運動方程式（第2法則）です。微小な流体要素のそれぞれが、押し付けてくる圧力と周囲の流体からの粘性摩擦とによって加速する様子を表しています。圧縮できない流体の場合、方程式は ∂u/∂t + (u·∇)u = −∇p/ρ + ν∇²u および ∇·u = 0 となります。ここで、u は速度、p は圧力、ρ は密度、ν は動粘度を指します。

各記号の意味

Re: レイノルズ数 UD/ν。流速 U、円柱の直径 D、動粘度 ν から計算されます。ステップ1で示すように、方程式に残る唯一の定数です。
St: ストローハル数 fD/U。U/D を単位とした渦放出周波数 f のことです。
C_D: 抗力係数。単位長さの円柱に働く抗力を ½ρU²D で割った値です。
A: 振動の振幅。円柱の3直径分後方における横方向速度の最大値と最小値の差の半分を、U に対する割合で表したものです。
ω: 渦度 ∂v/∂x − ∂u/∂y。流体要素が回転する速さの2倍に等しくなります。図では色として表現されています。

公式の導き方

なぜスライダーが一つで十分なのか。 長さを直径で、速度を U を単位として、時間を D/U を単位として測ってみましょう。つまり、u = U u′、x = D x′、t = (D/U) t′、p = ρU² p′ と置きます。すると、運動量方程式の左辺のすべての項は U²/D という係数を持ち、圧力項も同じ係数を持ちます。一方で粘性項は νU/D² を持ちます。全体を U²/D で割ると、方程式はダッシュ付きの変数で書かれた以前と全く同じ形になりますが、ν が ν/(UD) = 1/Re に置き換わっている点だけが異なります。Re が同じである二つの流れは、大きさと速度が異なるだけの全く同じ流れです。そのため、このページでは4つではなく一つだけパラメータを調整できればよいのです。
このページは、安定と減衰をどう見分けているか。 振動の振幅を A とします。振幅が小さいうちは指数関数的に成長または縮小し（dA/dt = σA）、σ は境界を超えれば正、下回れば負になります。そして境界付近では、σ は Re − Re_c に比例します。しかし、成長は永遠には続きません。上下反転させても同じに見えるような振動において、成長を止める効果を持つ最初の項は3次の項になります。これがランダウの方程式 dA/dt = σA − ℓA³ です。A で割ると、A² に対して直線の関係になる d ln A/dt = σ − ℓA² が得られます。このページでは、各実行の周期にその直線を当てはめて σ を読み取っています。Re 48 において、小さくなっていく振動が消えてしまうのではなく一定の振幅に落ち着こうとしていると判定できるのはこのためです。dA/dt = 0 と置けば、後流が落ち着く振幅 A∞² = σ/ℓ が求まります。これは Re に対して直線となり、Re_c でゼロに到達します。これが2番目のグラフに引かれた線の正体です。
粘性がなければ抗力はゼロ。 粘性が存在しない場合、半径 R の円柱周りの流れは ψ = U(r − R²/r) sin θ となり、表面上の速度は 2U sin θ になります。つまり、前面と背面でゼロ、側面では自由流れの2倍の速さです。流れに沿って p + ½ρv² が一定となるベルヌーイの定理により、これを圧力に変換すると C_p = (p − p∞)/(½ρU²) = 1 − 4 sin²θ となります。抗力とはこの圧力が表面を押す力であり、円周全体にわたって足し合わせることで求まります。計算すると C_D = ∫₀^π C_p cos θ dθ = ∫₀^π cos θ dθ − 4∫₀^π sin²θ cos θ dθ = 0 − [4 sin³θ/3]₀^π = 0 となります。前半部と後半部が完全に相殺されるわけです。このページではソルバ自身の圧力を同じように積分し、抗力のうち圧力によるものとして表示しています。
周波数について。 ステップ1のようにスケール変換すると、渦放出周波数 f は St = fD/U となります。そして、スケール変換された方程式の他の要素と同様に、これは Re のみに依存します。あるレイノルズ数において St が一度わかれば、同じ Re を共有するあらゆる流体、速度、直径に対して f = St·U/D が成り立ちます。Roshkoは Re 50 から 150 の範囲の風洞実験による測定値に St = 0.212(1 − 21.2/Re) を当てはめました。本ソルバの St は Re 55 ではそれより 4% 上回っていますが、Re 100 では 1% の差まで縮まります。

前提: 2次元モデルを仮定しています。無限に長い円柱であり、長手方向の振る舞いには一切の変化がありません。扱うのは圧縮できない流体ですが、ソルバは流体内部での音速の約6分の1で計算を走らせることでこの性質を近似しています。また、幅が直径の20倍ある水路の上面と下面が繋がっているものとしています。これは20直径分だけ離れた場所に別の円柱が並んでいる状況に等しく、その影響で抗力と周波数が数パーセント押し上げられています。
成り立たない場合: Re 188.5 を超えると、実際の後流には円柱に沿って波が発生し（Barkley and Henderson, 1996）、2次元の解では現実を描写できなくなります。また境界付近では、あらゆる変化が遅くなります。Re 47 においては 300 D/U が経過しても振動はまだ減衰の途中にあり、完全に消え去るまでにはさらに 2,000 D/U 以上の時間を要します。だからこそ、このページでは結果を待ち続けるのではなく、ランダウの方程式を用いてその後の傾向を読み取っているのです。

Re 47付近より下では後流は静止し、それを超えると永遠に振動する 🖖

Re 44 を押し、次に Re 48 を押してみてください。円柱は対称であり、そこへ向かう流れも対称です。そしてどちらのレイノルズ数においても、方程式には対称な解、つまり鏡合わせのまま静止した二つの渦が存在します。44 の場合、乱された後流はその対称な解へと戻りますが、48 では二度と戻りません。対称な解そのものは存在し続けているものの不安定になり、ほんのわずかなきっかけで流れがそこから外れてしまうのです。ここでは、設定を変えるたびに流入する流れに横方向の揺さぶりをかけることと、円柱を格子の中心線から1/3セルだけずらして配置していることが、その「きっかけ」となります。水路の実験であれば、どんな微小な振動でも同じ役割を果たします。トゥールーズのカミシェルらのグループはこの境界を 47 と測定しました。本ソルバの増幅率の符号は 47 と 48 の間（約 47.4）で反転し、2番目のグラフの当てはめからは 47.5 が得られます。

粘性がなければ円柱に抗力は働かない 🖖

方程式から粘性を取り除けば、円柱周りの流れは鉛筆と紙だけで解くことができます（下のステップ3）。背面にかかる圧力は前面にかかる圧力の鏡像となり、押し引きが相殺されて抗力は完全にゼロになります。ダランベールがこれを 1752 年に発表して以来、現実の物体には必ず抗力が働くことから「ダランベールのパラドックス」と呼ばれてきました。最初のグラフは、計算上のどこに現実とのズレが生じているかを示しています。Re 100 において、摩擦のない理想的な流れでは圧力が円柱の背面で +1 まで回復するのに対し、ソルバの計算結果は −0.94 で止まっています。これは、圧力が回復する前に流れが表面から剥離してしまうからです。この失われた圧力が抗力の4分の3以上を占めています。残りは流体が表面を引っ張る力であり、抗力全体の 22% (Re 100 の場合) あるいは 40% (Re 20 の場合) に相当します。

風鳴りの音程は一定ではない 🖖

St = fD/U とは、流れが円柱の1直径分を通過する時間を単位とした渦放出周波数のことです。この値はしばしば 0.2 と丸められます。しかしこのページでは、Re 50 で 0.128 だった値が 100 で 0.169、180 で 0.192 へと上昇しています。緑色の線で示されたRoshkoの風洞実験の近似値も同様に上昇し、Re 100 では 0.167 となっています。0.2 という数字は、もっと速い流れに適用されるべきものです。Roshkoが Re 300 から 2,000 の範囲で行った近似では 0.203 から 0.211 となりますが、そこでは後流が3次元になっており、2次元のソルバでは追跡できません。実際の後流は Re 188.5 で3次元になるため、このページのスライダーは 180 で止めてあります。

全プロセスの詳細解説

風の中で電線はどんな音程で鳴るのか？ 4 ステップ
太さ 3 mm の張られた電線が、8 m/s の風の中にあります。 渦は左右から交互に放出され、電線はその渦放出周波数でうなりを上げます。さて、周波数はいくつになるでしょうか。20 °C の空気の粘度は 1.813×10⁻⁵ Pa·s、密度は 1.204 kg/m³ です。

この問題をツールで開く
1. まず、動粘度 ν = η/ρ からレイノルズ数を求めます。
  
  $ν = \frac{η}{ρ} = \frac{1 . 813 \times 1 0 ^{- 5} Pa s}{1 . 204 kg/ m ^{3}} = 1.506 \times 1 0^{- 5} m^{2} /s, Re = \frac{U D}{ν} = \frac{8 m/s \times 0 . 003 m}{1 . 506 \times 1 0 ^{- 5} m ^{2} /s} = 1, 594$
2. これは、実際の後流が3次元に移行する 188.5 を優に超えているため、このページのソルバでは答えを出せません。しかし、Roshkoはまさにこの範囲で渦放出周波数を測定し、次のような近似式を導いています。
  
  $St = 0.212 (1 - \frac{12 . 7}{Re}) = 0.212 (1 - \frac{12 . 7}{1 , 594}) = 0.2103 (300 < Re < 2, 000)$
3. 次に、St の定義から f を計算します。
  
  $f = \frac{St U}{D} = \frac{0 . 2103 \times 8 m/s}{0 . 003 m} = 561 Hz$
4. このソルバが真価を発揮するのは、風洞実験の結果と直接比較できる Re 150 以下の領域です。Re 100 を押し、結果の表から両方の数値を読み取ってみてください。
  
  $Re = 100 : St_{solver} = 0.1690, St_{Roshko} = 0.212 (1 - \frac{21 . 2}{100}) = 0.1671$
解答

電線は約 561 Hz の音程で鳴ります。太さを半分にすると周波数は 1,113 Hz となります。2倍にわずかに届かないのは、Re も半分になり、St が 0.2103 から 0.2086 へと下がるためです。

$f \approx 561 Hz; D = 1.5 mm : Re = 797, St = 0.2086, f = 1, 113 Hz$

参考文献 (5)

Wind-tunnel shedding frequencies behind circular cylinders from Re 40 to 10,000, and the two fits this page plots and the worked problem uses: S = 0.212(1 − 21.2/R) for 50 < R < 150 and S = 0.212(1 − 12.7/R) for 300 < R < 2,000: A. Roshko, "On the Development of Turbulent Wakes from Vortex Streets", NACA Report 1191, National Advisory Committee for Aeronautics, 1954.
Where a two-dimensional solver stops describing a real wake: the periodic 2D wake becomes unstable to three-dimensional disturbances at Re = 188.5 ± 1.0, which is where this page's slider stops: D. Barkley and R. D. Henderson, "Three-dimensional Floquet stability analysis of the wake of a circular cylinder", Journal of Fluid Mechanics 322, 215–241 (1996).
The history of the instability: Bénard's experiments, the threshold Re_c = 47 that Camichel's group in Toulouse measured with circular cylinders, and Mathis, Provansal and Boyer's 1984 measurement that the oscillation's amplitude grows as (Re − Re_c)^½, the Landau law this page fits: J. E. Wesfreid, "Henri Bénard: Thermal convection and vortex shedding", Comptes Rendus Mécanique 345, 446–466 (2017).
The method this page solves the equations with: a lattice-Boltzmann fluid with a curved wall placed on the true circle by interpolating along each lattice link, and the drag from the momentum each link exchanges with the wall: M. Bouzidi, M. Firdaouss and P. Lallemand, "Momentum transfer of a Boltzmann-lattice fluid with boundaries", Physics of Fluids 13, 3452–3459 (2001).
Steady flow past a cylinder computed up to Re 100, with the drag, the surface pressure and the separation bubble, which first appears just below Re 7: S. C. R. Dennis and G.-Z. Chang, "Numerical solutions for steady flow past a circular cylinder at Reynolds numbers up to 100", Journal of Fluid Mechanics 42, 471–489 (1970).

42Math チーム · 公表された資料と照合しています。訂正のご指摘を歓迎します。