Navier-Stokes Simulation 渦度方程式

Navier-Stokes Simulation 基礎編 の続きです。

参考文献：Fluid Simulation for Computer Graphics (CG分野の教科書)

(追記) 渦度方程式を再帰的な積分方程式に見立てて、モンテカルロ法で確率的シミュレーションするアルゴリズムをCG分野の学会で論文発表しました
→ A Monte Carlo Method for Fluid Simulation (ACM SIGGRAPH Asia 2022)

Overview

非圧縮・非粘性流体のNavier-Stokes方程式:

\begin{eqnarray}
\frac{\partial\vec{u}}{\partial t} + (\vec{u}\cdot\nabla)\vec{u} 
= -\frac{1}{\rho}\nabla P + \vec{F}_{ext} \tag{1}
\end{eqnarray} 

を渦度$\vec\omega$

\begin{eqnarray}
\vec \omega &\equiv& \nabla \times \vec{u}
= \begin{pmatrix} \frac{\partial}{\partial x} \\ \frac{\partial}{\partial y} \\ \frac{\partial}{\partial z} \end{pmatrix} 
\times \begin{pmatrix} u_x \\ u_y \\ u_z \end{pmatrix}
= \begin{pmatrix} 
\frac{\partial u_z}{\partial y}-\frac{\partial u_y}{\partial z} \\
\frac{\partial u_x}{\partial z}-\frac{\partial u_z}{\partial x} \\
\frac{\partial u_y}{\partial x}-\frac{\partial u_x}{\partial y} \end{pmatrix} 
\end{eqnarray} 

によって書き換える。後述の通り、渦度方程式

\begin{eqnarray}
\frac{\partial\vec{\omega}}{\partial t} + (\vec u \cdot \nabla) \vec\omega = (\vec \omega \cdot \nabla) \vec u + \nabla \times \vec F_{ext} \tag{2}
\end{eqnarray}

のようになる。時刻$t$における速度場$\vec u_t$から定義通り渦度$\vec \omega_t$が計算出来、上の渦度方程式から$\vec \omega_{t+1}$が求まる。

\begin{eqnarray}
\vec u_t \,\,\,\,\,\, \rightarrow \,\,\,\,\,\, \vec \omega_t &\,\, \rightarrow \,\,& \vec\omega_{t+1} \,\,\,\,\,\, \rightarrow \,\,\,\,\,\, \vec u_{t+1} \tag{3} \\
\mathrm{Definition} \,\,\,\,\, &\mathrm{Vorticity\,eq.}& \,\,\,\,\,\,\,\,\,\,\,\,\, \mathrm{(later)}
\end{eqnarray}

そして最後の$\vec \omega_{t+1}$から$\vec u_{t+1}$の計算は少々複雑である(後述)。

Naiveに離散化して線形solverに投げることは難しく、渦度の定義にできるだけ近いような速度場$\vec u$を近似的な答えとして求めるような問題設定にすると、次のようなPoisson方程式が得られる。

\begin{eqnarray}
-\nabla^2 \vec u = \nabla \times \vec \omega \tag{4
}
\end{eqnarray} 

渦度方程式の導出 (Vorticity equation)

式$(1)$の両辺にrot演算子を作用させて、

\begin{eqnarray}
\nabla \times \left( \frac{\partial\vec{u}}{\partial t} + (\vec{u}\cdot\nabla)\vec{u} \right)
&=& \nabla \times \left( -\frac{1}{\rho}\nabla P + \vec{F}_{ext} \right) \\
\frac{\partial\vec{\omega}}{\partial t} + \nabla \times (\vec{u}\cdot\nabla)\vec{u} 
&=& -\frac{1}{\rho} \left(\nabla \times \nabla P\right) + \nabla \times \vec F_{ext} \tag{5


}
\end{eqnarray} 

         ・第$3$項について、

\begin{eqnarray}
\nabla \times \nabla &=& 
\begin{pmatrix} \frac{\partial}{\partial x} \\ \frac{\partial}{\partial y} \\ \frac{\partial}{\partial z} \end{pmatrix} \times
\begin{pmatrix} \frac{\partial}{\partial x} \\ \frac{\partial}{\partial y} \\ \frac{\partial}{\partial z} \end{pmatrix}
= \begin{pmatrix} \frac{\partial^2}{\partial y \partial z}-\frac{\partial^2}{\partial z \partial y} \\ \frac{\partial^2}{\partial x \partial z}-\frac{\partial^2}{\partial z \partial x} \\ \frac{\partial^2}{\partial x \partial y}-\frac{\partial^2}{\partial y \partial x} \end{pmatrix}
= \vec 0
\end{eqnarray} 

         ・第$2$項について、

\begin{eqnarray}
\nabla \times (\vec{u}\cdot\nabla)\vec{u} &=&
\begin{pmatrix} \frac{\partial}{\partial x} \\ \frac{\partial}{\partial y} \\ \frac{\partial}{\partial z} \end{pmatrix} \times
\left( u_x\frac{\partial}{\partial x} + u_y\frac{\partial}{\partial y} + u_z\frac{\partial}{\partial z} \right)
\begin{pmatrix} u_x \\ u_y \\ u_z \end{pmatrix} \\
&=& \cdots \\
&=& (\vec u \cdot \nabla) \vec\omega - (\vec \omega \cdot \nabla) \vec u
\end{eqnarray} 

(計算の詳細)

--> 第$1$成分のみについて計算するとして、

\begin{eqnarray}
&& \frac{\partial}{\partial y} \left( u_x\frac{\partial u_z}{\partial x} + u_y\frac{\partial u_z}{\partial y} + u_z\frac{\partial u_z}{\partial z} \right) - \frac{\partial}{\partial z} \left( u_x\frac{\partial u_y}{\partial x} + u_y\frac{\partial u_y}{\partial y} + u_z\frac{\partial u_y}{\partial z} \right) \\
&=& 
\frac{\partial u_x}{\partial y} \frac{\partial u_z}{\partial x} + u_x \frac{\partial^2 u_z}{\partial y \partial x} 
+ \frac{\partial u_y}{\partial y} \frac{\partial u_z}{\partial y} + u_y \frac{\partial^2 u_z}{\partial^2 y}
+ \frac{\partial u_z}{\partial y} \frac{\partial u_z}{\partial z} + u_z \frac{\partial^2 u_z}{\partial y \partial z} \\
&-& \frac{\partial u_x}{\partial z} \frac{\partial u_y}{\partial x} - u_x \frac{\partial^2 u_y}{\partial z \partial x} 
- \frac{\partial u_y}{\partial z} \frac{\partial u_y}{\partial y} - u_y \frac{\partial^2 u_y}{\partial z \partial y}
- \frac{\partial u_z}{\partial z} \frac{\partial u_y}{\partial z} - u_z \frac{\partial^2 u_y}{\partial^2 z}
\end{eqnarray}

-----> 第$2,4,6,8,10,12$項について、

\begin{eqnarray}
&& u_x \left( \frac{\partial^2 u_z}{\partial x \partial y}-\frac{\partial^2 u_y}{\partial x \partial z} \right) 
+ u_y \left( \frac{\partial^2 u_z}{\partial y \partial y}-\frac{\partial^2 u_y}{\partial y \partial z} \right)
+ u_z \left( \frac{\partial^2 u_z}{\partial y \partial z}-\frac{\partial^2 u_y}{\partial z \partial z} \right) \\
&=& u_x \frac{\partial}{\partial x} \left( \frac{\partial u_z}{\partial y}-\frac{\partial u_y}{\partial z} \right)
+ u_y \frac{\partial}{\partial y} \left( \frac{\partial u_z}{\partial y}-\frac{\partial u_y}{\partial z} \right)
+ u_z \frac{\partial}{\partial z} \left( \frac{\partial u_z}{\partial y}-\frac{\partial u_y}{\partial z} \right) \\
&=& \left( u_x \frac{\partial}{\partial x} + u_y \frac{\partial}{\partial y} + u_z \frac{\partial}{\partial z} \right)
\left( \frac{\partial u_z}{\partial y}-\frac{\partial u_y}{\partial z} \right) \\
&=& (\vec u \cdot \nabla) \omega_x
\end{eqnarray}

-----> 第$3,5,9,11$項について、

\begin{eqnarray}
&& \frac{\partial u_y}{\partial y} \left( \frac{\partial u_z}{\partial y} - \frac{\partial u_y}{\partial z} \right)
+ \frac{\partial u_z}{\partial z} \left( \frac{\partial u_z}{\partial y} - \frac{\partial u_y}{\partial z} \right) \\
&=& \omega_x \left( \frac{\partial u_y}{\partial y}+\frac{\partial u_z}{\partial z} \right) \\
&=& \omega_x \left( \nabla \vec u - \frac{\partial u_x}{\partial x} \right) \\
&=& \omega_x \nabla \vec u - \omega_x \frac{\partial u_x}{\partial x} \\
&=& - \omega_x \frac{\partial u_x}{\partial x}
\end{eqnarray}

-----> これに第$1,7$項を加えて、

\begin{eqnarray}
&& \frac{\partial u_x}{\partial y} \frac{\partial u_z}{\partial x} - \frac{\partial u_x}{\partial z} \frac{\partial u_y}{\partial x} - \omega_x \frac{\partial u_x}{\partial x} \\
&=& \left( \frac{\partial u_x}{\partial z}-\omega_y \right) \frac{\partial u_x}{\partial y} - \left( \frac{\partial u_x}{\partial y}+\omega_z \right) \frac{\partial u_x}{\partial z} - \omega_x \frac{\partial u_x}{\partial x} \\
&=& - \left( \omega_x\frac{\partial}{\partial x}+\omega_y\frac{\partial}{\partial y}+\omega_z\frac{\partial}{\partial z} \right) u_x \\
&=& - (\vec \omega \cdot \nabla) u_x
\end{eqnarray}

--> 第$1$成分は結局次のようになる:

\begin{eqnarray}
(\vec u \cdot \nabla) \omega_x - (\vec \omega \cdot \nabla) u_x 
\end{eqnarray}

--> 第$2$成分は第$1$成分を$x \rightarrow y, y \rightarrow z, z \rightarrow x$とズラしたものになる。よって、

\begin{eqnarray}
(\vec u \cdot \nabla) \omega_y - (\vec \omega \cdot \nabla) u_y 
\end{eqnarray}

--> 第$3$成分も同様になるので、結局

\begin{eqnarray}
(\vec u \cdot \nabla) \vec\omega - (\vec \omega \cdot \nabla) \vec u
\end{eqnarray}

---

すべてまとめると、Navier-Stokes方程式は

\begin{eqnarray}
\frac{\partial\vec{\omega}}{\partial t} + (\vec u \cdot \nabla) \vec\omega - (\vec \omega \cdot \nabla) \vec u = \nabla \times \vec F_{ext}
\end{eqnarray}

整理して、

\begin{eqnarray}
\frac{\partial\vec{\omega}}{\partial t} + (\vec u \cdot \nabla) \vec\omega = (\vec \omega \cdot \nabla) \vec u + \nabla \times \vec F_{ext} \tag{2}
\end{eqnarray}

左辺は渦度$\vec \omega$についての移流項になっている。物質微分$\frac{D}{Dt}\equiv\frac{\partial}{\partial t} + (\vec u \cdot \nabla)$を定義すると、

\begin{eqnarray}
\frac{D\vec\omega}{Dt} = (\vec \omega \cdot \nabla) \vec u + \nabla \times \vec F_{ext}
\end{eqnarray}

Velocity from Vorticity

渦度方程式$(2)$に基づいて時刻$t$における$\vec u_t$,$\vec\omega_t$から次の時刻ステップにおける渦度場$\vec\omega_{t+1}$を求めたとする。この$\vec\omega_{t+1}$から$\vec u_{t+1}$を求める方法を考える。

$\vec\omega$と$\vec u$の関係式は

\begin{eqnarray}
\vec \omega = \nabla \times \vec u = 
\begin{pmatrix}
\frac{\partial u_z}{\partial y}-\frac{\partial u_y}{\partial z} \\ 
\frac{\partial u_x}{\partial z}-\frac{\partial u_z}{\partial x} \\ 
\frac{\partial u_y}{\partial x}-\frac{\partial u_x}{\partial y}
\end{pmatrix}
\end{eqnarray} 

である。

・Naiveに格子法で空間方向に離散化して線形solverにする場合
Poisson方程式の場合ように疎な正定値対称行列にならず、前処理付共役勾配法(PCG)で解くことが出来ない。(？？？)

・また、求めた渦度場$\vec\omega$が渦度としての性質を満たす、つまり$\nabla \times \vec u=\vec\omega$を満たすような速度場$\vec u$が存在するとは限らない。例えば、$\nabla\cdot(\nabla\times\vec u)=0$の関係から、$\vec\omega$は$\nabla\cdot\vec\omega=0$のように発散$0$の条件を満たさなければならない(逆に発散$0$であれば$\vec u$の存在は保証され、grad$\Phi$の不定性を除いて一意に定まる:参考pdf)。しかし、格子法での空間分割による数値誤差等により条件が満たされる保証は無い。

         (離散化しなければOK? $\rightarrow$ 式$(4)$の両辺に$\nabla$を作用させると、

\begin{eqnarray}
\nabla \left(\frac{\partial\vec{\omega}}{\partial t}\right) + \nabla\cdot\left(\nabla \times (\vec{u}\cdot\nabla)\vec{u}\right) 
&=& \nabla\cdot\left(\nabla \times \vec F_{ext}\right)
\end{eqnarray} 

         第一項以外は$\nabla\cdot(\nabla\times\vec A)=0$になるので

\begin{eqnarray}
\frac{\partial(\nabla\cdot\vec{\omega})}{\partial t} &=& 0
\end{eqnarray} 

         よって計算誤差が無ければ$\nabla\cdot\vec\omega=0$は保持される。)

Poisson方程式

そこで、$\nabla \times \vec u$が$\vec\omega$になるべく近くなるような速度場$\vec u$を近似的な答えとみなすことにする。

問題設定としては、

\begin{eqnarray}
\vec u = \arg\min_{\vec u} \int \|\nabla\times\vec u - \vec \omega\|^2 dV \tag{6}
\end{eqnarray} 

となる$\vec u$を求めればよい。(明確な境界条件は無く、積分領域は全空間・無限遠では$\omega \rightarrow 0$であると仮定する。もし境界条件がある場合、その内部にのみ渦度があり外側では0であるので、上の仮定に基づく式変形は可能である。境界条件により発生する速度場への影響はその後に不定性を用いて修正することが出来る。 )

理想的な速度場$\vec u$について、ベクトル場$\vec p$(全領域で微分可能・境界付近で$\vec0$に収束する)を定義し、次のような関数を定義する。

\begin{eqnarray}
g(s) = \int\int\int \|\nabla\times (\vec u+s\vec p) - \vec \omega\|^2 dV
\end{eqnarray} 

$g(s)$は$s=0$において最小値を取る。つまり、微分値が$g(s=0)'=0$となるはずである。

\begin{eqnarray}
g(s) &=& \int \left(\nabla\times (\vec u+s\vec p) - \vec \omega\right) \left(\nabla\times (\vec u+s\vec p) - \vec \omega\right) dV \\
&=& \int \left[(\nabla\times\vec u - \vec \omega) + \nabla\times s \vec p\right] \left[(\nabla\times\vec u - \vec \omega) + \nabla\times s \vec p\right] dV \\
&=& \int \left[ \|\nabla\times\vec u - \vec \omega\|^2 + 2(\nabla\times\vec u - \vec \omega)(\nabla\times\vec p)s + \|\nabla\times\vec p\|^2 s^s \right] dV \\

g'(0) &=& \int (\nabla\times\vec u - \vec \omega)(\nabla\times\vec p) \, dV = 0
\end{eqnarray} 

ここで部分分数分解 (境界付近で$\vec p=\vec 0$になる仮定を使う)

\begin{eqnarray}
\int \vec A \cdot (\nabla \times \vec p) \,dV 
&=& \int A_x \left(\frac{\partial p_z}{\partial y} - \frac{\partial p_y}{\partial z}\right) 
+ A_y \left(\frac{\partial p_x}{\partial z} - \frac{\partial p_z}{\partial x}\right) 
+ A_z \left(\frac{\partial p_y}{\partial x} - \frac{\partial p_x}{\partial y}\right) \,dV \\
&=& \left( \left[ A_x p_z\right] - \int \frac{\partial A_x}{\partial y} p_z \,dV \right) - \left( \left[ A_x p_y\right] - \int \frac{\partial A_x}{\partial z} p_y \,dV \right) + \cdots \\
&=& - \int \left( \frac{\partial A_x}{\partial y} p_z - \frac{\partial A_x}{\partial z} p_y \right)
+ \left( \frac{\partial A_y}{\partial z} p_x - \frac{\partial A_y}{\partial x} p_z \right)
+ \left( \frac{\partial A_z}{\partial x} p_y - \frac{\partial A_z}{\partial y} p_x \right) \,dV \\
&=& - \int \left( \frac{\partial A_y}{\partial z} -\frac{\partial A_z}{\partial y} \right) p_x + \cdots \,dV \\
&=& \int (\nabla \times \vec A) \cdot \vec p \, dV
\end{eqnarray} 

により、

\begin{eqnarray}
g'(0) = \int \left( \nabla \times (\nabla \times \vec u - \vec \omega) \right) \cdot \vec p \, dV = 0
\end{eqnarray} 

ベクトル場$\vec p ({\bf r})$は任意に変化させられるので、結局

\begin{eqnarray}
\nabla \times (\nabla \times \vec u - \vec \omega) = \vec 0
\end{eqnarray} 

が全空間において成立することになる。書き直して、

\begin{eqnarray}
\nabla \times \nabla \times \vec u = \nabla \times \vec \omega \tag{7}
\end{eqnarray} 

左辺の$\nabla \times \nabla$をもう少し分解してみる。

\begin{eqnarray}
\nabla \times \nabla \times \vec u &=& 
\begin{pmatrix} \partial_x \\ \partial_y \\ \partial_z \end{pmatrix} \times
\begin{pmatrix} \partial_y u_z - \partial_z u_y \\ \partial_z u_x - \partial_x u_z \\ \partial_x u_y - \partial_y u_x \end{pmatrix} \\
&=& \begin{pmatrix} \partial_y\partial_x u_y - \partial_y^2u_x -\partial_z^2u_x + \partial_z\partial_xu_z \\ \cdots \\ \cdots \end{pmatrix} \\
&=& \begin{pmatrix} \partial_x (\partial_xu_x+\partial_yu_y+\partial_zu_z) - (\partial_x^2+\partial_y^2+\partial_z^2)u_x \\ \cdots \\ \cdots \end{pmatrix} \\
&=& \nabla\cdot(\nabla\vec u) - \nabla^2 \vec u
\end{eqnarray} 

第一項は速度場が理想的に発散$0$であれば消滅する(発散$0$になるとは限らない、後述の追加議論参照)。よって、

\begin{eqnarray}
-\nabla^2 \vec u = \nabla \times \vec \omega \tag{4}
\end{eqnarray} 

これは各成分ごとのPoisson方程式になっている。(式$(4)$を解いた結果の速度場$\vec u$が発散$0$を満たしていないと仮定が破綻する...)

*上記教科書通りの説明では速度場の発散$0$条件や不定性についての議論があいまいなままであるため詳細を追加する。ます、問題設定となる式$(6)$の時点で、速度場$\vec u({\bf r})$には任意のポテンシャル関数$\Phi({\bf r})$の勾配項だけ任意性が許されている。

\begin{eqnarray}
\vec u({\bf r}) \rightarrow \vec u({\bf r}) + \nabla \Phi({\bf r})
\end{eqnarray} 

その帰結となる式$(7)$の解についても同様の不定性が存在する。

よって式$(7)$を$\nabla \vec u=0$という仮定のもとに導出したPoisson方程式$(4)$を解いて得られる$\vec u$について、本当に$\nabla \vec u=0$を満たすようポテンシャル勾配項で調整する操作を行えば辻褄が合うことになる(もちろん勾配項を追加することでPoisson方程式$(4)$は満たされなくなるが、元の式$(7)$は保たれているのでOK)。つまり、

\begin{eqnarray}
\nabla \cdot \left( \vec u({\bf r}) + \nabla \Phi({\bf r}) \right) &=& 0 \\
\rightarrow - \nabla^2 \Phi &=& \nabla \vec u
\end{eqnarray} 

結局Poisson方程式$(4)$を解いて得られた速度場$\vec u$に対して、さらにこのPoisson方程式を解いてポテンシャル$\Phi$を決定する。$\vec u + \nabla \Phi$が発散$0$と渦度の定義を両方満たす解になっている。(下の球対称ポテンシャルを用いて発散0の解を書き下せるので、この議論は結局不必要となる。ただし、境界条件による影響を考える時に同様の不定性が使える。)

Poisson方程式の解

Poisson方程式$(4)$と発散$0$条件を同時に満たすような解$\vec u({\bf r})$は書き下す事ができる。

まず、次のような球対称関数$\Phi({\bf r})$を定義する。

\begin{eqnarray}
\Phi({\bf r}) &=& - \frac{1}{4\pi\|{\bf r}\|} \\
\nabla \Phi({\bf r}) &=& \frac{\bf r}{4\pi\|{\bf r}\|^3} \,\,\,({\bf r}\neq\vec 0) \\
\nabla^2 \Phi({\bf r}) &=& \frac{1}{4\pi} \nabla \begin{pmatrix} \frac{x}{(x^2+y^2+z^2)^{3/2}} \\ \frac{y}{(x^2+y^2+z^2)^{3/2}} \\ \frac{z}{(x^2+y^2+z^2)^{3/2}} \end{pmatrix} \\
&=& \frac{1}{4\pi} \left( \frac{(x^2+y^2+z^2)^{3/2} - x\frac{3}{2}(x^2+y^2+z^2)^{1/2} 2x}{(x^2+y^2+z^2)^3} + \cdots \right) \\
&=& \frac{1}{4\pi} \left( \frac{x^2+y^2+z^2 - 3x^2}{(x^2+y^2+z^2)^{5/2}} + \cdots \right) \\
&=& 0 \,\,\, ({\bf r \neq \vec0})\\
\int\int\int_{全空間} \nabla^2 \Phi({\bf r}) \,dV 
&=& \int\int\int_{単位球} \nabla^2 \Phi({\bf r}) \,dV \,\,\,(\nabla^2\Phiは{\bf r}=\vec0以外で0,積分領域を縮小) \\
&=& \int\int_{\partial S} \nabla\Phi({\bf r}) \, \vec n \,dS \,\,\,(ガウスの発散定理) \\
&=& \frac{1}{4\pi} \int\int_{\partial S} \frac{{\bf r}\cdot{\bf r}}{\|{\bf r}\|^3} \,dS \\
&=& 1
\end{eqnarray} 

$\nabla^2 \Phi({\bf r})$は原点で$\infty$かつ他で0、積分値が1なので定義通りのデルタ関数になっている。

(2Dの場合)

\begin{eqnarray}
\Phi({\bf r}) &=& \frac{1}{2\pi} \log\|{\bf r}\| \\
\nabla \Phi({\bf r}) &=& \frac{\bf r}{2\pi\|{\bf r}\|^2} \,\,\,({\bf r}\neq\vec 0) \\
\nabla^2 \Phi({\bf r}) &=& \frac{1}{2\pi} \nabla \begin{pmatrix} \frac{x}{(x^2+y^2)} \\ \frac{y}{(x^2+y^2)}  \end{pmatrix} \\
&=& 0 \,\,\, ({\bf r \neq \vec0})\\
\int\int\int_{全空間} \nabla^2 \Phi({\bf r}) \,dV 
&=& \int\int\int_{単位円} \nabla^2 \Phi({\bf r}) \,dV \,\,\,(\nabla^2\Phiは{\bf r}=\vec0以外で0,積分領域を縮小) \\
&=& \int\int_{\partial S} \nabla\Phi({\bf r}) \, \vec n \,dS \,\,\,(ガウスの発散定理) \\
&=& \frac{1}{2\pi} \int\int_{\partial S} \frac{{\bf r}\cdot{\bf r}}{\|{\bf r}\|^2} \,dS \\
&=& 1
\end{eqnarray} 

これを用いて、Poisson方程式の解は、

\begin{eqnarray}
\vec u({\bf r}) = -\int \Phi({\bf r}-{\bf p}) \nabla_{\bf p} \times \vec \omega({\bf p}) \, d{\bf p} \tag{7}
\end{eqnarray} 

のように書ける。

         (証明1) 両辺に$\nabla^2$を作用させると、式$(4)$と一致する:

\begin{eqnarray}
-\nabla^2_{\bf r} \vec u({\bf r}) &=& \nabla^2_{\bf r} \int \Phi({\bf r}-{\bf p}) \nabla_{\bf p} \times \vec \omega({\bf p}) \, d{\bf p} \\
&=& \int \delta({\bf r}-{\bf p}) \nabla_{\bf p} \times \vec \omega({\bf p}) \, d{\bf p} \\
&=& \nabla_{\bf r} \times \vec \omega({\bf r})
\end{eqnarray} 

         (証明2) 両辺に$\nabla$を作用させると、$\nabla\cdot\vec u=0$の発散$0$条件が満たされる:

\begin{eqnarray}
\nabla_{\bf r} \vec u({\bf r}) &=& -\nabla_{\bf r} \int \Phi({\bf r}-{\bf p}) \left(\nabla_{\bf p} \times \vec \omega({\bf p})\right) \, d{\bf p} \\
&=& - \int \left(\nabla_{\bf r}\Phi({\bf r}-{\bf p})\right) \left(\nabla_{\bf p} \times \vec \omega({\bf p})\right) \, d{\bf p} \\
&=& \int \left(\nabla_{\bf p}\Phi({\bf r}-{\bf p})\right) \left(\nabla_{\bf p} \times \vec \omega({\bf p})\right) \, d{\bf p} \\
&=& \left[ -\frac{1}{4\pi \|{\bf r}-{\bf p}\|} \left(\nabla_{\bf p} \times \vec \omega({\bf p})\right)\right]_{{\bf p}=-\infty}^{\infty} 
- \int \Phi({\bf r}-{\bf p}) \nabla_{\bf p}\cdot\left(\nabla_{\bf p} \times \vec \omega({\bf p})\right) \, d{\bf p} \\
&=& 0
\end{eqnarray} 

                      (無限遠において渦度場が0であるとして第一項は0、第二項は$\nabla\cdot\left(\nabla \times \vec \omega\right)=0$により$0$となる)

式$(7)$をさらに変形して、

\begin{eqnarray}
\vec u({\bf r}) &=& -\int \Phi({\bf r}-{\bf p}) \nabla_{\bf p} \times \vec \omega({\bf p}) \, d{\bf p} \\
&=& \int \nabla_{\bf p} \Phi({\bf r}-{\bf p}) \times \vec \omega({\bf p}) \, d{\bf p} \,\,\, (部分分数分解) \\
&=& -\int \nabla_{\bf r} \Phi({\bf r}-{\bf p}) \times \vec \omega({\bf p}) \, d{\bf p} \,\,\,\,\,\,(\leftarrow \mathrm{Biot–Savartの法則}) \tag{8} \\
&=& -\int \nabla_{\bf r} \times \left( \Phi({\bf r}-{\bf p}) \vec \omega({\bf p}) \right) \, d{\bf p} \\
&=& -\nabla_{\bf r} \times \int \Phi({\bf r}-{\bf p}) \vec \omega({\bf p}) \, d{\bf p} \tag{9}
\end{eqnarray} 

Streamfunction

ここで、Streamfunction $\vec \Psi({\bf r})$:

\begin{eqnarray}
\vec \Psi({\bf r}) \equiv -\int \Phi({\bf r}-{\bf p}) \vec \omega({\bf p}) \, d{\bf p}
\end{eqnarray} 

を定義すると、式$(9)$は

\begin{equation}
    \left\{
    \begin{aligned}
          \vec u({\bf r}) &= \nabla \times \vec \Psi({\bf r}) \\
          \nabla^2 \vec \Psi({\bf r}) &= - \vec \omega({\bf r})
    \end{aligned}
    \right.
\end{equation}

渦度場$\vec \omega$から下のポアソン方程式よりStreamfunction $\vec \Psi({\bf r})$が求まり、上式から速度場$\vec u({\bf r})$が求まる。($\vec u({\bf r})$がベクトルポテンシャルのrotで表されているので、$\nabla \cdot \vec u=0$は自動的に満たされている)

Bior-Savartの法則

・これまで速度場$\vec u({\bf r})$に関する次の条件

\begin{equation}
    \left\{
    \begin{aligned}
          \nabla \cdot \vec u &= 0 \\
          \nabla \times \vec u &= \vec \omega
    \end{aligned}
    \right.
\end{equation}

から、無限遠における境界条件$\vec \omega({\bf r}\rightarrow \infty)=0$の下で式$(8)$:

\begin{eqnarray}
\vec u({\bf r}) 
&=& -\int \nabla_{\bf r} \Phi({\bf r}-{\bf p}) \times \vec \omega({\bf p}) \, d{\bf p} \tag{8} \\
&=& -\int \frac{\bf r-p}{4\pi\|{\bf r-p}\|^3} \times \vec \omega({\bf p}) \, d{\bf p} \\
&=& \frac{1}{4\pi} \int \vec \omega({\bf p}) \times \frac{\bf r-p}{\|{\bf r-p}\|^3} \, d{\bf p} \\
\end{eqnarray} 

が導出された(3次元)。これは電磁気学におけるBiot-Savartの法則と全く同じ式である。

・電磁気学では、磁場$\vec B({\bf r})$に関するMaxwell方程式

\begin{equation}
    \left\{
    \begin{aligned}
          \nabla \cdot \vec B &= 0 \\
          \nabla \times \vec B &= \mu_0 \vec j
    \end{aligned}
    \right.
\end{equation}

から次のBior-Savartの法則:

\begin{eqnarray}
\vec B({\bf r}) 
&=& \frac{\mu_0}{4\pi} \int \vec j({\bf p}) \times \frac{\bf r-p}{\|{\bf r-p}\|^3} \, d{\bf p} \\
\end{eqnarray} 

が導かれる(参考pdf)。