初めに

これまでに開発してきた有限差分法コード（cf. こちらやこちら）を障害物周りの流れに拡張し，von Kármánの渦列をシミュレートする．

Methods

支配方程式

非圧縮性流体（Newtonian）を考える．無次元化しないまま話を進める．

\begin{alignat}{2}
    \nabla \cdot \boldsymbol{u}
    &= 0
    &&\quad \text{in} \quad \Omega \times \mathcal{I} \\
    \frac{\partial \boldsymbol{u}}{\partial t}
    + \left( \boldsymbol{u} \cdot \nabla \right) \boldsymbol{u}
    &= - \frac{1}{\rho} \nabla p
    + \nu \nabla^2 \boldsymbol{u}
    + \boldsymbol{f}
    &&\quad \text{in} \quad \Omega \times \mathcal{I} \\
\end{alignat}

ここで，$\boldsymbol{u}$は速度$[ \mathrm{m / s} ]$，$p$は圧力$[ \mathrm{Pa} ]$，$\rho$は密度$[ \mathrm{kg / m^3} ]$，$\nu$は動粘度$[ \mathrm{m^2 / s} ]$である（$\rho$と$\nu$は定数とする）．$\boldsymbol{f}$は外力$[ \mathrm{m / s^2} ]$であり，今回はゼロとする．

パラメータには，

密度　: $\rho = 1 \times 10^{0} \ [ \mathrm{kg / m^3} ]$
動粘度: $\nu = 1 \times 10^{-2} \ [ \mathrm{m^2 / s} ]$

を選ぶ（動粘度は場合によって変える）．このように設定しておくと，粘度$\mu \ [ \mathrm{Pa \cdot s} ]$と動粘度$\nu \ [ \mathrm{m^2 / s} ]$が同じ数値となり便利である（$\mu = \rho \nu$）．また，特徴長さと特徴流速に$L = 1 \ [ \mathrm{m} ], \ U = 1 \ [ \mathrm{m / s} ]$を取れば，$\mathrm{Re} = \nu^{-1} = 100$となるように，Reynolds数の設定が簡便になる（特徴長さと特徴流速を固定して動粘度だけ変えれば良い）．一見，「一体全体どのような流体を考えているのだ？」と感じるが，既に見たように，Reynoldsの相似則の範囲内では現象は本質的に同義であり，したがって，例えば，水の物性値（cf. こちら）を用いるときには，特徴長さや特徴流速を大きく設定して同じだけのReynolds数を達成すれば良い．

離散化

Chorinの射影法（Chorin1968）により速度と圧力をカップリングする．久しぶりに次元のある式を見つめるので，パラメータを忘れないよう注意深く進む．運動量の式を以下のように分離する．

\begin{align}
    \frac{\boldsymbol{u}^{(*)} - \boldsymbol{u}^{(n)}}{\tau}
    &= 
    - \left( \boldsymbol{u}^{(n)} \cdot \nabla \right) \boldsymbol{u}^{(n)}
    + \nu \nabla^2 \boldsymbol{u}^{(n)}
    + \boldsymbol{f}^{(n)} \\
    \frac{\boldsymbol{u}^{(n+1)} - \boldsymbol{u}^{(*)}}{\tau}
    &= - \frac{1}{\rho} \nabla p^{(n+1)} \\
\end{align}

ここで，$\tau$は時間刻みである．2式目で発散を取り，$\nabla \cdot \boldsymbol{u}^{(n+1)} = 0$を要求する．

\begin{align}
    \nabla \cdot \left(
        \frac{\boldsymbol{u}^{(n+1)} - \boldsymbol{u}^{(*)}}{\tau}
    \right)
    &= \nabla \cdot \left(
        - \frac{1}{\rho} \nabla p^{(n+1)}
    \right) \\
    \therefore
    \nabla^2 p^{(n+1)}
    &= \frac{\rho}{\tau} \nabla \cdot \boldsymbol{u}^{(*)} \\
\end{align}

即ち，

#	Description	Equation
1	速度の予測	$\boldsymbol{u}^{(*)} = \boldsymbol{u}^{(n)} + \tau \left( - \left( \boldsymbol{u}^{(n)} \cdot \nabla \right) \boldsymbol{u}^{(n)} + \nu \nabla^2 \boldsymbol{u}^{(n)} + \boldsymbol{f}^{(n)} \right)$
2	圧力の求解	$\nabla^2 p^{(n+1)} = \rho \cdot \tau^{-1} \cdot \nabla \cdot \boldsymbol{u}^{(*)}$
3	速度の修正	$\boldsymbol{u}^{(n+1)} = \boldsymbol{u}^{(*)} + \tau \left( - \rho^{-1} \cdot \nabla p^{(n+1)} \right)$

により時間進行する．

空間の作用素は，過去の投稿（こちらやこちら）に倣う．明言しない限り，渦粘性は導入しない．

Results

Hagen–Poiseuille流れ

先ずは簡単に，Hagen–Poiseuilleの流れをテストケースに選ぶ（とは言え，久しぶりすぎて破滅的に忘却していた，，）．

流入速度を規定する場合

$\Omega = (0, L_x) \times (0, L_y) = (0, 4 L) \times (0, L), \ \mathcal{I} = (0, T)$を考え，$L = 1 \ [ \mathrm{m} ], \ T = 10 \ [ \mathrm{s} ]$とする．$h = 5 \times 10^{-2} \ [ \mathrm{m} ], \ \tau = 1 \times 10^{-2} \ [ \mathrm{s} ]$で離散化する．

既に述べたように，$\nu = 1 \times 10^{-2} \ [ \mathrm{m^2 / s} ]$を選んでおり，いまの特徴長さは$L_y = 1 \ [ \mathrm{m} ]$である．$\boldsymbol{U}_{\infty} = [ 1, 0 ]^{\top} \ [ \mathrm{m / s} ]$なる一様流速で流入させる．$\mathrm{Re} \simeq 100$となる．

境界条件は，河野2022: 詳解流れの数値計算のp199辺りの数値計算を参考に，

左端: $\boldsymbol{u} = \boldsymbol{U}_{\infty}, \partial_n p = \boldsymbol{0}$
右端: $\partial_{n} \boldsymbol{u} = \boldsymbol{0}, \ p = 0$
上端: $\boldsymbol{u} = \boldsymbol{0}, \ \partial_n p = \boldsymbol{0}$
下端: $\boldsymbol{u} = \boldsymbol{0}, \ \partial_n p = \boldsymbol{0}$

とした．

初期の速度場を全くゼロとした場合と，乱数（$\sim \mathcal{U} (-0.5, 0.5)$）を用いて摂動を与えた場合の結果を示す．上から，圧力，流速ノルム，流速ベクトル，渦度である．

摂動無し	摂動有り

良く知られているように，十分発達したとき，速度プロファイルはparabolicになる．$(0, 4 L)$の区間では十分に発達しているとは若干言い難い気もするが（もう少しのrun-in区間を設ければ良かった），大らかにこの程度とする．なお，$\mathrm{Re} = 100$より，粘性の影響が中々に強いので，初期の摂動は素早く消散している．

圧力勾配で駆動する場合

ところで，この現象は本来，圧力勾配により駆動される流れである．これを考えると，計算上でも，流入速度を与えるより，圧力勾配を与える方が"自然"であろう．

fully-developedな流れを考え，parabolicな速度分布（$\boldsymbol{U} = \left[ U_{\mathrm{max}} y \left( L_y - y \right) / L_y^2, 0 \right]^{\top} \ [ \mathrm{m / s} ]$）をイメージする．このとき，平均流速$U_{\mathrm{mean}}$は，

\begin{align}
    U_{\mathrm{mean}}
    = \int_{0}^{L_y}
        U_{\mathrm{max}} \frac{y \left( L_y - y \right)}{L_y^2}
    = \frac{1}{6} U_{\mathrm{max}} L_y
    = \frac{1}{6} U_{\mathrm{max}} L
    = \frac{1}{6} U_{\mathrm{max}}
\end{align}

である．特徴流速として$U_{\mathrm{max}}$を採用する場合と$U_{\mathrm{mean}}$を採用する場合が考えられる．

特徴流速に$U_{\mathrm{max}}$を採用する場合
対応するReynolds数を

\begin{align}
    \mathrm{Re}_{\mathrm{max}}
    = \frac{U_{\mathrm{max}} L}{\nu}
\end{align}

と定めておく．$\mathrm{Re}_{\mathrm{max}} = 100$を狙うと，

\begin{align}
    100
    = \mathrm{Re}_{\mathrm{max}}
    = \frac{U_{\mathrm{max}} \cdot 1}{1 \times 10^{-2}}
    \Leftrightarrow
    U_{\mathrm{max}} = 1 \ \left[ \mathrm{m / s} \right]
\end{align}

より，最大流速が$1 \ [ \mathrm{m / s} ]$になるような圧力勾配を与えれば良い．解析解から，

\begin{align}
    U_{\mathrm{max}}
    &= - \frac{1}{8 \rho \nu} \frac{\partial p}{\partial x} L_y^2 \\
    \therefore
    \frac{\partial p}{\partial x}
    &= - 8 \rho \nu \frac{U_{\mathrm{max}}}{L_y^2}
     = - 0.08 \ \left[ \mathrm{Pa / m} \right] \\
\end{align}

なる圧力勾配を与えれば，$\mathrm{Re}_{\mathrm{max}} \simeq 100$．

特徴流速に$U_{\mathrm{mean}}$を採用する場合
対応するReynolds数を

\begin{align}
    \mathrm{Re}_{\mathrm{mean}}
    = \frac{U_{\mathrm{mean}} L}{\nu}
\end{align}

と定める．先程と同様に，

\begin{align}
    100
    = \mathrm{Re}_{\mathrm{mean}}
    = \frac{U_{\mathrm{mean}} \cdot 1}{1 \times 10^{-2}}
    \Leftrightarrow
    U_{\mathrm{max}} = 6 \ \left[ \mathrm{m / s} \right]
\end{align}

であるから，

\begin{align}
    U_{\mathrm{max}}
    &= - \frac{1}{8 \rho \nu} \frac{\partial p}{\partial x} L_y^2 \\
    \therefore
    \frac{\partial p}{\partial x}
    &= - 8 \rho \nu \frac{U_{\mathrm{max}}}{L_y^2}
     = - 0.48 \ \left[ \mathrm{Pa / m} \right] \\
\end{align}

なる圧力勾配を与えれば，$\mathrm{Re}_{\mathrm{mean}} \simeq 100$．

先程と同じ領域を考え，$T = 60 \ [ \mathrm{s} ]$までシミュレートする．$h = 5 \times 10^{-2} \ [ \mathrm{m} ], \ \tau = 2.5 \times 10^{-3} \ [ \mathrm{s} ]$とした．境界条件は，

左右端: 周期境界条件
上下端: 非滑り条件

とした．

$\mathrm{Re}_{\mathrm{max}} \simeq 100$	$\mathrm{Re}_{\mathrm{mean}} \simeq 100$

時刻$T = 60 \ [ \mathrm{s} ]$付近での，1ステップ当たりの速度の変化量や管路内の最大流速は，

$\mathrm{Re}_{\mathrm{max}} \simeq 100$
- 速度の変化量　　: $6.688 \times 10^{-7}$
- 管路内の最大流速: $0.997 \ [ \mathrm{m / s} ]$
- Courant数$C$，拡散数$D$: $C = 0.05, \ D = 0.02$
$\mathrm{Re}_{\mathrm{mean}} \simeq 100$
- 速度の変化量　　: $6.688 \times 10^{-7}$
- 管路内の最大流速: $5.983 \ [ \mathrm{m / s} ]$
- Courant数$C$，拡散数$D$: $C = 0.30, \ D = 0.02$

であり，十分に収束しているし，狙った速度も出ていると言って良いだろう．

Couette–Poiseuille流れ

Hagen–Poiseuille流れにおいて，上端の壁をある流速$U^{\prime} \ [ \mathrm{m / s} ]$で動かすと，pressure-driven（Hagen–Poiseuille）かつshear-driven（Couette）な流れ場になる．

先程のHagen–Poiseuilleの流れについて，parabolicな速度プロファイルの内から，その最大流速$U_{\mathrm{max}}$を特徴流速に採用する場合を思い出す．$\mathrm{Re}_{\mathrm{max}} = 100$を狙い，$\partial p / \partial x = - 0.08 \ \left[ \mathrm{Pa / m} \right]$だけの圧力勾配を与えるのであった．

同じ問題に対し，$U^{\prime} = 1 , -0.5 \ [ \mathrm{m / s} ]$とした結果を示す．

$U^{\prime} = 1 \ [ \mathrm{m / s} ]$	$U^{\prime} = -0.5 \ [ \mathrm{m / s} ]$

Hagen–Poiseuille流れとCouette流れが重なったような状態となった．2つ目の図は，$x = 0.5$での速度のスライスである．解析解（定常解，黒実線）と数値解（赤破線）が良く合致している．

Backstep流れ

打ち出された流れが，どの辺りで着地するかを見る問題らしい（backstep / backward-facing step flow）．この問題では，特徴流速に流入流速を，特徴長さに階段の高さを選ぶようである．

一様流速$\boldsymbol{U}_{\infty} = [ 1, 0 ]^{\top} \ [ \mathrm{m / s} ]$で流入させ，階段の高さ$H$には$H = 0.5 L = 0.5 \ [ \mathrm{m} ]$を取る（expansion ratio: 2）．動粘度を動かしてReynolds数を変化させる．$h = 2 \times 10^{-2} \ [ \mathrm{m} ], \ \tau = 2 \times 10^{-3} \ [ \mathrm{s} ]$で離散化し，$T = 20 \ [ \mathrm{s} ]$までシミュレートする（これは定常問題か？ううむ，良く知らない，，）．

$\mathrm{Re} \simeq 50$	$\mathrm{Re} \simeq 500$

少しだけ調べてみたが，これといった参照解が中々見つけられなかった．Biswas+2004と凡そ似ている，という，定性的かつ大雑把な結論に留める（もう少し横方向に長い方が良かったかもしれない）．しかし，激しい剥離が起きる箇所で負圧が生まれて吸い込みを引き起こすなど，これだけシンプルな形状なのに，複雑に渦を巻くというのは興味深い．

(2024/11追記)
暫く経って，最近は，負圧が吸い込んでいるというより，単に圧力勾配に乗っているという方が正しいかもしれない，と思ってきた．また，境界条件の与え方に誤りがあったので，修正した．コード修正のついでに，可視化に少し工夫を加え，圧力は無次元化圧力$p^{\prime} \left( = p / \left( 1/2 \rho U_{\infty} \right) \right)$を視ることにし，圧力勾配は運動方程式に登場する形式$- \nabla p / \rho$で表示することにした．また，速度はベクトルにノルムで色付けし，渦度はゼロ領域が見えやすいようにfilled contourにした．残念ながら，ファイルサイズ上限から，低解像度化する必要があったので，鮮やかな色使いではない．

$\mathrm{Re} \simeq 50$	$\mathrm{Re} \simeq 500$

また，これに合わせて少し調べたが，この問題はやはり定常解を見るようだ．例えば，Dick1982, Kim&Moin1985, Scott+1986, Rogers&Kwak1991，或いは，こちらがPIVの結果を公開してくれている．中山司2008: 流れ解析のための有限要素法入門の第8章も詳しい．気が向いたときに手元の計算と比較することにしておく．

von Kármán渦列

やっと，本題である（ここに至るまでの簡単な問題たちは，簡単だったものの，境界条件の取り扱いを考えるなどの点で大いに役立った）．気になっていた幾つかの事項に触れていく．

Initial perturbation

改めて，$\Omega = (0, L_x) \times (0, L_y) = (0, 4 L) \times (0, L), \ \mathcal{I} = (0, T)$を考える．流れ場の中央に1辺$D \left( = L / 5 \right) \ [ \mathrm{m} ]$の正方形のシリンダーを配置する（blockage ratio: 0.2, see, e.g., Kumar&Mittal2005, Prasanth+2006, Nguyen&Lei2021, Nvv&Chatterjee2022, Abdul-Salam&Tachie2024 for more about the effect of blockage ratio）．$L = 1 \ [ \mathrm{m} ], \ T = 40 \ [ \mathrm{s} ]$とする．最大流速が$1 \ [ \mathrm{m / s} ]$であるparabolicな流速で流入させる．その平均流速は$2/3 \ [ \mathrm{m / s} ]$である．以降，特徴流速に最大流速を，特徴長さにシリンダー径を選び，動粘度を動かしてReynolds数を変化させる¹．離散化パラメータは，$h = 5 \times 10^{-2} \ [ \mathrm{m} ], \ \tau = 5 \times 10^{-3} \ [ \mathrm{s} ]$とする（かなり大雑把だが，コードの確認のため）．境界条件の取り扱いは，こちらやこちらに準ずる（全く同じという訳ではない）．

特徴長さが$D = 0.2 \ [ \mathrm{m} ]$になっていることに注意して$\nu \ [ \mathrm{m^2 / s} ]$を選び，$\mathrm{Re} \simeq 40, \ 200, \ 2,000$の3つを設定する．渦を誘発するため，初期の速度場を乱数（$\sim \mathcal{U} (-0.5, 0.5)$）で乱しておいた（長時間待てば，きちんと渦が生まれるのだが，私はせっかちだ）．

	$\mathrm{Re} \simeq 40$	$\mathrm{Re} \simeq 200$	$\mathrm{Re} \simeq 2,000$
摂動無し
摂動有り

色の様子がこれまでと異なるのは，アップロード可能なファイル上限に抑えるよう圧縮したためである（が，レトロな雰囲気があって中々愛着が湧くではないか）．

渦が生まれるまでの間は，wake領域で流体から障害物向きの流れが発生し，流れ線が閉じている（教科書のPotential流れの章で学んだことが，計算上本当に起きている！）．また，摂動が無ければ，$t = 40 \ \left[ \mathrm{s} \right]$程度まで待たなければ渦が見えてこないが（40秒で支度しな！），摂動を与えておけば，渦が素早く成長する．安定に見えるモードが突然崩壊し，ある振動のパターンを呈する様子を見ると，何だかbifurcation diagramを思い出す（関係有るかな，関係無いかな，，）．

$\mathrm{Re} \simeq 40$では，粘性の影響が強く，緩やかな流れに落ち着いている．transientな状態に移るまでのReynolds数は$40 \sim 50$程度と知られているので，凡そ正しそうだと大雑把に片付ける．$\mathrm{Re} \simeq 200, \ 2,000$では，wake領域に渦が生成されているが，壁が渦の成長を強く制限しているので，壁の影響をもっと遠方に飛ばしてやろう．

Blockage ratio

昔受講した大学院の講義では，こういう問題では，blockage ratioを3%以下に抑えるのが望ましい，と教わった．$\mathrm{Re} \simeq 200, \ 2,000$の場合について，$\beta$を変化させてみる．空間解像度の議論は後回しにして，$h = 5 \times 10^{-3} \ [ \mathrm{m} ], \ \tau = 5 \times 10^{-4} \ [ \mathrm{s} ]$で離散化する．渦粘性として，Smagorinskyモデル（Smagorinsky1963）を導入し，SGS渦の特徴長さとグリッド幅の比$C_s$は$C_s = 0.1$とした（Lilly1966, Deardorff1970）．次項の議論で触れるため，特徴的な長さ$D$と空間解像度$h$との比率を併記しておく．

$\beta$	$D/h$	$\mathrm{Re} \simeq 200$	$\mathrm{Re} \simeq 2,000$
$0.2$	$40$
$0.1$	$20$
$0.05$	$10$
$0.03$	$6$

3回再生したら停止するようにした．ダンプアウト間隔をもう少し細かくしておくと良かったかもしれないが（$\beta \le 0.05$では，ダンプアウト間隔と"尾を振る"周期が殆ど同じだ），コストが高く，再計算する気力は起きない．

$\mathrm{Re} \simeq 200$について．$\beta = 0.2$では，既に見たように，渦が生まれて直ぐに壁の影響を受けており，その成長が阻害されている．$\beta = 0.1$では，その傾向は緩和されているが，障害物から離れるにつれ，流れが壁の影響を受けていないとは言い難くなっている（どれだけの長さの流路を考えるかにも依るが）．$\beta = 0.05, \ 0.03$では，壁の影響を殆ど受けておらず，渦がかなり自由に発達，移流している．また，$\beta = 0.05$の時点で既に壁の拘束は小さく，$\beta = 0.03$との相違はあまり無いように感じる．

$\mathrm{Re} \simeq 2,000$について．何れのblockage ratioに於いても，渦の成長は壁に制限されている．$\beta = 0.05, \ 0.03$ほどにしても，激しく放出された渦が，直ぐに壁に衝突，減衰，反射している．また，メインの話題からは逸れるが，同符号の渦同士が衝突して合体する様子も見ることができる．

Reynolds数の上昇に伴い，これより小さくすべきというblockage ratio（critical blockage ratioとでも言うべきか？）は変わってくるように感じる．教わった3%以下という指針も，恐らくユニバーサルではないのではないか，と感じる．しかし，やはり計算資源は有限なので，こういう指針の存在はコミュニティにとって非常に重要・有用であると思う．

Spatial resolution

障害物周りの流れをシミュレートするに当たり，必要な解像度はどの程度か気になる．最も特徴的な流れは障害物付近の流れであるから，障害物に対する相対的な解像度で以て議論すべきであろう．計算法が異なるが，例えば，Breuer+2000（LBM, FVM）やAsai+2023（SPH）は，障害物のすぐ傍の流れをどれだけ解像するか，という議論により数値解の収束を主張している．これに倣い，$\mathrm{Re} \simeq 200, \ 2,000$について，blockage ratioを$\beta = 0.05$に固定したまま（$D = 0.05 \ [ \mathrm{m} ]$），相対的な解像度$D/h$を少しずつ変化させたときの結果を示す．先程と同様に，時間分解能は$\tau = 5 \times 10^{-4} \ [ \mathrm{s} ]$とし，Smagorinskyの渦粘性（$C_s = 0.1$）を用いる．

$h \ [ \mathrm{m} ]$	$D/h$	$\mathrm{Re} \simeq 200$	$\mathrm{Re} \simeq 2,000$
$0.01$	$5$
$0.005$	$10$
$0.0025$	$20$

$\mathrm{Re} \simeq 200$では，$t = 5 \ [ \mathrm{s} ]$までの計算を実行した．$\mathrm{Re} \simeq 2,000, \ D/h = 20$では，Courant数が徐々に大きくなり，$t = 2.8 \sim 2.9 \ [ \mathrm{s} ]$周辺で爆発してしまったので，その少し前の$t = 2.0 \ [ \mathrm{s} ]$を見比べる（時間増分を小さくして再計算するには，あまりに計算コストが高い²）．

$\mathrm{Re} \simeq 200$について．$D/h = 5$と$D/h = 10$の間には，顕著な差がある．$D/h = 10$と$D/h = 20$の間には，違いが無いとは言い難い．やはり，$D/h = 10$の方が若干ぼやけていたり，diffusiveな印象がある．しかし，$D/h = 10$は凡そ$D/h = 20$に近いと言って良いだろう（少なくとも顕著な差がある訳ではない）．

$\mathrm{Re} \simeq 2,000$について．$D/h = 5$と$D/h = 10$の間には，やはり顕著な差があり，かなりなまっている．$D/h = 10$と$D/h = 20$を見比べても，放出される渦の概形は似ているが，両者には明確な差があると言えるだろう．ここまで局所的な渦構造を見る必要があるか？という疑問はあるが，Reynolds数に応じて，必要な解像度に差がありそうだということが分かる．

周辺の議論はcircular cylinderで展開されることが多いが，本稿同様にsquare cylinder周りの高Reynolds数での流れを対象に，これをDNSした例として，Trias+2015など．

Boundary conditions

これまで，壁には非すべり条件（no-slip condition）を与えて計算してきた．これは，粘着の条件とも呼ばれるように，静止した壁と粘性流体との間の境界条件としては，最も自然だと考える．一方で，壁が流体に対して相対的な速度を持っている場合や空気と海面の間の境界条件などを考慮する場合には，すべり条件（slip condition）を課す場合もある．また，十分に大きい計算領域を用意する場合や壁から十分に遠方の流れを考える場合などには，周期境界条件（priodic condition）を用いる事もある³．

上下端の境界に対して，

非すべり条件: $\boldsymbol{u} = \boldsymbol{0}, \ \partial_{\boldsymbol{n}} p = 0$
すべり条件　: $\partial_{\boldsymbol{n}} \boldsymbol{u} = \boldsymbol{0}, \ \partial_{\boldsymbol{n}} p = 0$
周期境界条件: $\boldsymbol{u}_u = \boldsymbol{u}_l, \ p_u = p_l$

を与え，流れの様子がどのように変わるかを見る（$\phi_u, \ \phi_l$は，上端，下端での物理量を示す）．$\mathrm{Re} \simeq 2,000, \ \beta = 0.05$の場合について，$h = 5 \times 10^{-3} \ [ \mathrm{m} ], \ \tau = 5 \times 10^{-4} \ [ \mathrm{s} ]$で離散化し（$D/h = 10$），Smagorinskyの渦粘性（$C_s = 0.1$）を用いる．流入条件は一様流速とした．先程まではparabolicな流入流速を与え，その最大流速を特徴流速としてReynolds数を概算していたので，Reynolds数の定義が少し違う．

先ず，$\mathrm{Re} \simeq 200$について．

$\mathrm{BC}$	$t = 1 \ [ \mathrm{s} ]$	$t = 4 \ [ \mathrm{s} ]$	$t = 7 \ [ \mathrm{s} ]$	$t = 10 \ [ \mathrm{s} ]$
$\mathrm{no \ slip}$
$\mathrm{slip}$
$\mathrm{periodic}$

続いて，$\mathrm{Re} \simeq 2,000$の場合について．

$\mathrm{BC}$	$t = 1 \ [ \mathrm{s} ]$	$t = 4 \ [ \mathrm{s} ]$	$t = 7 \ [ \mathrm{s} ]$	$t = 10 \ [ \mathrm{s} ]$
$\mathrm{no \ slip}$
$\mathrm{slip}$
$\mathrm{periodic}$

非すべり条件では流入境界の直ぐ傍の壁面からshearが起きているのに対し，すべり条件ではこれが起きない．但し，壁自体は存在するので，壁に衝突する渦もある．壁に衝突した渦は，非すべり条件では反射するが，すべり条件では壁の表面ですべり続ける．周期境界条件を課すと，壁の影響が陽に現れないので，shearが起きず，さらに，上端から飛び出した渦が下端から登場し（vice versa），wake regionにある他の渦に干渉する．周期境界を課してwake regionが十分に自由な挙動を示すようになるまでblockage ratioを小さくするなどの方針で，適切なblockage ratioを決めても良いかもしれない．$\mathrm{Re} \simeq 2,000$では，$t = 1 \ [ \mathrm{s} ]$では全ての条件で渦の放出パターンで同様であるのに対し，$t = 4 \ [ \mathrm{s} ]$の時点で既にズレが生じているのも面白い．

Advection scheme

特別断り無く，移流項には3次精度の風上近似（Kawamura&Kuwahara1984, Kawamura+1986）を適用してきたが，過去にそうしたように，こういった"それらしい"流れ場における近似スキームの影響を調べてみる．

移流スキームとして，

1次精度風上差分
QUICK（Leaonard1979）
KK（Kawamura&Kuwahara1984, Kawamura+1986）

を比較する（村上+1988がインテンシブに纏めている，こちらも詳しい）．移流項以外は全て同一の近似スキーム（2次精度中心差分）とする．過去のテストでは，以下のような収束の傾向を得ている．ところで，付録で議論しているように，私は，QUICKは2次精度，KKは3次精度であると考えている．そしてこの数値実験は，私の考えをサポートしている（この数値実験だけで結論を出すのは良くないだろうが，考えを黙っておくのも良くないだろうし，間違っていたら理解をsolidにする良い機会だ）．2次精度の方法には，Lax-Wendroff（Lax&Wendroff1960）も有名だが，空間離散化の影響にのみ集中したい（Lax-Wendroffは，時間・空間両者について2次精度）ので，これら3つを比較する．

昔，各スキームの収束性を調べたとき，以下のような結果を得た．

Gauss波の移流	各スキームの収束性

$\mathrm{Re} \simeq 200, \ 2,000, \ \beta = 0.05, \ h = 5 \times 10^{-3} \ [ \mathrm{m} ], \ \tau = 5 \times 10^{-4} \ [ \mathrm{s} ], \ C_s = 0.1$とする．上下端は非すべりし，parabolicな流入流速を与える．

$\mathrm{scheme}$	$\mathrm{Re} \simeq 200$	$\mathrm{Re} \simeq 2,000$
$\mathrm{1st}$
$\mathrm{QUICK}$
$\mathrm{KK}$

1次精度風上差分の結果を見ると，$\mathrm{Re} \simeq 200, \ 2,000$が殆ど同じ流れ場になっている．これは，1次精度スキームの数値粘性が，$\mathrm{Re} \simeq 2,000$の比較的激しい流れを酷く鈍らせていることに起因すると考える．実際，同スキームは，

\begin{align}
    u \frac{d \phi}{d x}
    \sim \mathcal{A}_{(1)} \left( u, \phi \right)
        - \mathcal{D}_{(2)}^{(2)} \left( |u| \frac{h}{2} , \phi \right)
\end{align}

と書くことができる．ここで，$\mathcal{A} (a, b)$は移流の近似（$a$の速度で$b$が移流する），$\mathcal{D} (a, b)$は拡散の近似（$a$の速さで$b$が拡散する）とした．下付き文字は近似精度，上付き文字は微分階数である．上式から，分子粘性$\nu$と同じオーダー（2階）の数値粘性$\bar{\nu}$が，$h (= 5 \times 10^{-3} \ [ \mathrm{m} ])$のオーダーで効いている．今回の問題設定では，

$\mathrm{Re} = 200$　: $\nu = 2.5 \times 10^{-4} \ [ \mathrm{m^2 / s} ] < \bar{\nu} \sim \mathcal{O} (10^{-3})$
$\mathrm{Re} = 2,000$: $\nu = 2.5 \times 10^{-5} \ [ \mathrm{m^2 / s} ] < \bar{\nu} \sim \mathcal{O} (10^{-3})$

であり，数値粘性が実質的な粘性（＝分子粘性＋数値粘性）の効果を殆ど支配してしまっている．これが，1次精度風上差分では，何れの$\mathrm{Re}$でもほぼ同様の計算結果となった理由であろうと考える．こちらでは，低次精度スキームは「計算が安定・収束性が良い」という旨の記述があるが，私は，分子粘性と同じオーダーの数値粘性が混入すると，本来とは異なる状態への収束が起きると考えている．即ち，（分子粘性と同程度の）数値粘性にバイアスされた状態への収束性が良い，という方が正しいと思う（accuracy vs. precision）．但し，数値粘性にバイアスされた解は，解空間の中で，ある尺度で，ある程度まで，本来の状態に近いと考えられるだろうから，手早く大雑把に現象を追跡したいときには有用だろうと考える（しかし，やはりある程度の範疇まで，，ここまで現象と乖離してしまうとね，，）．

QUICKは，1次精度風上差分よりシャープに渦を解像しているが，何だか違和感が拭えない．良い言葉が見当たらないが，$\mathrm{Re} \simeq 200$の時点で既に，流れが不自然に裂けている（渦が分散・振動する様子）．しかし，KKによる$\mathrm{Re} \simeq 200$の結果は自然に見える．或いは，私がいつもKKの結果を見るから，バイアスされているだけかもしれない．どちらが正しいかは分からない，ということを強調しつつ，相違点を調べてみる．両者は，形式的に，

\begin{alignat}{2}
    &\mathrm{QUICK:} \quad
    &&u \frac{d \phi}{d x}
    \sim \mathcal{A}_{(2)} \left( u, \phi \right)
        + \mathcal{D}_{(2)}^{(4)} \left( |u| \frac{h^3}{16}, \phi \right) \\
    &\mathrm{KK:} \quad
    &&u \frac{d \phi}{d x}
    \sim \mathcal{A}_{(3)} \left( u, \phi \right)
        + \mathcal{D}_{(2)}^{(4)} \left( |u| \frac{h^3}{4}, \phi \right)
\end{alignat}

と書ける．相違点は，移流項の近似精度と数値粘性の大きさである．ううむ，移流項の近似精度が2次か3次かで結果にここまで大きな影響は出るだろうか？拡散項と圧力勾配項を2次精度で近似しているので，スキーム全体で空間に関して2次精度なのだ（勿論，inertiaが強い問題であると認めた上で）．数値粘性の強さが差異の原因ではなかろうかと当たりを付けてみる．両者をたすきがけして，

\begin{alignat}{2}
    &\mathrm{QUICK-KK:} \quad
    &&u \frac{d \phi}{d x}
    \sim \mathcal{A}_{(2)} \left( u, \phi \right)
        + \mathcal{D}_{(2)}^{(4)} \left( |u| \frac{h^3}{4}, \phi \right) \\
    &\mathrm{KK-QUICK:} \quad
    &&u \frac{d \phi}{d x}
    \sim \mathcal{A}_{(3)} \left( u, \phi \right)
        + \mathcal{D}_{(2)}^{(4)} \left( |u| \frac{h^3}{16}, \phi \right) 
\end{alignat}

を用意し⁴，再計算してみる．私の予想は合っているだろうか．

$\mathrm{scheme}$	$\mathrm{Re} \simeq 200$	$\mathrm{Re} \simeq 2,000$
$\mathrm{QUICK-KK}$
$\mathrm{KK-QUICK}$

観測結果として，QUICK-KK（移流項の2次近似 + 大きい4階数値粘性）は，流れ場が安定しており，KKと酷似した結果を与えている．一方，KK-QUICK（移流項の3次近似 + 小さい4階数値粘性）は，流れ場が不安定であり，QUICKと非常によく似ている．以上から，恐らく私の予想は正しかったと考える．村上+1988によると，河村は数値粘性にさらに調整可能な係数を乗じたことがあるそうだが（原論文を見つけられなかった），数値粘性を極限まで小さくしてゼロにしてしまったりすると，スキーム自体が不安定化するから（実際，QUICKはKKと同じ階数の数値粘性であるのに，非常に不安定），そこまでしなくても良いと思う（そのままのKKで十分だろう，という意味）．

本項での結論は，QUICKは数値粘性が小さすぎてスキームが不安定化しつつある，ということである．KKは安定であるが，しかし本当に現象を正しく追従できているかについては，未だ議論の余地がある．次項で定量的な評価を試みる．

Strouhal数

Strouhal数

\begin{align}
    \mathrm{St}
    &= \frac{f D}{U}
\end{align}

は，流れ場の特徴を示す無次元量であり，前項で触れた移流スキームの定量的な比較の指標となり得る（$f$は現象の周波数$[ \mathrm{ / s} ]$，$D$は特徴長さ$[ \mathrm{m} ]$，$U$は特徴流速$[ \mathrm{m / s} ]$）．

blockage ratioを$\beta = 0.05$に固定し，境界条件の議論から上下端に周期条件を課した流れ場に対して，空間解像度を少しずつ上昇させつつ，Strouhal数を計測する．定量的な評価のため，先行研究（Norberg1993, Sharma&Eswaran2004, Sen+2011, Mashhadi+2021）を参照し，手元の計算と比較する．また，Williamson1988は，circular cylinderのlaminar vortex shedding regime（Williamson1988は$60 \le \mathrm{Re} \le 180$程度を調べている）において，Strouhal数とReynolds数の間には，

\begin{align}
    \mathrm{St}
    &= \frac{1}{\mathrm{Re}} \left( c_0 + c_1 \mathrm{Re} + c_2 \mathrm{Re}^2 \right) \\
    \mathrm{where} \
    \left( c_0, c_1, c_2 \right) 
    &= \left( -3.3265, 0.1816, 1.6 \times 10^{-4} \right)
\end{align}

なる関係があると報告している．障害物の形状が違うのでsolidには言い難いが，参考程度に併記した．$\mathrm{St}$の計測に当たり，周波数$f \ [ \mathrm{/s} ]$をどうやって測るのが良いか悩んだが，簡単に，$(x, y) = (1.5, 0.5)$における鉛直流速を指標にした．つまり，1秒間に何回鉛直流速が振動するか，である．しかし，これはあまり良い方針じゃないかもしれない．

$\mathrm{QUICK}$	$\mathrm{KK}$

先ず，circular cylinderに関する議論であるWilliamson1988と，square cylinderの議論であるNorberg1993, Sharma&Eswaran2004, Sen+2011, Mashhadi+2021の間には大きな乖離がある．Williamson1988の式を$y$軸方向に平行移動させた破線は，square cylinderに対してもtendencyを捉えているかもしれないが，多分，もう全く違う現象なのではないかと思ってしまう．

続いて，手元の計算結果について．$\mathrm{Re} = 100$では，KKは数値解らしく参照点に近づいている．QUICKはそこそこである（中間の解像度が爆発しているのは変だが）．$\mathrm{Re} = 150$も調べたが，QUICKもKKも参照点からは激しく乖離し，もう自信が無い（数値解らしい収束の傾向も無い）．QUICKは散付きが顕著だが，KKもそう優れているとは言い辛い結果になってしまった．色々，当初の目論見から大きく外れ，満足できる結果を得なかった．残念．

== 雑談ここから ==
こちらによると，円柱の場合，流れ場のStrouhal数は，$200 \le \mathrm{Re} \le 200,000$の広い範囲で，$\mathrm{St} \approx 0.2$程度に留まることが知られているらしい．一方，こちらによると，角柱では$\mathrm{St} \approx 0.16$程度らしい．激しい流れになると，3次元的な渦の広がりが卓越するので，現在のコードではこれだけ広い範囲のReynolds数を調べることは叶わない．計算コスト的にも，かなり難しい．

また，こちらにもあるように，例えばバスなどの長方形らしい形状をした車両（bluff body）は，そうでない車両より，背圧による抗力が大きい．同サイトには，角柱を含め，様々な形状に於ける抗力係数も載っている．
== 雑談ここまで ==

終わりに

本当は，blockage ratioについて調べたり，壁の影響を遠くにしてやったりなど，やりたいことは沢山あるのだが，本稿だけで今月のファイルアップロード上限に達してしまい，もうこれ以上の数値実験を報告できないので，一旦はここまでとする．目に楽しい結果たちであった．

(2024/11追記)
数値実験の章で議論が広がりすぎてしまった感は否めない．が，色々とやっている内に気になることが増えてしまい，気になったことはやってみないと気が済まない性分なので，まあそういう訳で，全部上手く行った訳ではないが，やりたかったことは一通り済ませたつもりだ．

直交グリッドしか考慮していないので，DFG benchmarkなどと比較できないのが心残りだ（空間を変換すれば可能かもしれないが，そういうことは，現状，考えていない）．

Appendix

無次元量

本文に登場した無次元量について記しておく．

Reynolds数

慣性力と粘性力の比．

\begin{align}
    \mathrm{Re}
    &= \frac{U D}{\nu}
\end{align}

ここで，$\nu$は動粘度$\nu \ [ \mathrm{m^2 / s} ]$，$U$は特徴流速$[ \mathrm{m / s} ]$，$D$は特徴長さ$[ \mathrm{m} ]$である．

Strouhal数

時間変化によって生じる力と流れの慣性力の比．

\begin{align}
    \mathrm{St}
    &= \frac{f D}{U}
\end{align}

ここで，$f$は現象の周波数$[ \mathrm{ / s} ]$，$D$は特徴長さ$[ \mathrm{m} ]$，$U$は特徴流速$[ \mathrm{m / s} ]$である．

Expansion ratio

backstep流れの話題では，次のexpansion ratioがしばしば話題に上がる．

\begin{align}
    \alpha
    &= \frac{L}{H}
\end{align}

ここで，$L$は流入口の高さ$[ \mathrm{m} ]$，$H$は階段の高さ$[ \mathrm{m} ]$である．

Blockage ratio

von Kármán渦列の話題では，次のblockage ratioが良く話題に上がる．

\begin{align}
    \beta
    &= \frac{D}{H}
\end{align}

ここで，$D$は障害物の特徴長さ$[ \mathrm{m} ]$，$H$は流入口の高さ$[ \mathrm{m} ]$である．

解析解

簡単な問題には，その解析的な解が分かる場合がある．

Hagen–Poiseuille流れ

流れが1次元的であるとして，

\begin{align}
    \frac{\partial u}{\partial x}
    &= 0 \\
    \frac{\partial u}{\partial t}
    + u \frac{\partial u}{\partial x}
    &= - \frac{1}{\rho} \frac{\partial p}{\partial x}
    + \nu \left(
        \frac{\partial^2 u}{\partial x^2} + \frac{\partial^2 u}{\partial y^2}
    \right)
\end{align}

を考える．十分な時間が経ち，fully-developedであるとすると，

\begin{align}
    \frac{\partial u}{\partial t}
    = 0, \quad
    \frac{\partial u}{\partial x}
    = 0, \quad
    \frac{\partial p}{\partial x}
    = \mathrm{(const.)}
\end{align}

より，

\begin{align}
    0
    &= - \frac{1}{\rho} \frac{\partial p}{\partial x}
    + \nu \frac{\partial^2 u}{\partial y^2} \\
    \therefore
    u
    &= \frac{1}{2 \rho \nu} \frac{\partial p}{\partial x} y^2 + c_1 y + c_2 \\
\end{align}

境界条件に，$u (0) = u (L_y) = 0$を与えると，

\begin{align}
    c_1
    = - \frac{1}{2 \rho \nu} \frac{\partial p}{\partial x} L_y, \quad
    c_2
    = 0 
\end{align}

より，

\begin{align}
    u
    &= \frac{1}{2 \rho \nu} \frac{\partial p}{\partial x} y \left( y - L_y \right) \\
\end{align}

が解析解である．これは綺麗なparabolaなので，適当に微分したり平方完成したりして，

\begin{align}
    u \left( y = \frac{1}{2} L_y \right)
    &= - \frac{1}{8 \rho \nu} \frac{\partial p}{\partial x} L_y^2 \\
\end{align}

が最大流速であると分かる．

Couette–Poiseuille流れ

Hagen–Poiseuilleのときと同様の手続きにより，

\begin{align}
    u
    &= \frac{1}{2 \rho \nu} \frac{\partial p}{\partial x} y^2 + c_1 y + c_2 \\
\end{align}

という一般解に，先程とは異なる境界条件（$u (0) = 0, u (L_y) = U^{\prime}$）を与え，

\begin{align}
    c_1
    = - \frac{1}{2 \rho \nu} \frac{\partial p}{\partial x} L_y + \frac{U^{\prime}}{L_y}, \quad
    c_2
    = 0 
\end{align}

より，

\begin{align}
    u
    &= \frac{1}{2 \rho \nu} \frac{\partial p}{\partial x} y \left( y - L_y \right) 
        + \frac{U^{\prime}}{L} y \\
\end{align}

が解析解である．歪んではいるが，2次式なので，$\partial u / \partial y = 0$に最大流速がある．

移流項の有限差分近似

何度でも勉強し直そう．過去の議論との差分として，打ち切りによる誤差項をより明示的に残しながら考える．

1次精度風上差分

風上方向からの移流を意識して，

\begin{align}
    u \frac{d \phi}{d x}
    =
    \begin{cases}
        u \left(\phi_{i}   - \phi_{i-1}\right) / h + \mathcal{O} \left[ h \right] & \left( u \ge 0 \right) \\
        u \left(\phi_{i+1} - \phi_{i}  \right) / h + \mathcal{O} \left[ h \right] & \left( u \lt 0 \right)
    \end{cases}
\end{align}

を基礎に考える．移流の向きによる場合分けを，

\begin{align}
    u 
    = \frac{u + \lvert u \rvert}{2} + \frac{u - \lvert u \rvert}{2}
\end{align}

の関係に吸収させれば，

\begin{align}
    u \frac{d \phi}{d x}
    &= \frac{u + \lvert u \rvert}{2} \frac{\phi_{i} - \phi_{i-1}}{h}
        + \frac{u - \lvert u \rvert}{2} \frac{\phi_{i+1} - \phi_{i}}{h}
        + \mathcal{O} \left[ h \right] \\
    &= u \frac{\phi_{i+1} - \phi_{i-1}}{h}
        - |u| \frac{h}{2} \frac{\phi_{i+1} - 2 \phi_{i} + \phi_{i-1}}{h^2}
        + \mathcal{O} \left[ h \right] \\
\end{align}

を得る．スキーム全体で1次精度である．

QUICK

$\phi_{i \pm 1/2}$を$\phi_{i}$周りでTaylor級数展開し，差分を取る．

\begin{align}
    \phi_{i \pm 1/2}
    &= \phi_{i} 
        \pm \frac{h}{2} \left. \frac{d \phi}{d x} \right|_{i}
        + \frac{1}{2!} \left( \frac{h}{2} \right)^2 \left. \frac{d^2 \phi}{d x^2} \right|_{i}
        + \mathcal{O} \left[ h^3 \right] \\
    \phi_{i + 1/2} - \phi_{i - 1/2}
    &= 2 \cdot \frac{h}{2} \left. \frac{d \phi}{d x} \right|_{i}
        + \mathcal{O} \left[ h^3 \right] \\
    \therefore
    \left. \frac{d \phi}{d x} \right|_{i}
    &= \frac{\phi_{i + 1/2} - \phi_{i - 1/2}}{h}
        + \mathcal{O} \left[ h^2 \right] \\
\end{align}

続いて，$\phi_{i + 1/2}$を$\phi_{i}$周りで，$\phi_{i - 1/2}$を$\phi_{i-1}$周りで，それぞれTaylor級数展開する（途中で2次精度の差分近似を含むので，実際の打ち切り誤差は2次まで含んでしまうだろうと思っている，どちらにせよ，後から2次の打ち切り誤差が支配するが）．

\begin{align}
    \phi_{i + 1/2}
    &= \phi_{i} 
        + \frac{h}{2} \left. \frac{d \phi}{d x} \right|_{i}
        + \frac{1}{2} \left( \frac{h}{2} \right)^2 \left. \frac{d^2 \phi}{d x^2} \right|_{i}
        + \mathcal{O} \left[ h^3 \right] \\
    &= \phi_{i} 
        + \frac{h}{2} \frac{\phi_{i+1} - \phi_{i-1}}{2 h}
        + \frac{1}{2} \left( \frac{h}{2} \right)^2 \frac{\phi_{i+1} - 2 \phi_{i} + \phi_{i-1}}{h^2}
        + \mathcal{O} \left[ h^2 \right] \\
    &= \frac{3 \phi_{i+1} + 6 \phi_{i} - \phi_{i-1}}{8}
        + \mathcal{O} \left[ h^2 \right] \\
    \phi_{i - 1/2}
    &= \phi_{i-1} 
        + \frac{h}{2} \left. \frac{d \phi}{d x} \right|_{i-1}
        + \frac{1}{2} \left( \frac{h}{2} \right)^2 \left. \frac{d^2 \phi}{d x^2} \right|_{i-1}
        + \mathcal{O} \left[ h^3 \right] \\
    &= \frac{3 \phi_{i} + 6 \phi_{i-1} - \phi_{i-2}}{8}
        + \mathcal{O} \left[ h^2 \right] \\
\end{align}

したがって，元の式は，

\begin{align}
    \left. \frac{d \phi}{d x} \right|_{i}
    &= \frac{\phi_{i + 1/2} - \phi_{i - 1/2}}{h}
        + \mathcal{O} \left[ h^2 \right] \\
    &= \frac{3 \phi_{i+1} + 3 \phi_{i} - 7 \phi_{i-1} + \phi_{i-2}}{8 h}
        + \mathcal{O} \left[ h^2 \right] \\
\end{align}

即ち，

\begin{align}
    u \left. \frac{d \phi}{d x} \right|_{i}
    &= \frac{u + |u|}{2} \frac{3 \phi_{i+1} + 3 \phi_{i} - 7 \phi_{i-1} + \phi_{i-2}}{8 h} \\
    &+ \frac{u - |u|}{2}\frac{- 3 \phi_{i-1} - 3 \phi_{i} + 7 \phi_{i+1} - \phi_{i+2}}{8 h}
     + \mathcal{O} \left[ h^2 \right] \\
    &= u \frac{ - \phi_{i+2} + 10 \phi_{i+1} - 10 \phi_{i-1} + \phi_{i-2}}{16 h} \\
    &+ |u| \frac{h^3}{16} \frac{\phi_{i+2} - 4 \phi_{i+1} + 6 \phi_{i} - 4 \phi_{i-1} + \phi_{i-2}}{h^4}
     + \mathcal{O} \left[ h^2 \right] \\
\end{align}

を得る．Leaonardは，QUICKは3次精度であると主張しているが（Leaonard1979），そしてそのように言われたりもするようだが（e.g., Versteeg&Malalasekera1995），実際は2次精度しかないのではないか，と思う．数値実験からも，QUICKは2次程度の速度で収束することが確認できる．

KK

Kawamura+1986が詳しい．

$\phi_{i-1}$と$\phi_{i-2}$を$\phi_{i}$周りでTaylor級数展開する．

\begin{align}
    \phi_{i-1} 
    &= \phi_{i} 
    - h \left. \frac{d \phi}{d x} \right|_{i}
    + \frac{h^{2}}{2!} \left. \frac{d^2 \phi}{d x^2} \right|_{i}
    + \mathcal{O} \left[ h^3 \right] \\
    \phi_{i-2} 
    &= \phi_{i} 
    - 2 h \left. \frac{d \phi}{d x} \right|_{i}
    + \frac{(2 h)^{2}}{2!} \left. \frac{d^2 \phi}{d x^2} \right|_{i}
    + \mathcal{O} \left[ h^3 \right] \\
\end{align}

両者の差を取り，2階微分を除去し，1階微分の差分近似を考える．

\begin{align}
    4 \phi_{i-1} - \phi_{i-2}
    &= 3 \phi_{i} 
    - 2 h \left. \frac{d \phi}{d x} \right|_{i}
    + \mathcal{O} \left[ h^3 \right] \\
    \therefore 
    \left. \frac{d \phi}{d x} \right|_{i}
    &= \frac{3 \phi_{i} - 4 \phi_{i-1} + \phi_{i-2}}{2 h}
    + \mathcal{O} \left[ h^2 \right] \\
\end{align}

より，

\begin{align}
    u \left. \frac{d \phi}{d x} \right|_{i}
    &= \frac{u + |u|}{2} \frac{3 \phi_{i} - 4 \phi_{i-1} + \phi_{i-2}}{2 h} \\
    &+ \frac{u - |u|}{2} \frac{- 3 \phi_{i} + 4 \phi_{i+1} - \phi_{i+2}}{2 h}
    + \mathcal{O} \left[ h^2 \right] \\
    &= u \frac{- \phi_{i+2} + 4 \phi_{i+1} - 4 \phi_{i-1} + \phi_{i-2}}{4 h} \\
    &+ |u| \frac{h^3}{4} \frac{\phi_{i+2} - 4 \phi_{i+1} + 6 \phi_{i} - 4 \phi_{i-1} + \phi_{i-2}}{h^4}
    + \mathcal{O} \left[ h^2 \right] \\
\end{align}

ここで，第1項を見つめる．先ずは，$\phi_{i \pm 1}$と$\phi_{i \pm 2}$のTaylor級数展開を記す．

\begin{align}
    \phi_{i \pm 1}
    &= \phi_{i} 
        \pm h \left. \frac{d \phi}{d x} \right|_{i}
        +   \frac{h^2}{2!} \left. \frac{d^2 \phi}{d x^2} \right|_{i}
        \pm \frac{h^3}{3!} \left. \frac{d^3 \phi}{d x^3} \right|_{i}
        + \mathcal{O} \left[ h^4 \right] \\
    \phi_{i \pm 2}
    &= \phi_{i} 
        \pm 2h \left. \frac{d \phi}{d x} \right|_{i}
        +   \frac{(2h)^2}{2!} \left. \frac{d^2 \phi}{d x^2} \right|_{i}
        \pm \frac{(2h)^3}{3!} \left. \frac{d^3 \phi}{d x^3} \right|_{i}
        + \mathcal{O} \left[ h^4 \right] \\
\end{align}

差分から，

\begin{align}
    \phi_{i+1} - \phi_{i-1}
    &= 2h \left. \frac{d \phi}{d x} \right|_{i}
        + 2 \cdot \frac{h^3}{3!} \left. \frac{d^3 \phi}{d x^3} \right|_{i}
        + \mathcal{O} \left[ h^4 \right] \\
    \phi_{i+2} - \phi_{i-2}
    &= 2 \cdot 2h \left. \frac{d \phi}{d x} \right|_{i}
        + 2 \cdot \frac{(2h)^3}{3!} \left. \frac{d^3 \phi}{d x^3} \right|_{i}
        + \mathcal{O} \left[ h^4 \right] \\
\end{align}

より，

\begin{align}
    \left. \frac{d \phi}{d x} \right|_{i}
    &= \frac{- \phi_{i+2} + 4 \phi_{i+1} - 4 \phi_{i-1} + \phi_{i-2}}{4h}
        + \mathcal{O} \left[ h^2 \right] \\
\end{align}

であるから，上式の第1項は，2次程度の精度であることが分かる．一方，係数を上手く選ぶと3次の項を取り除くことができ，これにより3次の打ち切り誤差を除去することができる．

\begin{align}
    \left. \frac{d \phi}{d x} \right|_{i}
    &= \frac{- \phi_{i+2} + 8 \phi_{i+1} - 8 \phi_{i-1} + \phi_{i-2}}{12h}
        + \mathcal{O} \left[ h^3 \right] \\
\end{align}

したがって，先程の，移流項の2次精度近似を，新たに得た（係数を上手く選んだ）3次精度の近似で置き換えてしまおう，というのがKKスキームである．

\begin{align}
    u \left. \frac{d \phi}{d x} \right|_{i}
    &= u \frac{- \phi_{i+2} + 8 \phi_{i+1} - 8 \phi_{i-1} + \phi_{i-2}}{12 h} \\
    &+ |u| \frac{h^3}{4} \frac{\phi_{i+2} - 4 \phi_{i+1} + 6 \phi_{i} - 4 \phi_{i-1} + \phi_{i-2}}{h^4}
    + \mathcal{O} \left[ h^3 \right] \\
\end{align}

重みと除算を少し書き換えただけだが，3次精度まで上昇している．Gauss波を移流させる数値実験からも，KKは3次の速度があることが確認できる．

したがって，ここでの$\mathrm{Re}$は，$\mathrm{Re}{\mathrm{max}}$の意味で定義している．$\mathrm{Re}{\mathrm{mean}}$の意味で解釈するには，最大流速と平均流速の関係から，それぞれ2/3倍して考えれば良い． ↩
陽解法は，こういうとき（失敗するとき）に，ショックが大きい．簡単で速い陽解法は素敵だが，安定性の高い陰解法も魅力的に感じる． ↩
対称な境界というものもあるらしい．cf. こちら．が，正直に白状すると，良く分かっていない． ↩
村上+1988によると，LeonardによるUTOPIA（Uniformly Third Order Polynomial Interpolation Algorithm）と呼ばれるスキームがあるらしく（昔受講した大学院の講義でも，UTOPIAスキームに少しだけ触れていたが，この原論文を見つけられなかった），それは，ここでのKK-QUICKによく似ている． ↩