Go to Qiita Advent Calendar 2024 Top
LoginSignup
5
5

Delete article

Deleted articles cannot be recovered.

Draft of this article would be also deleted.

Are you sure you want to delete this article?

@ShotaDeguchi(Shota DEGUCHI)

キャビティ流れの有限差分法シミュレーション

Last updated at Posted at 2023-11-21

初めに

有限差分法の基礎的な勉強の後，もう少しだけ応用（らしい）問題を取り上げようと思った．

有名な2次元キャビティ流れのシミュレーションを，有限差分法によって実行することにした．

間違いを含む可能性があるので，鵜呑みにしないよう慎重に読んでいただき，その場合はご指摘いただきたい．

コードはここにある．

Methodology

支配方程式

水などの非圧縮性のニュートン流体を考える．支配方程式は，質量保存則と運動量保存則を記す以下の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}

ここで，$\mathbf{u}$は速度，$p$は圧力，$\mathbf{f}$は外力である（空間次元は$M$次元とする i.e. $\mathbf{x} \in \Omega \subset \mathbb{R}^M$）．ただし，キャビティ流れは外力が作用しない問題設定なので，$\mathbf{f} = \mathbf{0}$とする．また，$\mathrm{Re}$はReynolds数と呼ばれる無次元量であり，慣性と粘性の比を示す．$\mathrm{Re}$が大きいとき，その流れ場では慣性が卓越しており，流体が激しく動く．$\mathrm{Re}$が小さいとき，粘性が支配的となっており，流体粒子たちが規則正しく並んで動く．

空間方向の離散化

移流，拡散，圧力勾配を何らかの差分スキームで近似する．

ここでは，移流を1次精度風上差分，QUICKEST，KKの3つのスキームで近似し（各手法については，以前の投稿を参照），拡散，圧力勾配については2次精度中心差分で以て表現する．

例えば，以下のように近似する．

\begin{align}
&\mathrm{advection:} \; c \frac{\partial \phi}{\partial x} 
    \approx c \frac{\phi_{i+1} - \phi_{i-1}}{2 \Delta x} 
    - \frac{|c| \Delta x}{2} \frac{\phi_{i+1} - 2 \phi_{i} + \phi_{i-1}}{(\Delta x)^2} \\
&\mathrm{diffusion:} \;\; \frac{\partial^2 \phi}{\partial x^2} 
    \approx \frac{\phi_{i+1} - 2 \phi_{i} + \phi_{i-1}}{(\Delta x)^2} \\
&\mathrm{gradient:} \;\;\;\;\ \frac{\partial \phi}{\partial x} 
    \approx \frac{\phi_{i+1} - \phi_{i-1}}{2 \Delta x}
\end{align}

ここで，$\phi_{i}$は節点$i$に与えられた物理量，$\Delta x$は空間解像度である．

時間方向の離散化

Chorinの射影法（Chorin1968）を用いる．これは，分離型解法，演算子分割法などと呼ばれることもある1．

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

ここで，$\mathbf{u}^{(*)}$は新たに導入された中間速度である（仮の速度などとも呼ばれる）．

中間速度$\mathbf{u}^{(*)}$が連続の式を満たしていなくても目を瞑ることにし，$\mathbf{u}^{(n+1)}$でこれが満たされると考える．すなわち，$\nabla \cdot \mathbf{u}^{(n+1)} = 0$を要求し，以下の圧力Poisson方程式（PPE: Pressure Poisson Equation）を得る．

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

これは，$\nabla \cdot \mathbf{u}^{(n+1)} = 0$を満たすような圧力場を定めておき，そのような速度場に向けて中間速度を圧力勾配によって修正する，という操作だと解釈できる．

纏めると，時間積分を以下の流れで実行する：

# Description Equation
1 中間の速度場を得る． $\mathbf{u}^{(*)} = \mathbf{u}^{(n)} + \Delta t \left( - \left( \mathbf{u}^{(n)} \cdot \nabla \right) \mathbf{u}^{(n)} + \nabla^2 \mathbf{u}^{(n)} / \mathrm{Re} + \mathbf{f}^{(n)} \right)$
2 PPEを解き次の圧力場を得る． $\nabla^2 p^{(n+1)} = \nabla \cdot \mathbf{u}^{(*)} / \Delta t$
3 次の速度場へ修正する． $\mathbf{u}^{(n+1)} = \mathbf{u}^{(*)} + \Delta t \left( - \nabla p^{(n+1)} \right)$

全く違う話題だが，そう云えば，と思い，mermaid記法を使ってみる．が，醜い（こういうのを伝えるのには向いていないか）．

空間解像度

流体内部では絶え間なく渦が発生，消散している．大きい渦は小さい渦を形成し，小さい渦がさらに小さい渦を形成し，さらにさらに，，，といった具合である．この過程でエネルギーは大きい渦から小さい渦へ受け渡され，最終的には極めて小さい（Kolmogorovの長さスケール）渦となり，粘性（摩擦）の影響により熱エネルギーとして散逸する．

渦の最小単位とされるKolmogorovの長さスケール$\eta$は，次元解析から概算できる．単位時間当たりの乱れエネルギー供給率$\delta \ [ \mathrm{m^2/s^3} ]$は，エネルギー$U^2$を時間$U / L$で除すことにより，

\delta \sim \mathcal{O} \left[ \frac{U^{3}}{L} \right]

と得られる．一方，乱れエネルギーの散逸率$\epsilon \ [ \mathrm{m^{2}/s^{3}} ]$は，

\epsilon \sim \mathcal{O} \left[ \frac{\nu^{3}}{\eta^{4}} \right]

程度である．供給と散逸が同程度に発生しているとすると，特徴長さ$L$と渦の最小長さ$\eta$の比が得られる．

\frac{L}{\eta} \sim \mathcal{O} \left[ \mathrm{Re}^{3/4} \right]

すなわち，高Reynolds数の流れ場では，エネルギーを供給するスケールに対して，エネルギーが散逸するスケールが非常に小さいということである．計算の都合上の話をすれば，1軸方向において，特徴長さ$L$と渦の最小長さ$\eta$の比がReynolds数の冪乗に比例するということであり，渦の解像に必要となる格子点の数が分かる．特に，空間$M$次元の広がりを持つ問題では，

\left( \frac{L}{\eta} \right)^{M}
\sim \mathcal{O} \left[ \left( \mathrm{Re}^{3/4} \right)^{M} \right]

となる．例えば，$\mathrm{Re} = 1,000$の流れ場をシミュレートしたいのならば，1軸方向に
$\left( L / \eta \right)^{M} \sim 180^{M}$程度の格子点が要求される．2次元空間では$\left( L / \eta \right)^{2} \sim 32,000$，3次元空間では$\left( L / \eta \right)^{3} \sim 5,600,000$というように，要求される計算資源は指数的に増加する．また，実際にシミュレートしたい現象のReynolds数は数万スケール以上であることが多いので2，それはそれは莫大な計算資源が要求されることになる．

キャビティ流れでは，特徴長さ$L$を$L = 1$として，渦の最小長さ$\eta$が凡そ分かる．すなわち，$\eta$より小さい解像度を用いれば，（凡そ）全ての渦が解像できると考えられる．再度，$\mathrm{Re} = 1,000$とすれば$\eta \approx 5.5 \times 10^{-3}$なので，$\Delta x = 5 \times 10^{-3}$などとすれば，最も小さい渦まで解像できると考えられる（はずである）．

ただし，先に述べたように，全ての渦を解像しようとすると，要求される計算資源は実に膨大となる（DNS: Direct Numerical Simulation）．そのような場合には，何らかの乱流モデルを導入することで，解像し切れない渦による熱散逸を表現する（LES: Large Eddy Simulation）3．

Smagorinskyモデル

ここでの議論は，この動画が勉強になる．

渦の最小長さより大きい計算格子を使ってシミュレートする場合，計算格子より大きい（グリッドスケール，GS）渦はきちんと解像できるが，それより小さい（サブグリッドスケール，SGS）渦は解像し切れないので，熱へのエネルギー散逸が表現できない．ここで，SGSの渦の効果をモデル化して計算するのがLESである．先に述べたように，小さい渦はエネルギーを熱に散逸させるため，粘性の効果としてエネルギー消散を表現することにする．ここで，分子粘性に加えて「渦粘性」を導入する．すなわち，

\nu \nabla^2 \mathbf{u} \rightarrow \left( \nu + \nu_{T} \right) \nabla^2 \mathbf{u}

として，渦によるエネルギー消散を粘性によるエネルギー消散に吸収させる．ここで，$\nu_{T}$は渦粘性である．

Smagorinskyモデル（Smagorinsky1963）は，広く普及している渦粘性モデルの一つである．エネルギーカスケードにおける順カスケードを表現する一方，逆カスケードは考慮していない点に注意する必要があるが，数値安定性などのメリットが大きい．

\begin{split}
    \nu_{T} &= \left( C_{S} \Delta \right)^{2} | S | \\
            &= \left( C_{S} \Delta \right)^{2} \left( 2 S_{ij} S_{ij} \right)^{1/2}
\end{split}

ここで，$C_{S} \left( \in (0, 1) \right)$はSmagorinsky定数，$\Delta$はグリッドサイズ，$S_{ij}$はひずみ速度テンソルである．

Smagorinsky定数$C_{S}$が唯一のモデル定数である4．Lilly（Lilly1966）は，一様かつ等方的な乱流に対して理論値を与えており，$C_{S} = 0.173$と結論付けている（そのため，Smagorinsky-Lillyモデルと呼ばれることもある）．

経験的に，$C_{S} = 0.173$では渦粘性の効果が強すぎることが報告されており，もう少し小さい値にすることが多いらしい．例えば，Deardorff1970は$C_{S} = 0.1$が実験と良く合致すると報告しており，Ansys Fluentも広い問題で実験と合致する値として$C_{S} = 0.1$をデフォルト値に採用している．本当は，事前にパラメトリックスタディを行うべきなのだろうが，今回は諸々の問題に目を瞑り，$C_{S} = 0.1$とする．さらに本当は，壁面付近で渦が形成されないよう（壁面付近で渦粘性が入らないよう）減衰関数（例えば，van Driestのwall damping function）を導入すべきであるが，現状，実装していない．

時間増分

Chorinの射影法は速度を陽的に更新するため，時間増分を慎重に選択する必要がある．

Courant数の条件

Courant+1967は，移流に関する時間増分の条件を与えている．

その主張は，以下で定義されるCourant数$C$は$1$より大きくなってはならない，というものである．これにより，移流速度$c$の移流を安定にシミュレートするための時間増分$\Delta t_{C}$を得る．

\begin{align}
C = \sum_{j}^{M} \frac{| c_{j} | \Delta t_{C}}{\Delta x_{j}} &\le 1 \\
\therefore \Delta t_{C} &\le \left( \sum_{j}^{M} \frac{ |c_{j}| }{\Delta x_{j}} \right)^{-1}
\end{align}

ここで，$j$は空間の軸を示すインデックスであり，$M$は問題の空間次元である．各軸方向に同じ解像度，同じ移流速度が与えられているとすれば，

\Delta t_{C} \le \frac{\Delta x}{|c| M}

が要求される．

拡散数の条件

Neumann&Richtmyer1950は，拡散に関する時間増分の条件を与えている．

その主張は，以下で定義される拡散数$D$は$1 / 2$より大きくなってはならない，というものである．これにより，拡散係数$d$の拡散を安定にシミュレートするための時間増分$\Delta t_{D}$を得る．

\begin{align}
D = \sum_{j}^{M} \frac{d \Delta t_{D}}{(\Delta x_{j})^2} &\le \frac{1}{2} \\
\therefore \Delta t_{D} &\le \frac{1}{2} \left( \sum_{j}^{M} \frac{d}{(\Delta x_{j})^2} \right)^{-1}
\end{align}

ここで，$j$，$M$には先ほどと同じ意味を与えた．また，各軸方向に同じだけの解像度を選択すれば，

\Delta t_{D} \le \frac{1}{2} \frac{(\Delta x)^2}{d M}

を満たせば良いことになる．高粘性流体のシミュレーションを実行する際には，Courant数の条件より拡散数の条件が卓越し，かつ拡散数の条件は空間解像度の2乗の速度で抑えつけてくるので，陰解法が好まれたりする．

時間増分の選択

Courant数の条件，拡散数の条件を満たしても，外力などが作用する場合には，その効果によってさらに時間増分が制限される．全ての条件を事前に検討するのは大変なので，Courant数の条件，拡散数の条件の両者を満たす数値の内，なるべく大きいものを時間増分として選択することが一般的である．

\Delta t = f_{\mathrm{safety}} \min \left( \Delta t_{C}, \Delta t_{D} \right)

ここで，$f_{\mathrm{safety}}$は安全率である．今回は$f_{\mathrm{safety}} = 0.4$とした．

グリッド

差分法では，グリッドの選択（ここでは，格子の大きさではなく，変数の配置の意味）が重要となる（Arakawa&Lamb1977）．

説明は省くが，安定性と実装の容易さのバランスから，本稿ではArakawa-B型グリッドを採用する．最も原始的なArakawa-A型グリッドも選択肢の一つであるが，このグリッドは圧力勾配を中心差分近似したとき，checkerboard patternの圧力振動が発生し易い．

Results

キャビティ流れは，Reynolds数が比較的低い範囲では，定常問題として近似できる．本稿では，以下が満たされたとき，定常状態に到達したと判定することとした．

\left( \frac{\sum_{i, j} | u_{i,j}^{(n+1)} - u_{i,j}^{(n)} |^{2}}{\sum_{i, j} | u_{i,j}^{(n)} |^{2}} \right)^{1/2} \le \varepsilon

ここで，$\varepsilon = 10^{-8}$とした．

領域は，単位矩形領域である（$\Omega = (0, 1)^2$）．初期条件は，全ての変数をゼロで初期化した．境界条件は，速度のDirichlet条件として上壁でのみ$\mathbf{u} = [1, 0]^{\top}$を与え，それ以外の壁で非すべり条件を与えた．圧力の条件として，全ての壁で斉次Neumann条件を与えた．矩形領域の左下の頂点で斉次Dirichlet条件を与えた．

移流項の離散化スキームの比較

格子系の数値計算法では，移流項の取り扱いが重要となる．なぜなら，移流項を差分近似した結果，「数値粘性」と呼ばれる，分子粘性とは異なる粘性 / 拡散の効果が現れるためである．

本節では，

  • A: 1次精度風上差分: 2次精度中心差分 + 2階の数値粘性
  • B: QUICKEST: 4次精度中心差分 + 4階の数値粘性
  • C: KK: 4次精度中心差分 + 4階の数値粘性

の3者の比較を行う．

空間解像度は$\Delta x = 5 \times 10^{-3}$，時間増分は$\Delta t = 1 \times 10^{-3}$とした．十分細かい解像度を用いている（つもり）なので，渦粘性は導入していない．

速度場のアニメーション，および定常状態での各変数を示す．

A B C
$| \mathbf{u} |$ ezgif.com-gif-maker.gif ezgif.com-gif-maker.gif ezgif.com-gif-maker.gif
$| \mathbf{u} |$ velocity_norm.png velocity_norm.png velocity_norm.png
$p$ pressure.png pressure.png pressure.png
$u(y=0.5)$ u.png u.png u.png
$v(x=0.5)$ v.png v.png v.png

中央部スライスでの水平，鉛直速度は，Ghia+1982，Erturk+2005と比較している．

以前の投稿でも述べた通り，1次精度風上差分は数値粘性の影響が大きく，本来より粘り気の強い流体をシミュレートしたような結果になった．QUICKEST，KKの違いは殆ど分からないが，注意深く見れば，KKの方が僅かに流速が大きい．

渦粘性モデルの有無の比較

Smagorinskyの渦粘性モデルを導入し，その有無による効果を比較する．

上記のシミュレーション（$\mathrm{Re} = 1,000$）では$\eta \sim 5.6 \times 10^{-3}$であり，$\Delta x = 5 \times 10^{-3}$の解像度で凡そ全てのスケールの渦が解像できているはずである．同じ問題を，もう少し解像度を下げ，$\Delta x = 1 \times 10^{-2}$ほどの解像度で計算したい，というシナリオを考える．

本節では，

  • a: $\Delta x = 5 \times 10^{-3}$，Smagorinskyの渦粘性: 無し
  • b: $\Delta x = 1 \times 10^{-2}$，Smagorinskyの渦粘性: 無し
  • c: $\Delta x = 1 \times 10^{-2}$，Smagorinskyの渦粘性: 有り

の3者の比較を行う．全ての計算でKKスキームを用い，$\Delta t = 1 \times 10^{-3}$とする．

a b c
$| \mathbf{u} |$ velocity_norm.png velocity_norm.png velocity_norm.png
$p$ pressure.png pressure.png pressure.png
$u(y=0.5)$ u.png u.png u.png
$v(x=0.5)$ v.png v.png v.png

渦粘性を導入せずに解像度が異なる計算結果（aとb）を比較すると，解像度が低い場合の方が，参照解（Ghia+1982，Erturk+2005）より少し大きい速度を得る結果となった．

同じ解像度で渦粘性の有無が異なる計算結果（bとc）を比較すると，目に見える差異は殆ど無い．違いと言えば，定常状態に入るまでの時間（無次元化時間）が1（秒？）異なるほどであった．本当は壁でのdampingを入れないといけないので，その辺りに問題がある可能性はある．

終わりに

幾つかの気になっていた事項を文字にして纏めた．

Appendix A: 無次元化

非圧縮性NS式の無次元化をメモしておく．

元の（次元が有る）方程式は，

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

である．

特徴長さ$L$，特徴速度$U$を用いて，有次元の物理量$\phi$と無次元の物理量$\bar{\phi}$を関係づける．

\mathbf{x} = L \bar{\mathbf{x}}, \
\mathbf{u} = U \bar{\mathbf{u}}, \
t = \frac{L}{U} \bar{t}

質量保存則

\begin{align}
    \nabla \cdot \mathbf{u}
    &= \frac{\partial u_j}{\partial x_j} \\
    &= \frac{\partial U \bar{u}_j}{\partial L \bar{x}_j} \\
    &= \frac{U}{L} \frac{\partial \bar{u}_j}{\partial \bar{x}_j} \\
    &= 0
\end{align}

$U/L \neq 0$ならば，

\frac{\partial \bar{u}_j}{\partial \bar{x}_j} = 0

として，無次元化された質量保存則を得る．

運動量保存則

\begin{align}
    \frac{\partial u_j}{\partial t}
    + u_i \frac{\partial}{\partial x_i} u_j
    &= - \frac{1}{\rho} \frac{\partial}{\partial x_j} p
    + \nu \frac{\partial}{\partial x_i} \frac{\partial}{\partial x_i} u_j
    + f_j \\
    (\mathrm{LHS})
    &= \frac{\partial U \bar{u}_j}{\partial L/U \bar{t}}
    + U \bar{u}_i \frac{\partial}{\partial L \bar{x}_i} U \bar{u}_j \\
    &= \frac{U^2}{L} \left( \frac{\partial \bar{u}_j}{\partial \bar{t}}
    + \bar{u}_i \frac{\partial}{\partial \bar{x}_i} \bar{u}_j \right) \\
    (\mathrm{RHS})
    &= - \frac{1}{\rho} \frac{\partial}{\partial L \bar{x}_j} p
    + \nu \frac{\partial}{\partial L \bar{x}_i} \frac{\partial}{\partial L \bar{x}_i} U \bar{u}_j
    + f_j \\
    &= - \frac{1}{\rho L} \frac{\partial}{\partial \bar{x}_j} p
    + \nu \frac{U}{L^2} \frac{\partial}{\partial \bar{x}_i} \frac{\partial}{\partial \bar{x}_i} \bar{u}_j
    + f_j \\
    \therefore
    \frac{U^2}{L} \left( \frac{\partial \bar{u}_j}{\partial \bar{t}}
    + \bar{u}_i \frac{\partial}{\partial \bar{x}_i} \bar{u}_j \right)
    &= - \frac{1}{\rho L} \frac{\partial}{\partial \bar{x}_j} p
    + \nu \frac{U}{L^2} \frac{\partial}{\partial \bar{x}_i} \frac{\partial}{\partial \bar{x}_i} \bar{u}_j
    + f_j \\
\end{align}

$U^2/L \neq 0$ならば，

\begin{align}
    \frac{\partial \bar{u}_j}{\partial \bar{t}}
    + \bar{u}_i \frac{\partial}{\partial \bar{x}_i} \bar{u}_j
    &= - \frac{1}{\rho U^2} \frac{\partial}{\partial \bar{x}_j} p
    + \frac{\nu}{UL} \frac{\partial}{\partial \bar{x}_i} \frac{\partial}{\partial \bar{x}_i} \bar{u}_j
    + \frac{L}{U^2} f_j \\
    &= - \frac{\partial}{\partial \bar{x}_j} \bar{p}
    + \frac{1}{\mathrm{Re}} \frac{\partial}{\partial \bar{x}_i} \frac{\partial}{\partial \bar{x}_i} \bar{u}_j
    + \frac{L}{U^2} f_j \\
\end{align}

ここで，$\bar{p} = p / \rho U^2$を導入した．また，$f_j = g_j + h_j$として，重力加速度$g_j$を取り出す．

\begin{align}
    \frac{L}{U^2} f_j
    &= \frac{L}{U^2} (g_j + h_j) \\
    &= \frac{1}{\mathrm{Fr}_j^2} + \frac{L}{U^2} h_j \\
\end{align}

より，

\begin{align}
    \frac{\partial \bar{u}_j}{\partial \bar{t}}
    + \bar{u}_i \frac{\partial}{\partial \bar{x}_i} \bar{u}_j
    &= - \frac{\partial}{\partial \bar{x}_j} \bar{p}
    + \frac{1}{\mathrm{Re}} \frac{\partial}{\partial \bar{x}_i} \frac{\partial}{\partial \bar{x}_i} \bar{u}_j
    + \frac{1}{\mathrm{Fr}_j^2} + \frac{L}{U^2} h_j \\
\end{align}

ここで，無次元量

  • Reynolds数（$\mathrm{Re} = UL / \nu$）
  • Froude数（$\mathrm{Fr}_j = U / \sqrt{L g_j}$）

を導入した．

Appendix B: Reynoldsの相似則

学部3年次の水理学の講義では，Reynolds数が等しい2つの流れ場は，本質的に等価であると学んだ．当時は「そうなのか，覚えておこう」と疑いもせず鵜呑みにしていたが，実際に確認したことがないことを思い出した．学部で学んだことは本当なのだろうか．

3つのパラメータ（代表流速$U \ [\mathrm{m/s}]$，代表長さ$L \ [\mathrm{m}]$，動粘性係数$\nu \ [\mathrm{m^2/s}]$）を動かして，$\mathrm{Re} = 1,000$の流れ場を作る．

$\Delta x / L = 5 \times 10^{-3}$とし，時間増分はCourant数の条件，拡散数の条件の両者を満たすものに安全率$f_{\mathrm{safety}} = 0.4$を乗じたものを採用した．

1 2 3
$U$ $1,000$ $1$ $1$
$L$ $1$ $1,000$ $1$
$\nu$ $1$ $1$ $0.001$
$| \mathbf{u} |$ vel_035000.png vel_035000.png vel_035000.png

カラーバーのレンジや1辺の長さは違うが，分布として同じ流れ場に収束した（定常状態に収束するまでに要した時間積分の回数は同一であった）．過去に学んだことは本当だったようだ．

Appendix C: 移流項の非線形性

(2024/09追記)
こちらに，大変興味深い議論を見つけた．

いま，流速$\mathbf{u}$をFourier級数和で表すことができるとし，

\begin{align}
    \mathbf{u}_k
    &= U (t) 
    \left[
        \begin{array}{c}
            \sin (k_1 x_1) \\
            \sin (k_2 x_2) \\
            \sin (k_3 x_3)
        \end{array}
    \right]
\end{align}

だけ取り出して考える（例えば，von Neumannの線形安定性解析などにも登場する，よくあるやり方）．すると，移流項$(\mathbf{u}_k \cdot \nabla) \mathbf{u}_k$は，

\begin{align}
    \left[ (\mathbf{u}_k \cdot \nabla) \mathbf{u}_k \right]_1
    &= \frac{U^2 k_1}{2} \sin (2 k_1 x_1)
\end{align}

と書け，移流項が，波数$k_1$の波から，波数$2 k_1$の波を生み出している．

Appendix D: 高Reynolds数の流れ

(2024/10追記)
時間が経ったから，別の文章に纏めようかとも思ったが，手法は同じなので，ここに追記することにした．

移流スキームの比較を，もう少し高Reynolds数の流れでやってみようと思った．初めにこの記事を書いたときは，QUICKESTとKKを比較したが，QUICKESTには時間の項も入っているので，最近は，空間の作要素だけ比較すべきではないか，と思っている．ここでは，高Reynolds数の流れ場を対象に，QUICKとKKを比較する．$v (x=0.5)$だけ見る．

$\mathrm{Re} = 5,000$ $\mathrm{Re} = 10,000$
$\mathrm{QUICK}$ Re5e3_QUICK_v.gif Re1e4_QUICK_v.gif
$\mathrm{KK}$ Re5e3_KK_v.gif Re1e4_KK_v.gif

何れも，同様の結果を得た．この問題では，移流スキームがある程度高精度であれば，それらの差は殆ど見え辛いものであろう．まあしかし，なんだ，正しい計算ができたという自信はないし，いま手元にある計算資源と計算スキームではこの結果をえた．という程度である．

(2024/12追記)
John2016に依れば，キャビティ流れには，$\mathrm{Re} = 8,000$前後までは定常解が存在し，高Reynolds数帯で定常解を得るには，何らかの安定化が必要であるとされている（Cazemier+1998, Tiesinga+2002, Bruneau&Saad2006などに基づいている（bifurcation!））．

以下は渦粘性を入れていない結果（離散化パラメータ: $h = 5 \times 10^{-3}, \ \tau = 1 \times 10^{-3}$）．

$\mathrm{Re}$ $t$ $u$ $v$ $\psi$ $p$ $\mathrm{residual}$
$ 1,000$ $92.0$ u.png v.png streamline.png pressure_bar.png residual.png
$ 5,000$ $455.2$ u.png v.png streamline.png pressure_bar.png residual.png
$ 7,500$ $600.0$ u.png v.png streamline.png pressure_bar.png residual.png
$ 8,500$ $600.0$ u.png v.png streamline.png pressure_bar.png residual.png
$10,000$ $600.0$ u.png v.png streamline.png pressure_bar.png residual.png

$t$は収束までに要した無次元化時間（最大で$t=600$まで）．圧力は定数を動かす自由度のもとで一意に定まるので，平均値がゼロになるよう平行移動させている．$|| \cdot ||$は2のノルム．$\mathrm{Re} = 8,500$も十分な時間が経つと収束しそうだ．$\mathrm{Re} = 10,000$は驚くほど収束しない．見た目の変化は案外小さく見えるのだけれど，，

以下はSmagorinsky-Lillyモデルを入れた結果（モデル定数: $C_S = 0.1$）．

$\mathrm{Re}$ $t$ $u$ $v$ $\psi$ $p$ $\mathrm{residual}$
$ 1,000$ $91.9$ u.png v.png streamline.png pressure_bar.png residual.png
$ 5,000$ $455.3$ u.png v.png streamline.png pressure_bar.png residual.png
$ 7,500$ $600.0$ u.png v.png streamline.png pressure_bar.png residual.png
$ 8,500$ $600.0$ u.png v.png streamline.png pressure_bar.png residual.png
$10,000$ $600.0$ u.png v.png streamline.png pressure_bar.png residual.png

案外，収束しない．

無次元化時間$t=1,200$までやってみた．

$C_s$ $\text{point-wise residual}$ $\text{mean residual}$
$0.0$ vel_res_vis.png residual.png
$0.1$ vel_res_vis.png residual.png
$0.2$ vel_res_vis.png residual.png

point-wise residualは，1ステップ間の速度の変化

    \left( \frac{| u^{(n+1)} - u^{(n)} |^2}{| u^{(n)} |^2} \right)^{1/2}
    + \left( \frac{| v^{(n+1)} - v^{(n)} |^2}{| v^{(n)} |^2} \right)^{1/2}

を節点ごとに計算したもの．時間平均などしておらず，$t=1,200$でのスナップショットだが，中央部が暗く，角付近が明るい色になる傾向は，ずっと継続していたことを強調しておく．この図から，角付近のsecondary vortexが激しく変動していて，この辺りが全体の残差の低減を妨げているようだ．全体の残差は，$\mathcal{O} (10^{-6})$程度の振動をしながら$\mathcal{O} (10^{-4})$で停滞していた．

より強い安定化を導入しなければ定常解に収束しないような気がするが，果たしてそれは，本来のReynolds数の流れをシミュレートしたことになるのだろうか．

  1. 現状の私の理解を記しておく．非圧縮性流体シミュレーションの解法として導入されたのがChorinの射影法である．これは，Helmholtz分解に基づき，速度場を2つの直交成分（発散フリーの場と回転フリーの場）に分解する．似た方法に，Fractional Step法というものがある．これはより広い枠組みで（恐らく非圧縮性流体に限らず），演算子分割をより多段階にして（例えば4段階）時間積分を実行するものである．ただし，演算子分割を2段階にしたFractional Step法が最も良く普及しており，これを非圧縮性流体シミュレーションへ適用するとChorinの射影法に似通った形になる． ↩

  2. 航空機の主翼まわりでのReynolds数は$10^{6} \sim 10^{8}$ほどらしい． ↩

  3. 名前のメモだけ残しておく．乱流モデルを導入する必要が無いほどに細かい解像度を用いて直接シミュレートする方法をDNS（Direct Numerical Simulation）と呼ぶ．全ての渦を解像することは諦め，解像し切れない渦による熱散逸をモデル化する方法をLES（Large Eddy Simulation）と呼ぶ．この他に，RANS（Reynolds-Averaged Navier-Stokes），URANS（Unsteady RANS）などもある． ↩

  4. Smagorinsky定数$C_{S}$を動的に変化する変数とするDynamic Smagorinskyモデルもある． ↩

5
5
0

Register as a new user and use Qiita more conveniently

  1. You get articles that match your needs
  2. You can efficiently read back useful information
  3. You can use dark theme
What you can do with signing up
Sign upLogin
5
5

Delete article

Deleted articles cannot be recovered.

Draft of this article would be also deleted.

Are you sure you want to delete this article?