Euler型粒体モデル — 理論と支配方程式

カテゴリ: 流体解析 | 2026-01-15
eulerian-granular-theory
理論と物理の世界へ

概要

🧑‍🎓

先生、Euler型粒体モデルって何ですか?粒子を連続体として扱うんですか?


🎓

その通り。Euler型粒体モデル(Eulerian Granular Model)は、粉体や顆粒の集合体を「粒体相」という擬似連続体として扱い、気体相とともにEuler-Euler法の枠組みで解く手法だ。流動層、空気輸送、サイクロン、粉体混合などの固気二相流で広く使われている。


🧑‍🎓

DEM-CFDとは何が違うんですか?


🎓

DEM-CFDは個々の粒子を離散的に追跡するが、粒子数が数百万〜数十億になる工業スケールでは計算コストが膨大になる。Euler型粒体モデルは粒子群を連続体として扱うため、大規模な系でも現実的な計算時間で解ける。ただし個々の粒子の接触力や形状効果は直接扱えない。


支配方程式

🧑‍🎓

粒体相の方程式はどうなりますか?


🎓

粒体相の運動を記述するために、KTGF(Kinetic Theory of Granular Flow:粒体の運動論)を使う。固相の連続の式は通常のEuler-Eulerと同じだが、固相の応力テンソルにKTGFから導かれる構成則を用いるのが特徴だ。


$$ \frac{\partial (\alpha_s \rho_s)}{\partial t} + \nabla \cdot (\alpha_s \rho_s \mathbf{u}_s) = 0 $$

🧑‍🎓

KTGFって気体分子運動論の粒子版ということですか?


🎓

まさにその通り。粒子の速度揺動をGranular Temperature $\Theta_s$ で特徴づける。


$$ \Theta_s = \frac{1}{3} \langle \mathbf{u}_s' \cdot \mathbf{u}_s' \rangle $$

🎓

Granular Temperatureの輸送方程式は次の通りだ。


$$ \frac{3}{2} \left[ \frac{\partial (\alpha_s \rho_s \Theta_s)}{\partial t} + \nabla \cdot (\alpha_s \rho_s \mathbf{u}_s \Theta_s) \right] = (-p_s \mathbf{I} + \boldsymbol{\tau}_s) : \nabla \mathbf{u}_s + \nabla \cdot (\kappa_s \nabla \Theta_s) - \gamma_s + \phi_{gs} $$

🎓

右辺は順に、せん断による生成、拡散($\kappa_s$は拡散係数)、非弾性衝突による散逸($\gamma_s$)、気体-固体間のエネルギー交換($\phi_{gs}$)だ。


固相の構成則

🧑‍🎓

固相圧力とか粘性はどう決まるんですか?


🎓

固相圧力 $p_s$ はGranular Temperatureと体積分率から求まる。Lun et al.(1984)のモデルが代表的だ。


$$ p_s = \alpha_s \rho_s \Theta_s + 2 \rho_s (1 + e_{ss}) \alpha_s^2 g_0 \Theta_s $$

🎓

ここで $e_{ss}$ は粒子間の反発係数、$g_0$ は動径分布関数だ。$g_0$ は体積分率が最密充填率 $\alpha_{s,max}$(約0.63)に近づくと急増し、粒子が密集した状態での接触圧力を表現する。


構成則モデル例物理量
固相圧力Lun, Syamlal-O'Brien$p_s(\alpha_s, \Theta_s)$
固相粘性Gidaspow, Syamlal$\mu_s(\alpha_s, \Theta_s)$
固相体積粘性Lun et al.$\lambda_s$
摩擦応力Schaeffer, Johnson-Jackson密充填領域の応力
Coffee Break よもやま話

レイノルズの実験(1883年)——乱流発見の瞬間

オズボーン・レイノルズは、管内の水にインクを流す実験で「層流から乱流への遷移」を発見しました。流速を上げていくと、インクの線がある瞬間にグチャグチャに乱れる。この劇的な瞬間を、レイノルズは数学的に $Re = \rho uD/\mu$ という無次元数で表現した。100年以上経った今も、CFDエンジニアが最初に確認するのはこのレイノルズ数です。

各項の物理的意味
  • 時間項 $\partial(\rho\phi)/\partial t$:蛇口をひねった瞬間を思い浮かべてください。最初は水がバタバタと不安定に出て、しばらくすると安定した流れになりますよね? この「変化している最中」を記述するのが時間項です。心臓の拍動で血流が脈打つのも、エンジンのバルブが開閉するたびに流れが変動するのも、すべて非定常現象。では定常解析とは? 「十分時間が経って流れが落ち着いた後」だけを見る——つまりこの項をゼロにする。計算コストが大幅に下がるため、まず定常で解いてみるのがCFDの基本戦略です。
  • 対流項 $\nabla \cdot (\rho \mathbf{u} \phi)$:川に落ち葉を落としたらどうなりますか? 流れに乗って下流に運ばれますよね。これが「対流」——流体の動きが物を運ぶ効果です。暖房の温風が部屋の端まで届くのも、空気という「運び屋」が熱を対流で輸送しているから。ここが面白いところ——この項は「速度×速度」を含むため非線形です。つまり、流れが速くなるとこの項が急激に強くなり、制御が難しくなる。これが乱流の根本原因です。よくある勘違い:「対流と伝導は同じようなもの」→ 全然違います! 対流は流れが運ぶ、伝導は分子が伝える。桁違いの効率差があります。
  • 拡散項 $\nabla \cdot (\Gamma \nabla \phi)$:コーヒーにミルクを入れて放置したことはありますか? かき混ぜなくても、しばらく経つと自然に混ざりますよね。あれが分子拡散です。では次の質問——ハチミツとお水、どちらが流しやすいですか? 当然お水ですよね。ハチミツは粘性($\mu$)が高いから流れにくい。粘性が大きいと拡散項が強くなり、流体は「もったりした」動きになります。レイノルズ数が小さい流れ(ゆっくり、ドロドロ)では拡散が支配的。逆にRe数が大きい流れでは対流が圧倒し、拡散は脇役になります。
  • 圧力項 $-\nabla p$:注射器のピストンを押すと、液体が針先から勢いよく出ますよね? なぜでしょう? ピストン側が高圧、針先が低圧——この圧力差が流体を押す力になるからです。ダムの放水も同じ原理。天気図で等圧線がギュッと密になっている場所では? そう、強風が吹きます。「圧力差があるところに流れが生まれる」——これがナビエ-ストークス方程式の圧力項の物理的意味。ここでの勘違いポイント:CFDの「圧力」は絶対圧ではなくゲージ圧のことが多い。圧縮性解析に切り替えたとたんに結果がおかしくなる場合、絶対圧/ゲージ圧の混同が原因かもしれません。
  • ソース項 $S_\phi$:暖められた空気が上に昇る——なぜでしょう? 周囲より軽く(密度が低く)なったから、浮力で押し上げられるのです。この浮力はソース項として方程式に追加されます。他にも、ガスコンロの炎で化学反応熱が発生する、工場の電磁ポンプで金属溶湯にローレンツ力がかかる…これらはすべて「外部から流体にエネルギーや力を注入する」作用であり、ソース項で表現します。ソース項を忘れるとどうなるか? 自然対流の解析で浮力を入れ忘れると、流体は一切動かない——冬の部屋で暖房をつけたのに暖かい空気が上に行かない、という物理的にありえない結果になります。
仮定条件と適用限界
  • 連続体仮定:クヌッセン数 Kn < 0.01(分子平均自由行程 ≪ 代表長さ)で成立
  • ニュートン流体仮定:せん断応力と歪み速度が線形関係(非ニュートン流体では粘度モデルが必要)
  • 非圧縮性仮定(Ma < 0.3の場合):密度を一定として扱う。マッハ数0.3以上では圧縮性効果を考慮
  • ブシネスク近似(自然対流):密度変化を浮力項のみで考慮し、他の項では一定密度を使用
  • 適用外ケース:希薄気体(Kn > 0.1)、超音速・極超音速流れ(衝撃波捕捉が必要)、自由表面流れ(VOF/Level Set等が必要)
次元解析と単位系
変数SI単位注意点・換算メモ
速度 $u$m/s入口条件で体積流量から換算する際、断面積の単位に注意
圧力 $p$Paゲージ圧と絶対圧の区別。圧縮性解析では絶対圧を使用
密度 $\rho$kg/m³空気: 約1.225 kg/m³@20°C、水: 約998 kg/m³@20°C
粘性係数 $\mu$Pa·s動粘性係数 $\nu = \mu/\rho$ [m²/s] との混同に注意
レイノルズ数 $Re$無次元$Re = \rho u L / \mu$。層流/乱流遷移の判定指標
CFL数無次元$CFL = u \Delta t / \Delta x$。時間刻みの安定性に直結

数値例:円管内層流(d=10mm, L=1m, 水μ=0.001Pa·s, Q=0.1L/min)

Re = ρuD/μ = 998×0.021×0.01/0.001 ≈ 212層流) 圧力損失 ΔP = 128μLQ/(πd⁴) ≈ 68.2 Pa

乱流モデル別の精度比較(後向きステップ、再付着長さ):

k-ε標準5.8h(実験6.1h)
-4.9%
k-ω SST6.0h
-1.6%
RSM6.05h
-0.8%
LES6.12h
+0.3%
実験値6.1h
基準

k-ω SSTは精度とコストのバランスが良く、多くの実務で最初の選択肢になります。

簡易計算ツール:流体力学基礎

レイノルズ数 Re = ρuL/μ を計算し、層流/乱流の判定を行います。

CFL数 = u·Δt/Δx を計算し、時間刻みの安定性を確認します。

円管内の層流ハーゲン-ポアズイユ流れの圧力損失 ΔP = 128μLQ/(πd⁴)

CFDメッシュの品質管理や乱流モデルの選定に悩む時間を、もっと創造的な設計作業に使えたら。 — Project NovaSolverはそんな実務者の声から生まれました。

Project NovaSolver — CAE実務の課題に向き合う研究開発

「Euler型粒体モデルをもっと効率的に解析できないか?」——私たちは実務者の声に耳を傾け、既存ワークフローの改善を目指す次世代CAEプロジェクトに取り組んでいます。具体的な機能はまだ公開前ですが、開発の進捗をお届けします。

進捗通知を受け取る →