Helmholtzの定理とChorinの射影法

初めに

Helmholtzの定理とChorinの射影法（Chorin1968）には関連がありそうだ，ということが，射影法のWikiやChorinの原論文からから窺える．しかし，私にはあれだけの説明では良く分からない．もう少しだけで良いから両者の繋がりを明確に解釈したい．

著者は工学部の出身で，数学をきちんと学んできた訳ではない．本稿は，間違いを含む可能性が大いにあるし，寧ろ間違っているだろう，というくらいには注意深く読んでいただき，その場合はご指摘いただきたい．

A mathematical introduction to fluid mechanics (Chorin&Marsden1993) などを参考にしながら記す．

(2024/05追記)
儀我1982をきちんと勉強すれば，もう少し詳しいことが分かりそうだ，というメモ．

Helmholtzの定理

Helmholtzの定理，Helmholtzの分解の定理，Helmholtz(-Hodge)分解などと呼ばれる．「ベクトル解析の基本定理」と呼ばれるほどに重要な定理らしいが，著者が大学で履修したベクトル解析の講義では学ばなかった（工学部だったからだろうか）．

任意のベクトル場$\mathbf{v}$が，回転無し場$\mathbf{v}_L$と発散無し場$\mathbf{v}_T$に分解できるという定理である．

\begin{align}
    \mathbf{v}
    &= \mathbf{v}_L + \mathbf{v}_T \\
    &= \nabla \phi + \nabla \times \mathbf{\Psi} \\
    &= \mathrm{grad} \ \phi + \mathrm{curl} \ \mathbf{\Psi} \\
\end{align}

ここで，$\mathbf{v}_L$は回転無し場（longitudinal / irrotational / curl-free field），$\mathbf{v}_T$は発散無し場（transverse / solenoidal / divergence-free field）である．また，$\phi$は$\mathbf{v}_L$のスカラーポテンシャル，$\mathbf{\Psi}$は$\mathbf{v}_T$のベクトルポテンシャルである．

任意のベクトル場$\mathbf{v}$が，，と書いたが，一般化しすぎると話が難しいので（かつ，それは本稿のスコープではないので），特に非圧縮性Navier-Stokes方程式との関係を探ることにする（この定理は電磁気学で良く顔を出すらしいが，そちらはあまり考えない）．

一意性？

分解が存在することは取り敢えず認めて「回転無し場と発散無し場は一意に定まるか？」というのが疑問として湧いてくる．

回答（解答ではない）は「適切な境界条件の下で，かつ，一定の自由度の下で，一意に定まる」だと思う．

ここからは，非圧縮性Navier-Stokes方程式との関連付けを明確化するために，先に記したHelmholtzの定理を改めて記す．

$\Omega \subset \mathbb{R}^{d} \ (d = 2, 3)$上のベクトル場$\mathbf{v}$は，回転無し場$\mathbf{v}_L = \nabla p$と発散無し場$\mathbf{v}_T = \mathbf{u}$に分解できる（明らかに非圧縮性NS式の変数を意識した命名だが，未だ圧力場とか速度場とか呼ぶことは避け，取り敢えずスカラー場$p$，ベクトル場$\mathbf{u}$と呼ぶこととする）．

\begin{align}
    \mathbf{v}
    &= \mathbf{v}_L + \mathbf{v}_T \\
    &= \nabla p + \mathbf{u} \\
\end{align}

ここで，$p$は$\mathbf{v}_L$を与えるスカラーポテンシャルであり，$\Omega$上で定義されるものとする．かつ，境界条件として，

\mathbf{n} \cdot \mathbf{u} = 0

を課す．ここで，$\mathbf{n}$は境界$\partial \Omega$上で定義される外向き単位法線ベクトルとする．

では，一意性を考える．スカラー場$p$とベクトル場$\mathbf{u}$に対して，

\begin{align}
    \nabla \cdot \left( p \mathbf{u} \right)
    &= \nabla p \cdot \mathbf{u}
    + p \nabla \cdot \mathbf{u} \\
    &= \nabla p \cdot \mathbf{u}
\end{align}

$\Omega$上で積分すると，

\begin{align}
    \int_{\Omega} \nabla \cdot \left( p \mathbf{u} \right) dV
    &= \int_{\Omega} \nabla p \cdot \mathbf{u} dV \\
    &= \int_{\partial \Omega} p \mathbf{u} \cdot \mathbf{n} dS \\
    &= 0
\end{align}

すなわち，回転無し場$\nabla p$と発散無し場$\mathbf{u}$は（積分の上で）直交していると云える．

さて，元のベクトル場$\mathbf{v}$を分解するベクトル場として，$( \nabla p, \mathbf{u} )$と$( \nabla p^{\prime}, \mathbf{u}^{\prime} )$があるとすると，

\begin{align}
    \mathbf{0}
    &= \mathbf{v} - \mathbf{v} \\
    &= \left( \nabla p + \mathbf{u} \right)
     - \left( \nabla p^{\prime} + \mathbf{u}^{\prime} \right) \\
    &= \left( \nabla \left( p - p^{\prime} \right) \right)
     + \left( \mathbf{u} - \mathbf{u}^{\prime} \right) \\
\end{align}

であるから，あとは，$\nabla \left( p - p^{\prime} \right) = \mathbf{0}$，$\mathbf{u} - \mathbf{u}^{\prime} = \mathbf{0}$のどちらか一方が示せれば良さそうだ．

上式で$\mathbf{u} - \mathbf{u}^{\prime}$との内積を取り，$\Omega$上で積分すると，

\begin{align}
    0
    &= \int_{\Omega} \nabla \left( p - p^{\prime} \right) \cdot \left( \mathbf{u} - \mathbf{u}^{\prime} \right) d\Omega
    + \int_{\Omega} \| \mathbf{u} - \mathbf{u}^{\prime} \|^2 d\Omega \\
    &= \int_{\Omega} \| \mathbf{u} - \mathbf{u}^{\prime} \|^2 d\Omega \\
\end{align}

より，$\mathbf{u} = \mathbf{u}^{\prime}$である．よって，分解の式から$\nabla p = \nabla p^{\prime}$を得る．勾配が一致するので，適当に平行移動する範囲内で一意に定まると云える（$\nabla p^{\prime} = \nabla (p + c) = \nabla p \ \text{for some} \ c \in \mathbb{R}$）．

Neumann問題

Helmholtzの定理$\mathbf{v} = \nabla p + \mathbf{u}$に対して，発散を取る．

\begin{align}
    \nabla \cdot \mathbf{v}
    &= \nabla^2 p + \nabla \cdot \mathbf{u} \\
    &= \nabla^2 p
\end{align}

また，$\partial \Omega$上の単位法線ベクトル$\mathbf{n}$との内積を取る．

\begin{align}
    \mathbf{n} \cdot \mathbf{v}
    &= \mathbf{n} \cdot \nabla p + \mathbf{n} \cdot \mathbf{u} \\
    &= \mathbf{n} \cdot \nabla p
\end{align}

上記の2式から，スカラー場$p$に関するNeumann問題を得る．

\begin{align}
    \nabla^2 p &= \nabla \cdot \mathbf{v} \ \text{in} \ \Omega \\
    \mathbf{n} \cdot \nabla p &= \mathbf{n} \cdot \mathbf{v} \ \text{on} \ \partial \Omega \\
\end{align}

この問題に対する解$p$は，先ほどと同様に，定数を適当に動かす範囲内で一意に定まることが知られている（らしい, cf. Courant&Hilbert1953, Courant&Hilbert1961）．すなわち，幾つかの解（$p_1, p_2; p_1 \neq p_2$）が見つかっても，$\nabla p_1 = \nabla p_2$なので，Helmholtz分解$\mathbf{v} = \nabla p + \mathbf{u}$には特に影響を与えない．

解として得られる$p$は，

\begin{align}
    \mathbf{u} &= \mathbf{v} - \nabla p \\
    \nabla \cdot \mathbf{u}
    &= \nabla \cdot \mathbf{v} - \nabla^2 p = 0 \\
\end{align}

より，$\mathbf{u}$に要求される発散無しの条件を満たす．また，上式から，$\mathbf{v}$と$\nabla p$が既知のとき，$\nabla p$は$\mathbf{v}$を発散無しの条件を満たす空間に引き戻す作用を持っていると解釈できる．

Chorinの射影法

まずは，非圧縮性のNavier-Stokes方程式を記す．

\begin{align}
    \nabla \cdot \mathbf{u}
    &= 0 \\
    \frac{\partial \mathbf{u}}{\partial t}
    + \left( \mathbf{u} \cdot \nabla \right) \mathbf{u}
    &= - \nabla p
       + \frac{1}{\mathrm{Re}} \nabla^2 \mathbf{u}
       + \mathbf{f} \\
\end{align}

ここで，$p$は圧力場，$\mathbf{u}$は速度場，$\mathbf{f}$は外力，$\mathrm{Re}$はReynolds数である．

Chorin（Chorin1968）は，運動量保存則を次のように並び替え，Helmholtzの定理を適用する．

\begin{align}
    \nabla p
    + \frac{\partial \mathbf{u}}{\partial t}
    &= - \left( \mathbf{u} \cdot \nabla \right) \mathbf{u}
       + \frac{1}{\mathrm{Re}} \nabla^2 \mathbf{u}
       + \mathbf{f} \\
    \mathbf{w}_L + \mathbf{w}_T = \mathbf{w}
    :&= \nabla p + \frac{\partial \mathbf{u}}{\partial t} \\
    \mathbf{v}
    :&= - \left( \mathbf{u} \cdot \nabla \right) \mathbf{u}
        + \frac{1}{\mathrm{Re}} \nabla^2 \mathbf{u}
        + \mathbf{f} \\    
\end{align}

ここで，$\nabla p$をスカラー場$p$が定める回転無し場と解釈する．$\partial \mathbf{u} / \partial t$に対して発散を作用させれば，$\nabla \cdot (\partial \mathbf{u} / \partial t) = \partial / \partial t \ (\nabla \cdot \mathbf{u}) = 0$であるから，これは発散無し場である．

作用素$\mathcal{P} (\cdot)$を定義する．$\mathcal{P} (\cdot)$は，適当なベクトル場に作用して，その発散無し場を取り出すような作用素とする．

\begin{align}
    \mathcal{P} ( \mathbf{w} ) &= \mathbf{w}_T \\
    \therefore \mathbf{w} &= \nabla p + \mathcal{P} ( \mathbf{w} )
\end{align}

発散無し場$\mathbf{w}_T$に対して$\mathcal{P} (\cdot)$を作用させれば，$\mathbf{w}_T$から発散無し場を取り出すのだから，当然そのものが返される．

\mathcal{P} ( \mathbf{w}_T ) = \mathbf{w}_T

回転無し場$\nabla p = \mathbf{w}_L$は，元のベクトル場$\mathbf{w}$から発散無し場$\mathbf{w}_T$を取り除いだものであるから，$\nabla p$に$\mathcal{P} (\cdot)$を作用させるとゼロである（ゼロベクトルの場が取り出される）．

\mathcal{P} ( \nabla p ) = \mathbf{0}

これを，非圧縮性NS式の運動量保存則に作用させる．

\begin{align}
    \mathcal{P}
    \left(
        \nabla p
        + \frac{\partial \mathbf{u}}{\partial t}
    \right)
    &= \mathcal{P} \left( \mathbf{v} \right) \\
    \therefore
    \frac{\partial \mathbf{u}}{\partial t}
    &= \mathcal{P} \left( \mathbf{v} \right) \\
\end{align}

すなわち，非圧縮性の速度場$\mathbf{u}$の時間発展は，移流・拡散・外力の効果を発散無し場に射影したもので記述される，という主張に見える．

ここまでの式展開は，凡そChorin&Marsden1993の主張通りに記したつもりである．この後，圧力は右辺（$\mathbf{v} = - \left( \mathbf{u} \cdot \nabla \right) \mathbf{u} + \nabla^2 \mathbf{u} / \mathrm{Re} + \mathbf{f}$）の勾配として回復できると記してあるが，どういう意味か良く理解できない．

The pressure can then be recovered as the gradient part of $- \mathbf{u} \cdot \nabla \mathbf{u} + 1/\mathrm{R} \Delta \mathbf{u}$.
Chorin&Marsden1993, ch1.3 The Navier-Stokes Equations, p39
（Chorin&Marsden1993では，Laplace作用素を$\Delta$と，Reynolds数を$\mathrm{R}$と書くことにしている．また，外力はゼロとしている．）

やはり自分なりに解釈したいので，もう少し考えてみる（つまり，以降の議論は私が勝手に云っているだけである）．

そもそも，$\partial \mathbf{u} / \partial t = \mathcal{P} \left( \mathbf{v} \right)$の右辺はどうやって計算するのだろう？Euler陽解法とすれば，少なくとも，$\mathcal{P} ( \cdot )$の中身は計算できそうだ．

\frac{\mathbf{u}^{(n+1)} - \mathbf{u}^{(n)}}{\Delta t}
= \mathcal{P} \left( \mathbf{v}^{(n)} \right) 

ここで，$\mathbf{v}^{(n)} := - \left( \mathbf{u}^{(n)} \cdot \nabla \right) \mathbf{u}^{(n)} + \nabla^2 \mathbf{u}^{(n)} / \mathrm{Re} + \mathbf{f}^{(n)}$とした．えいやっと，$\mathcal{P} ( \cdot )$を恒等写像$\mathcal{I} ( \cdot )$に置き換えてみる．

\frac{\mathbf{u}^{(n+1)} - \mathbf{u}^{(n)}}{\Delta t}
= \mathcal{P} \left( \mathbf{v}^{(n)} \right) 
\simeq \mathcal{I} \left( \mathbf{v}^{(n)} \right) 

これで，Euler陽解法として時間積分が実行できる．めでたしめでたし．．で良いのだろうか？$\mathcal{P} \left( \mathbf{v}^{(n)} \right)$は，$\mathbf{v}^{(n)}$の発散無し場を取り出してくるものだったから，$\left( \mathbf{u}^{(n+1)} - \mathbf{u}^{(n)} \right) / \Delta t = \mathcal{P} \left( \mathbf{v}^{(n)} \right)$の時間積分は$\mathbf{u}^{(n)}$を発散無し場$\mathbf{u}^{(n+1)}$に安全に写してくれそうな気がする．

今回は，$\mathcal{P} \left( \mathbf{v}^{(n)} \right) \simeq \mathcal{I} \left( \mathbf{v}^{(n)} \right)$という乱暴な近似（のようなもの）を経由した．この時間積分は，$\mathbf{u}^{(n)}$を発散無し場$\mathbf{u}^{(n+1)}$に写してくれるだろうか？

回答は，否である（数値的にも，中間速度場の発散は$\Omega$上の至る所でゼロ，とはなっていない¹）．すなわち，上式の時間積分で得られるのは，何らかの仮の速度場に過ぎない．本稿ではこれを，中間速度場$\mathbf{u}^{(*)}$と呼ぶこととし，$\mathbf{u}^{(n)}$，$\mathbf{u}^{(n+1)}$と区別する．

\begin{align}
    \frac{\mathbf{u}^{(n+1)} - \mathbf{u}^{(n)}}{\Delta t}
    &= \mathcal{P} \left( \mathbf{v}^{(n)} \right) \\
    \rightarrow
    \frac{\mathbf{u}^{(*)} - \mathbf{u}^{(n)}}{\Delta t}
    &= \mathcal{I} \left( \mathbf{v}^{(n)} \right)
\end{align}

本当は発散無し場のみを作用させるはずだったのに（$\mathcal{P} \left( \mathbf{v}^{(n)} \right)$），回転無し場まで作用させてしまった（$\mathcal{I} \left( \mathbf{v}^{(n)} \right)$）．つまり，回転無し場$\nabla p$を引けば，発散無し場に戻ってくれるはずだろう．先に記したように，「$\mathbf{v}$と$\nabla p$が既知のとき，$\nabla p$は$\mathbf{v}$を発散無しの条件を満たす空間に引き戻す作用を持っている」ような気がしてくる（Kim&Moin1985にも同様の記述がある）．

In application of this method to the Navier-Stokes equations, one can interpret the role of pressure in the momentum equations as a projection operator which projects an arbitrary vector field into a divergence-free vector field.
Kim&Moin1985: Application of a fractional-step method to incompressible Navier-Stokes equations

\begin{align}
    \mathrm{(div \text{-} free \ field)}
    =:
    \Delta \mathbf{u}
    + \frac{\mathbf{u}^{(*)} - \mathbf{u}^{(n)}}{\Delta t}
    &= \mathcal{I} \left( \mathbf{v}^{(n)} \right)
    - \nabla p
\end{align}

ここで，$\Delta \mathbf{u}$は速度場の修正量である（$\Delta$はLaplace作用素ではなく，何らかの差分）．これは，本来の非圧縮性のNS式に一致し，また，左辺は発散無し場に戻るはずであるから，$\Delta \mathbf{u} = (\mathbf{u}^{(n+1)} - \mathbf{u}^{(*)}) / \Delta t$であろう．結局，元の非圧縮性NS式は，

\begin{align}
    \frac{\mathbf{u}^{(*)} - \mathbf{u}^{(n)}}{\Delta t}
    &= \mathcal{I} \left( \mathbf{v}^{(n)} \right) \\
    \frac{\mathbf{u}^{(n+1)} - \mathbf{u}^{(*)}}{\Delta t}
    &= - \nabla p^{(n+1)}
\end{align}

の2式に分解されていたのだと解釈する．なお，$(n+1)$ステップ目で発散無し場に引き戻すため，$\nabla p$には$\nabla p^{(n+1)}$を選んだ．また，第2式の発散から圧力Poisson方程式を得るから，その解が為すベクトル場で以て速度場の修正を行うものと解釈する．

上記を図化すると，以下のようになるだろうか．

終わりに

何となく尤もらしい，恰も正しそうなことを書いた．

間違いなく間違いがあるだろう（普通に怖い）から，公開しようか迷ったが，詳しい方からのご指摘を期待し，まずは公開することにした．自分でも，理解を深化して加筆したい．

また，Chorinの射影法は結局のところEuler陽解法なので，時間方向に1次の精度しか無い．時間積分法についても，周辺の方法（例えば，Crank–NicolsonやRunge–Kutta）と合わせて纏めてみたい．

Appendix A: 回転無し，発散無しか？

計算の復習を兼ね，$\mathbf{v} = \mathbf{v}_L + \mathbf{v}_T = \nabla \phi + \nabla \times \mathbf{\Psi}$と分解される$\mathbf{v}_L$，$\mathbf{v}_T$がそれぞれが回転無し，発散無しになっているのか確かめてみる（ベクトル解析の知識から自明な性質であるが，，）．

以降，$\partial_{j} \left( \cdot \right) = \partial \left( \cdot \right) / \partial x_{j}$とした．

回転無しか？

$\mathbf{v}_L$の回転を計算してみると，

\begin{align}
    \mathrm{curl} \left( \mathbf{v}_L \right)
    &= \mathrm{curl} \left( \mathrm{grad} \ \phi \right) \\
    &= \nabla \times \left( \nabla \phi \right) \\
    &= \nabla \times \left(
        \begin{array}{c}
            \partial_{1} \phi \\
            \partial_{2} \phi \\
            \partial_{3} \phi
        \end{array}
    \right) \\
    &= \left(
        \begin{array}{c}
            \partial_{2} \left( \partial_{3} \phi \right) - \partial_{3} \left( \partial_{2} \phi \right) \\
            \partial_{3} \left( \partial_{1} \phi \right) - \partial_{1} \left( \partial_{3} \phi \right) \\
            \partial_{1} \left( \partial_{2} \phi \right) - \partial_{2} \left( \partial_{1} \phi \right) \\
        \end{array}
    \right) \\
    &= \mathbf{0} \\
\end{align}

となり，確かに回転無しである．

発散無しか？

$\mathbf{v}_T$の発散を計算してみると，

\begin{align}
    \mathrm{div} \left( \mathbf{v}_T \right)
    &= \mathrm{div} \left( \mathrm{curl} \ \mathbf{\Psi} \right) \\
    &= \nabla \cdot \left( \nabla \times \mathbf{\Psi} \right) \\
    &= \nabla \cdot \left(
        \begin{array}{c}
            \partial_{2} \Psi_{3} - \partial_{3} \Psi_{2} \\
            \partial_{3} \Psi_{1} - \partial_{1} \Psi_{3} \\
            \partial_{1} \Psi_{2} - \partial_{2} \Psi_{1}
        \end{array}
    \right) \\
    &= \partial_{1} \left( \partial_{2} \Psi_{3} - \partial_{3} \Psi_{2} \right) 
     + \partial_{2} \left( \partial_{3} \Psi_{1} - \partial_{1} \Psi_{3} \right) 
     + \partial_{3} \left( \partial_{1} \Psi_{2} - \partial_{2} \Psi_{1} \right) \\
    &= 0 \\
\end{align}

となり，確かに発散無しである．

1. 例えば，Cavity流れの差分法シミュレーションでは，これが中間速度場の発散であり，これが次の速度場の発散である． ↩