0
0

Delete article

Deleted articles cannot be recovered.

Draft of this article would be also deleted.

Are you sure you want to delete this article?

@NovaSolver

流動層シミュレーターで解き明かす「沸騰する砂」の流体力学

0
Posted at

みなさん、こんな経験ありませんか？コーヒーの粉にスプーンで空気を送り込むと、粉がふわっと舞い上がる。あれ、実は超本格的な工学現象なんです。

「流動層」 ─ 粒子が流体によって浮遊・混合される状態。化学プラントから製薬、食品乾燥まで、あらゆる産業で使われているこの現象、実はたった1つの数式で核心が説明できてしまいます。

今日は、ブラウザで動く流動層シミュレーターを題材に、「粒子が浮き始める瞬間」を数式とコードで完全に理解していきましょう。

ざっくり本質：粒子は「エレベーター」と「エスカレーター」の間を行き来する

流動層のイメージ、こんなふうに捉えてみてください。

  • 固定床（流速が小さい） → 粒子は「階段」。エレベーターは動かない。
  • 最小流動化速度（Umf） → エレベーターがちょうど浮き始める瞬間。
  • 流動層（Umf以上） → 粒子が浮遊しながらぐるぐる。エレベーターが動き出した状態。
  • 終末速度（Ut） → エレベーターの天井を突き破って粒子が飛び出す限界。

「流動化とは、流体の抗力が粒子の見かけ重量を上回った瞬間に始まる」

この一言に尽きます。

数式で理解する：Ergun式と最小流動化速度

流動層の心臓部は、Ergun式という1949年に提案された実験式です。固定床における圧力損失 $\Delta P$ を、流速 $U$ の関数として表します。

Ergun式（固定床の圧力損失）

\frac{\Delta P}{L} = 150 \frac{(1-\varepsilon)^2}{\varepsilon^3} \frac{\mu U}{(\phi_s d_p)^2} + 1.75 \frac{1-\varepsilon}{\varepsilon^3} \frac{\rho_f U^2}{\phi_s d_p}

各項の意味：

  • 第1項（粘性項）：流体の粘性による摩擦損失。流速に比例。
  • 第2項（慣性項）：流体の運動エネルギーによる損失。流速の2乗に比例。
  • $\varepsilon$：空隙率（粒子間の隙間の割合。0〜1）
  • $\phi_s$：球形度（1が完全球。現実は0.6〜0.9）
  • $d_p$：粒子径 [m]
  • $\mu$：流体の粘度 [Pa·s]
  • $\rho_f$：流体の密度 [kg/m³]

最小流動化速度 Umf の導出

流動化が始まる条件は、「圧力損失による上向きの力」が「粒子層の見かけ重量」と釣り合うことです：

\Delta P \cdot A = (\rho_p - \rho_f) g (1-\varepsilon) A L

ここで $\rho_p$ は粒子密度。この式とErgun式を連立させ、Umf について解くと：

\frac{1.75}{\phi_s \varepsilon_{mf}^3} \left( \frac{d_p U_{mf} \rho_f}{\mu} \right)^2 + \frac{150 (1-\varepsilon_{mf})}{\phi_s^2 \varepsilon_{mf}^3} \left( \frac{d_p U_{mf} \rho_f}{\mu} \right) = \frac{d_p^3 \rho_f (\rho_p - \rho_f) g}{\mu^2}

右辺こそが アルキメデス数 (Ar) です：

Ar = \frac{d_p^3 \rho_f (\rho_p - \rho_f) g}{\mu^2}

アルキメデス数は、「浮力」と「粘性力」の比。この値が 小さい（Ar < 10⁴）** と粘性支配でUmfは粒子径の2乗に比例。**大きい（Ar > 10⁵） と慣性支配でUmfは粒子径の平方根に比例します。

実用的に、よく使われる簡略式（Wen & Yu の相関式）では：

Re_{mf} = \frac{d_p U_{mf} \rho_f}{\mu} = \sqrt{27.2^2 + 0.0408 Ar} - 27.2

これでUmfが一発計算できます。

コードで実装する：JavaScriptで流動層計算機を作る

では、実際に動くコードを書いてみましょう。シミュレーターのコア計算を再現します。

/**
 * 流動層 核心計算モジュール
 * 入力: 粒子径, 粒子密度, 流体密度, 流体粘度, 空隙率, 球形度
 * 出力: 最小流動化速度Umf, 終末速度Ut, アルキメデス数Ar, 圧力損失曲線
 */
class FluidizedBedCalculator {
  constructor(d_p, rho_p, rho_f, mu, epsilon, phi_s = 1.0, g = 9.81) {
    this.d_p = d_p;           // 粒子径 [m]
    this.rho_p = rho_p;       // 粒子密度 [kg/m^3]
    this.rho_f = rho_f;       // 流体密度 [kg/m^3]
    this.mu = mu;             // 流体粘度 [Pa·s]
    this.epsilon = epsilon;   // 空隙率 [-]
    this.phi_s = phi_s;       // 球形度 [-]
    this.g = g;               // 重力加速度 [m/s^2]
  }

  // アルキメデス数
  archimedesNumber() {
    const { d_p, rho_f, rho_p, mu, g } = this;
    return (d_p ** 3 * rho_f * (rho_p - rho_f) * g) / (mu ** 2);
  }

  // 最小流動化速度 Umf [m/s] (Wen & Yu 相関式)
  minimumFluidizationVelocity() {
    const Ar = this.archimedesNumber();
    const Re_mf = Math.sqrt(27.2 ** 2 + 0.0408 * Ar) - 27.2;
    return (Re_mf * this.mu) / (this.d_p * this.rho_f);
  }

  // 終末速度 Ut [m/s] (球体の自由沈降、Allenの式)
  terminalVelocity() {
    const { d_p, rho_p, rho_f, mu, g } = this;
    const Ar = this.archimedesNumber();
    let Re_t;
    if (Ar < 3.6) { // ストークス域
      Re_t = Ar / 18;
    } else if (Ar < 10**5) { // 中間域 (Allen)
      Re_t = (Ar / 18) ** 0.8; // 近似式
    } else { // ニュートン域
      Re_t = Math.sqrt(Ar / 0.33);
    }
    return (Re_t * mu) / (d_p * rho_f);
  }

  // Ergun式による圧力損失 [Pa/m]
  pressureDrop(U) {
    const { d_p, rho_f, mu, epsilon, phi_s } = this;
    const term1 = 150 * ((1 - epsilon) ** 2) / (epsilon ** 3) * (mu * U) / ((phi_s * d_p) ** 2);
    const term2 = 1.75 * (1 - epsilon) / (epsilon ** 3) * (rho_f * U ** 2) / (phi_s * d_p);
    return term1 + term2;
  }

  // 圧力損失曲線データ生成（0〜Utまで50点）
  generatePressureCurve(points = 50) {
    const Ut = this.terminalVelocity();
    const Umf = this.minimumFluidizationVelocity();
    const data = [];
    for (let i = 0; i <= points; i++) {
      const U = (Ut / points) * i;
      const dP = this.pressureDrop(U);
      data.push({ U, dP });
    }
    return { data, Umf, Ut, Ar: this.archimedesNumber() };
  }
}

// 使用例：シリカサンド（d_p=200μm, ρ_p=2650, 空気中）
const calc = new FluidizedBedCalculator(
  200e-6,  // 0.2mm
  2650,    // シリカ密度
  1.2,     // 空気密度
  1.8e-5,  // 空気粘度
  0.4,     // 空隙率
  0.8      // 球形度
);

const result = calc.generatePressureCurve();
console.log(`Ar = ${result.Ar.toExponential(3)}`);
console.log(`Umf = ${(result.Umf * 1000).toFixed(2)} mm/s`);
console.log(`Ut = ${(result.Ut * 1000).toFixed(2)} mm/s`);

このコード、実際にNode.jsやブラウザのコンソールで動かせます。たった50行足らずで、流動層の核心計算が再現できるんです。

数値例で確かめる：シリカサンドを空気で流動化

実際に数値を入れて、計算過程を追ってみましょう。

条件：

  • 粒子：シリカサンド（$d_p = 200,\mu m$, $\rho_p = 2650,kg/m^3$）
  • 流体：空気（$\rho_f = 1.2,kg/m^3$, $\mu = 1.8\times 10^{-5},Pa\cdot s$）
  • 空隙率 $\varepsilon = 0.4$, 球形度 $\phi_s = 0.8$

ステップ1：アルキメデス数を計算

Ar = \frac{(200\times10^{-6})^3 \times 1.2 \times (2650 - 1.2) \times 9.81}{(1.8\times10^{-5})^2}

分子：$8.0\times10^{-15} \times 1.2 \times 2648.8 \times 9.81 = 2.49\times10^{-10}$
分母：$3.24\times10^{-10}$

Ar \approx 0.769

ステップ2：レイノルズ数（Umf時）を計算

Re_{mf} = \sqrt{27.2^2 + 0.0408 \times 0.769} - 27.2
= \sqrt{739.84 + 0.0314} - 27.2
= 27.2006 - 27.2 \approx 0.0006

ステップ3：Umfを求める

U_{mf} = \frac{0.0006 \times 1.8\times10^{-5}}{200\times10^{-6} \times 1.2}
= \frac{1.08\times10^{-8}}{2.4\times10^{-4}} = 4.5\times10^{-5}\,m/s

結果：Umf = 0.045 mm/s（非常に小さい！）

これは粒径200μmと小さいため、わずかな流速で流動化が始まることを示しています。実際の装置では、この数十倍の速度で運転します。

シミュレーターで遊ぶ：3つの実験

流動層シミュレーターを開いて、実際に試してみましょう。

実験1：粒子径を変えてみる

  • 初期値：$d_p = 200,\mu m$, $\rho_p = 2500,kg/m^3$（ガラスビーズ）
  • 操作：粒子径スライダーを100μm → 500μm → 1000μmと動かす

結果：

粒子径 Umf アルキメデス数 気づき
100μm 0.01 mm/s 0.096 ほとんど浮く
500μm 0.28 mm/s 12.0 10倍以上に
1000μm 1.12 mm/s 96.1 さらに4倍

なぜ？ Umfは粒子径の約2乗に比例します。径が10倍になるとUmfは約100倍に。粒子が重くなれば、浮かせるのに強い風が必要なのは直感的ですよね。

実験2：粒子密度を変えてみる

  • 固定：$d_p = 300,\mu m$
  • 操作：密度を $\rho_p = 1000$（水と同等）→ 2500（ガラス）→ 7800（鉄）

結果：

  • 密度1000：Umf = 0.02 mm/s（ほとんど浮く）
  • 密度2500：Umf = 0.10 mm/s
  • 密度7800：Umf = 0.31 mm/s

ポイント：密度差 $(\rho_p - \rho_f)$ がUmfに効きます。軽い粒子ほど小さな流速で流動化。これ、浮力の概念そのものです。

実験3：空隙率を変えて圧力損失曲線を見る

  • 固定：$d_p = 200,\mu m$, $\rho_p = 2500$
  • 操作：空隙率 $\varepsilon$ を 0.35 → 0.45 → 0.55

結果：

  • $\varepsilon = 0.35$：圧力損失曲線が急峻。固定床領域の勾配が大きい。
  • $\varepsilon = 0.45$：中程度。
  • $\varepsilon = 0.55$：曲線が寝る。同じ流速でも圧損が小さい。

なぜ？ Ergun式の分母に $\varepsilon^3$ が入っています。空隙率が小さくなると、粒子間の流路が狭くなり、圧力損失が急増。「粒子が密に詰まっているほど、風を通すのに力がいる」 という直感と一致します。

現場でハマるポイント：理論と現実のギャップ

このシミュレーターは教育用として極めて優秀ですが、実機設計にそのまま使うと痛い目を見ます。

落とし穴1：「Umf = 運転速度」の誤解

シミュレーターで計算されたUmfは「流動が始まる最低速度」。実際の装置では、粒子を活発に混合・熱交換させるため、Umfの2〜10倍で運転します。

例えば、FCC触媒（$d_p \approx 70,\mu m$）のUmfは約0.002 m/sですが、実機の流動床反応器では0.1〜0.3 m/sで運転。「最低速度」と「適正速度」は全く別物です。

落とし穴2：「均一球体」という仮定

シミュレーターは完全な球体を仮定しますが、現実の粒子は：

  • 粒径分布を持つ（例：50〜300μmまでバラバラ）
  • 形状がいびつ（球形度 $\phi_s < 0.7$）
  • 湿潤時に凝集する

このため、計算値はあくまでオーダー推定。実機設計では、サウター平均径やメディアン径を使い、必ず実験データと突き合わせます。

落とし穴3：圧力損失の「一定値」神話

「流動化後、圧力損失は一定」と教科書には書いてあります。しかし、大型装置では：

  • 気泡の発生による圧力変動
  • 粒子の偏析（重い粒子は下、軽い粒子は上）
  • チャネリング（ガスが一部の経路だけを通る）

これらにより、圧力損失は常に変動します。シミュレーターの美しいカーブは理論的な骨格であり、実機では計測値とのフィッティングが不可欠です。

まとめ：流動層は「見える化」が全て

今回のポイントを3つに絞ります。

  1. 流動化の本質は「抗力 vs 重量」の釣り合い。アルキメデス数がその指標。
  2. Ergun式 + Wen & Yu 相関式で、Umf・Ut・圧力損失が計算できる。コードは50行。
  3. 理論値は「目安」。実機設計では2〜10倍の安全率、粒度分布、凝集を考慮する。

そして何より、この流動層シミュレーターは、「粒子が浮く瞬間」をリアルタイムで可視化してくれる稀有なツールです。パラメータを動かしながら、数式の意味を体感してみてください。

「百聞は一見に如かず」 ─ 工学も同じ。数式をコードで動かし、シミュレーターで確かめる。このサイクルこそが、真の理解への近道です。

▶ 流動層シミュレーター — ブラウザで即動作、登録不要
NovaSolverでは700以上の工学シミュレーターを無料公開中 👉 一覧はこちら

0
0
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
0
0

Delete article

Deleted articles cannot be recovered.

Draft of this article would be also deleted.

Are you sure you want to delete this article?