Digital Reactor
数値計算・最適化

流体力学シミュレーションによる材料構造評価 - 絶対浸透率と固有浸透テンソル、曲路率

流体力学シミュレーションによる材料構造評価 - 絶対浸透率と固有浸透テンソル、曲路率

はじめに

多孔質材料の透水性や通気性を実物の透過試験で測ろうとすると、試料の切り出しや装置の準備に手間がかかるうえ、同じ試料で条件を変えて測り直すのも容易ではありません。一方、X線CTで内部構造の3次元データが手元にあれば、その画像の上で流体シミュレーションを実行し、絶対浸透率・固有浸透テンソル・曲路率といった指標を数値として取り出せます。この記事では、その計算手順と背後にある理論、結果の読み方を順に説明します。

対象読者:

  • 材料科学の研究者・技術者
  • 流体力学シミュレーションに興味がある方
  • 多孔質材料の特性評価に関心がある方

記事のポイント:

  • 流体力学シミュレーションを用いた材料の透過性評価の原理と手順を理解できる
  • 絶対浸透率、固有浸透テンソル、曲路率といった指標の計算方法と、それらが示す材料特性について学べる
  • シミュレーション結果の可視化による、材料内部の流れの挙動の理解を深められる

処理の手順

セグメンテーション画像の準備

シミュレーションの精度は、入力となるセグメンテーション画像の品質でほぼ決まります。X線CTやMRIで得た3次元画像データを、材料の固体部分と空隙部分に分離しておきます。この分離が粗いと空隙のつながり方が実物と変わってしまい、流れの経路が変わるため、後段で求める浸透率にそのまま影響します。

流体力学シミュレーションの実行

支配方程式と仮定

絶対浸透率を数値的に推定するには、ナビエ・ストークス方程式を解きます。粘性率が一定の流れは以下の式で記述できます。

∂u∂t+(u⋅∇)u=−1ρ∇p+ν∇2u+af\frac{\partial \mathbf{u}}{\partial t} + (\mathbf{u} \cdot \nabla) \mathbf{u} = -\frac{1}{\rho} \nabla p + \nu \nabla^2 \mathbf{u} + \mathbf{a_f}

ここで、

  • ∇\nabla はナブラ演算子（勾配）
  • ∇2\nabla^2 はラプラシアン演算子（勾配の発散）
  • u\mathbf{u} は流体の速度 (m/s)
  • ν\nu は流体の動粘性率 (m²/s) （粘性率 (Pa・s) / 密度 (kg/m³)）
  • pp は流体の圧力 (Pa)
  • ρ\rho は流体の密度 (kg/m³)
  • af\mathbf{a_f} は外力による加速度 (m/s²)

流体の速度が遅く、レイノルズ数が小さい場合、非線型項である対流項 (u⋅∇)u(\mathbf{u} \cdot \nabla) \mathbf{u} は無視でき、ストークス方程式が適用できます。

∂u∂t=−1ρ∇p+ν∇2u+af\frac{\partial \mathbf{u}}{\partial t} = -\frac{1}{\rho} \nabla p + \nu \nabla^2 \mathbf{u} + \mathbf{a_f}

非圧縮性流体を扱うため、以下の制約も考慮します。

∇⋅u=0\nabla \cdot \mathbf{u} = 0

材料の浸透性評価では、多くの場合、流速が遅いと仮定してストークス方程式を用います。シミュレーションでは、以下の条件を想定します。

  • 非圧縮性流体：流体の密度は一定
  • ニュートン流体：流体の動粘性率は一定
  • 層流：流れは層状であり、乱流は発生しない

境界条件

周期的境界条件を適用すると、固有浸透テンソルを計算できます。フーリエ変換を利用した解法も使えるようになり、計算を大幅に高速化できます。周期的境界条件では、固定の圧力や流量を直接設定せず、流体の全粒子に均等な力を加えて流れを発生させます。

周期境界条件: 3次元的な繰り返し構造を仮想的に生成してシミュレーションを実行する

収束条件

流れの速度が時間経過とともに変化しなくなる（定常状態になる）まで計算を続けます。流れの速度や圧力の変化が一定範囲内に収まるまで計算を繰り返し、残差の大きさや収束の速さなどを基に収束を判定します。

出力

シミュレーションの主な出力は、定常状態での流体の速度場と圧力分布のデータです。以降の分析では、この2つから各指標を計算します。

速度場と圧力分布から材料特性を求める

絶対浸透率 (Absolute Permeability) の計算

絶対浸透率は、多孔質材料が単相流体を透過させる能力を示す指標です。SI単位は平方メートル (m²) ですが、平方マイクロメートル (μm²) やダルシー (d) もよく用いられます（1ダルシー ≈ 9.869233 × 10⁻¹³ m² ≈ 9.869233 μm²）。絶対浸透率は、材料固有の特性であり、外部条件に依存しません。

絶対浸透率 kk (m²) は、ダルシーの法則に基づき、以下の式で計算されます。

k=QμLSΔPk = \frac{Q \mu L}{S \Delta P}
  • QQ は流量 (m³/s)
  • μ\mu は流体の粘性率 (Pa・s)
  • LL は流れ方向のサンプルの長さ (m)
  • SS はサンプルの断面積 (m²)
  • ΔP\Delta P はサンプル両端の圧力差 (Pa)

周期的境界条件を適用した場合、明示的な圧力差は設定されません。一定の外力をかけ続けたときに、流量が収束する値を観察して、ダルシーの法則から計算します。 ここで、外力と圧力差の関係について説明します。 簡単のため、流体として水を想定します。流体部分に外力による加速度として、重力加速度と同じ 9.8m/s29.8 m/s^2 が与えられている時を考えます。 水の密度は、1000kg/m31000 kg/m^3ですので、 高さ1m1 mあたり、面積1m21 m^2あたり、9800kg⋅m/s2=9800N9800 kg \cdot m/s^2 = 9800 Nの外力がかかっていることになります。 1Pa=1N/m21 Pa = 1 N / m^2ですので、高さ1m1 mにつき、9800Pa9800 Paの圧力差がある場合に相当することになります。

ダルシーの式は圧力差を必要としますが、周期境界条件を適用したシミュレーションでは圧力差ではなく外力で直接流体を駆動しています。通常は「圧力差があるから流体が加速する」と考えるところを、ここでは向きを逆にして「この外力で加速するなら、どれだけの圧力差を与えたときと同じ力を受けているか」と読み替えます。

計算例:

1辺0.1 mの立方体形状の多孔質サンプルに水を浸透させ、重力加速度をかけたとします。 この時、ΔP=9800Pa\Delta P = 9800 Pa となります。 定常状態での流量が Q=50ml/s=0.00005m3/sQ = 50 ml/s = 0.00005 m^3/s であった場合、水の粘性率 0.001Pa⋅s0.001 Pa \cdot s を用いて、絶対浸透率 kk は以下のように計算できます。

k=(0.00005m3/s)×(0.001Pa⋅s)×(0.1m)(0.01m2)×(9800Pa)=5.10×10−11m2k = \frac{(0.00005 m^3/s) \times (0.001 Pa \cdot s) \times (0.1 m)}{(0.01 m^2) \times (9800 Pa)} = 5.10 \times 10^{-11} m^2

参考：砂利の絶対浸透率は 10−7∼10−10m210^{-7} \sim 10^{-10} m^2 程度、粘土は 10−11∼10−15m210^{-11} \sim 10^{-15} m^2 程度です。

固有浸透テンソル (Intrinsic Permeability Tensor) の計算

固有浸透テンソルを使うと、空間の任意の方向の浸透率が分かり、多孔質媒体の異方性を評価できます。

テンソルは、スカラー、ベクトル、行列を一般化した概念です。3次元空間における多孔質材料の浸透率は、以下の3x3行列（浸透テンソル）で表されます。

K=(kxxkxykxzkyxkyykyzkzxkzykzz)\mathbf{K} = \begin{pmatrix} k_{xx} & k_{xy} & k_{xz} \\ k_{yx} & k_{yy} & k_{yz} \\ k_{zx} & k_{zy} & k_{zz} \end{pmatrix}

ここで、kijk_{ij} は、i方向の圧力勾配がj方向の流れを生成する効果を表します。

固有浸透テンソルを計算するには、3つの直交する方向（x, y, z）で絶対浸透率を計測します。つまり、外力の方向を3パターン用意し、3回のシミュレーションを行います。kijk_{ij} (i≠j) を計算する場合、流量Qには、外力の方向とは異なる方向の流量を指定します。例えば、kxyk_{xy} を計算するには、外力の方向を (x, y, z) = (1, 0, 0) とし、y平面を通過する流量を用います。

正しく計算されていれば、この行列は対称性を持ち、kij=kjik_{ij} = k_{ji} となります。

次に、多孔質媒体の異方性を評価するため、主成分分析 (PCA) を用いて固有ベクトルと固有値を計算します。これにより、テンソルが対角化され、主成分方向とそれに対応する浸透率が特定できます。

固有ベクトルは材料内の主要な流れの方向を示し、固有値は各固有ベクトルに対応する浸透率を表します。これらの情報から、透過性に関する異方性を評価できます。例えば、3つの固有値が大きく異なる場合、透過性に強い異方性があることを意味します。これは、材料内部に特定の方向に沿った繊維や平板構造が存在する場合によく見られます。

曲路率 (Tortuosity) の計算

曲路率は、流体が材料内を通過する際の経路の複雑さを示す指標です。材料の多孔質構造が流体の流れをどれだけ妨げるかを定量的に評価します。

曲路率は、流路の実際の長さと直線距離の比として定義されることもありますが、シミュレーション結果からは、以下の式で計算することもできます。

tortuosity=∑∣u∣∑u⋅xtortuosity = \frac{\sum |\mathbf{u}|}{\sum \mathbf{u} \cdot \mathbf{x} }

ここで、分子は流体部分の速度ベクトルの大きさの合計を示します。分母は、関心のある軸方向x\mathbf{x}へ速度ベクトルを射影したものの合計です。

典型的な流路の可視化

シミュレーション結果から流体の典型的な流路を可視化すると、材料内部の流れの挙動を直感的につかめます。流れが集中している箇所や、ほとんど流れに寄与していない空隙が見えるため、材料の欠陥や改良点の検討にも役立ちます。可視化には、計算に用いたソフトウェアのほか、ParaViewのようなオープンソースソフトウェアが使えます。

ParaViewによる可視化: 速度が大きいところの色を変えて表示している

まとめ

セグメンテーション画像さえ用意できれば、流体シミュレーションで絶対浸透率・固有浸透テンソル・曲路率を数値として取り出せます。押さえておくべき点は、周期的境界条件では圧力差を直接与えないため、外力を等価な圧力差に読み替えてダルシーの法則に当てはめるところです。あとは定常状態まで計算を進め、得られた速度場を集計するだけです。

実物の透過試験と違って、同じ試料で外力の方向を変えた再計算が容易なので、異方性の評価まで一度の準備で進められます。手元にCTデータがあるなら、まず固体と空隙の分離を丁寧に行い、1方向の絶対浸透率を計算して、砂利や粘土の文献値と桁が合うかを確かめるのが最初の確認になります。

関連記事

← 技術ブログ一覧へ