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

はじめに

本記事では、材料の特性評価における流体力学シミュレーションの応用について解説します。X線CTなどで取得したデータをもとに、デジタル空間上で流体の透過性をシミュレーションする方法を紹介し、その理論的背景、計算手順、および得られる洞察について詳しく説明します。

この記事を書いたひと

デジタルリアクタ合同会社 代表
望月 優輝 / Yuki Mochizuki   

機械学習・統計、数値計算などの領域を軸としたソフトウェアエンジニアリングを専門としています。スタートアップからグローバル企業まで、さまざまなスケールの企業にて、事業価値に直結する計算システムを設計・構築してきました。主に機械学習の応用分野において、出版・発表、特許取得等の実績があり、また、IT戦略やデータベース技術等の情報処理に関する専門領域の国家資格を複数保有しています。九州大学理学部卒業。日本ITストラテジスト協会正会員。

対象読者:

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

記事のポイント:

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

処理の手順

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

シミュレーションの精度は、セグメンテーション画像の品質に大きく左右されます。X線CTやMRIなどの技術で得られた3次元画像データを、材料の固体部分と空隙部分に正確に分離することが重要です。このセグメンテーション処理により、材料内部の微細構造が明確になり、流体の流れをより正確にシミュレートできます。

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

理論

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

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

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

\frac{\partial \mathbf{u}}{\partial t} = -\frac{1}{\rho} \nabla p + \nu \nabla^2 \mathbf{u} + \mathbf{a_f}

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

\nabla \cdot \mathbf{u} = 0

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

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

境界条件

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

収束条件

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

出力

シミュレーションの主な出力は、定常状態における流体の速度場と圧力分布のデータです。これらのデータを解析することで、材料の特性を詳細に評価できます。

結果の分析

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

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

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

k = \frac{Q \mu L}{S \Delta P}

• Q は流量 (m³/s)
• \mu は流体の粘性率 (Pa・s)
• L は流れ方向のサンプルの長さ (m)
• S はサンプルの断面積 (m²)
• \Delta P はサンプル両端の圧力差 (Pa)

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

計算例:

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

k = \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} \sim 10^{-10} m^2 程度、粘土は 10^{-11} \sim 10^{-15} m^2 程度です。

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

固有浸透テンソルは、空間の任意の方向における浸透率の情報を提供し、多孔質媒体の異方性を評価できます。

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

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

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

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

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

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

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

曲路率 (Tortuosity) の計算

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

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

tortuosity = \frac{\sum |\mathbf{u}|}{\sum \mathbf{u} \cdot \mathbf{x} }

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

典型的な流路の可視化

シミュレーション結果から得られる流体の典型的な流路を可視化することで、材料内部の流れの挙動を直感的に理解できます。これにより、材料の特性、欠陥、改良点などを視覚的に確認できます。可視化には、計算に用いたソフトウェアや、ParaViewなどの専用のオープンソースソフトウェアが適しています。

まとめ

本記事では、流体力学シミュレーションを用いて材料の透過性評価を行う方法について解説しました。シミュレーションにより、絶対浸透率、固有浸透テンソル、曲路率などの指標を計算し、材料の異方性を評価できることを示しました。

流体力学シミュレーションは、環境工学や医療分野など、幅広い材料科学分野で応用が期待されています。特に、微細構造の解析や新材料の開発において、迅速かつ再現性の高い評価ツールとして重要性を増していくでしょう。

関連投稿:

確率的最適化アルゴリズムの実装と比較：遺伝的アルゴリズムと勾配降下法 準モンテカルロ法で使用される低不一致列の解説、次元数による有効性の違い 次元削減手法の比較検証：定量的・定性的評価によるPCA、t-SNE、UMAPの特徴分析 数値積分の実装と比較：台形法からガウス求積法まで