[Python] オイラー的渦識別法

はじめに

流体力学において「渦」の定義は存在しません．速度場$\boldsymbol{u} = (u, v)$において渦度
$$
\begin{aligned}
\omega = \frac{\partial v}{\partial x} - \frac{\partial u}{\partial y}
\end{aligned}
$$
が定義され，渦度の高い領域が渦としてしばしば考えられますが，渦の存在しない平行せん断流においても渦度の高い領域は存在します．そのため客観的に渦を識別する様々な方法が提案されています．ここではオイラー的，つまり速度場の瞬時場（snapshot）から渦を識別する方法を紹介したいと思います．

オイラー的渦識別法

オイラー的渦識別法は数多提案されていますが，ここでは速度勾配テンソル$\mathbf{D}$に基づいた$Q$-criterion，$\lambda_2$-criterion，$\lambda_\rm{ci}$-criterionの３つを紹介します．

簡単のため2次元の速度場について考えます．速度場 $\boldsymbol{u} = (u ,v)$ の速度勾配テンソル $\mathbf{D}$ は
$$
\begin{aligned}
\mathbf{D} = \nabla \boldsymbol{u} =
\begin{pmatrix}
\frac{\partial u}{\partial x} & \frac{\partial u}{\partial y} \\
\frac{\partial v}{\partial x} & \frac{\partial v}{\partial y} \\
\end{pmatrix}.
\end{aligned}
$$

・$\boldsymbol{Q}$-criterion
Hunt et al. (1995)は渦識別法として$Q$-criterionを提案しました．現在渦の可視化にもっともよく使われている方法だと思います．
変形速度テンソル$\mathbf{S} = \frac12 [\mathbf{D} + \mathbf{D}^T]$，渦度テンソル$\mathbf{\Omega} = \frac12 [\mathbf{D} - \mathbf{D}^T]$としたときに，$Q$-criterionでは
$$
Q = \frac12 [|\mathbf{\Omega}|^2 - |\mathbf{S}|^2] > 0
$$
の領域が渦として同定されます．
2次元の場合についてはHuntらよりも先にOkubo (1970)とWeiss (1991)がそれぞれ独立に導いているためOkubo-Weiss criterionとも呼ばれます（呼ぶべきです）．

・$\boldsymbol{\lambda_2}$-criterion
Jeong and Hussain (2006)によって提案された$\lambda_2$-criterionでは圧力が最小となる領域を渦として識別するようにデザインされています．
$$
\lambda_2(\mathbf{\Omega}^2 + \mathbf{S}^2) < 0
$$
が渦として定義されます．ここで$\lambda_2$は2番めの固有値を意味します．

・$\boldsymbol{\lambda_\rm{ci}}$-criterion (swirling-strength)
Zhou et al. (1999)が提案した$\lambda_\rm{ci}$-criterionでは，速度勾配テンソル$\mathbf{D}$の固有値$\lambda_k$の虚部が渦として識別されます．
$$
\lambda_\rm{ci} = \rm{Img}\{ \lambda_k \} > 0
$$

コード

PIV（粒子画像流速測定方法）によって得られた2次元翼周りの流れについて適用してみます．
位置ベクトルの2次元配列X,Yと速度ベクトルの2次元配列u,vがあるとします．xn $\times$ ynの要素数があるとします．

import numpy as np
from numpy import sqrt, log, pi, sin, cos, abs, cross
import matplotlib.pyplot as plt
from numpy.linalg import norm, eig
import matplotlib.pyplot as plt

# 速度勾配の計算
du_dx = (u[1:-1, 2:] - u[1:-1, 0:-2])/(2*h)
du_dy = (u[0:-2, 1:-1] - u[2:, 1:-1])/(2*h)
dv_dx = (v[1:-1, 2:] - v[1:-1, 0:-2])/(2*h)
dv_dy = (v[0:-2, 1:-1] - v[2:, 1:-1])/(2*h)
du_dx = du_dx.ravel() # .ravel()で2次元配列を1次元に
du_dy = du_dy.ravel()
dv_dx = dv_dx.ravel()
dv_dy = dv_dy.ravel()

Q = []
LAMBDA2 = []
SWIRL = []
for k in np.arange(len(du_dx)):
    D = np.array([du_dx[k], du_dy[k], dv_dx[k], dv_dy[k]])
    D = D.reshape((2, 2))  # 2x2の速度勾配テンソルに
    S = 0.5*(D + D.T)  # 変形速度テンソル 
    OMEGA = 0.5*(D - D.T)  # 渦度テンソル
    q = 0.5*(norm(OMEGA)** 2 - norm(S)** 2)  # Q-criterion
    eigvals2, eigvecs2 = eig(OMEGA**2 + S**2)  # lambda2
    lambda2 = np.min(eigvals2)
    eigvals, eigvecs = eig(D)  # 速度勾配テンソルの固有値問題
    swirl = np.imag(eigvals)  # 速度勾配テンソルの固有値の虚部
    Q.append(q)
    LAMBDA2.append(lambda2)
    SWIRL.append(swirl[0])
Q = np.array(Q)
Q_2D = Q.reshape((xn - 2, yn - 2))  # .reshape(())で2次元配列に戻す
LAMBDA2 = np.array(LAMBDA2)
LAMBDA2_2D = LAMBDA2.reshape((xn - 2, yn - 2))
SWIRL = np.array(SWIRL)
SWIRL_2D = SWIRL.reshape((xn - 2, yn - 2))

プロットして比較してみましょう．

$\lambda_2$-criterionでは圧力が低い部分のみを抽出できています．$Q$-criterionと$\lambda_\rm{ci}$-criterionの結果は似ていますが，$\lambda_\rm{ci}$-criterionの方が後縁渦も抽出できています．そのため私は$\lambda_\rm{ci}$-criterionを好んで使用しています．

次回はラグランジュ的渦識別法について書きたいと思います．