（様式 20）成果報告書

TSUBAME 共同利用 平成 25 年度 学術利用 成果報告書

利用課題名 多相流 LBM シミュレーションの大規模間隙モデルへの適用
英文： Multiphase LBM simulation for large‐scale porous media
辻 健
Takeshi TSUJI
九州大学 カーボンニュートラル・エネルギー国際研究所 (I2CNER)
International Institute for Carbon‐Neutral Energy Research (I2CNER),
Kyushu University
http://i2cner.kyushu‐u.ac.jp/~tsuji/
格子ボルツマン法の計算コードを、Graphics Processing Unit（GPU）に実装し、岩石の間隙内部の CO2の挙動を計算した。GPU を用いることにより、マイクロ X-CT によって取得された巨大なデジタルロック（グリッドサイズ：10003）に対して、多相流 LBM シミュレーションを行うことが可能となった。計算領域のサイズアップを行うことで、これまで一般的であったポアスケール（μm スケール）から、コアスケール（mmスケール）へ計算領域を拡張することができた。それにより、LBM シミュレーション結果を、実験結果（コアスケール）と比較することが可能となった。
We investigated fluid flow characteristics inside 3D porous media by a highly efficient multi-phase lattice
Boltzmann method (LBM). By using a graphics processing unit (GPU), we can conduct two-phase LBM simulation forthe largest digital rocks (grid size) in the world (10003). Multi-phase flows were simulated for the digital rock model at
extremely high resolution (~2 μm) reconstructed from micro-CT scanned image. These efforts enable us to enlarge thesize of digital rock from pore-scale (μm scale) to core-scale (mm scale). Therefore, we can validate the LBM simulation
results with laboratory experiments (i.e., core-scale measurements).

Keywords: 多相流 LBM，GPGPU，世界最高サイズのデジタルロック，二酸化炭素の地中貯留(CCS)

背景と目的多孔質形状内の流れの計算に適している。一方で、二酸化炭素の地中貯留（CCS）は、岩石の間隙に計算機やイメージング技術などの制限から、従来はCO2 を注入することで永久的に地下に閉じ込める技

Finney pack モデルといった単純な間隙モデルに対術で、大気中の CO2 の削減に貢献できる技術としてして、LBM シミュレーションが適用される場合が多注目されている。CCS では、岩石の間隙内部の CO2かった（例えば Tsuji et al. 2012）。の挙動を調べることが重要であり、多くの実験が実

近年、Micro-CT などのイメージング技術が発達に施されている。さらに流体シミュレーションを用いより、実際の岩石の間隙形状（デジタルロック）をることで、CO2 を圧入する貯留層の条件（例えば界抽出することができるようになった。さらに計算機

面張力や粘性）を変えて、圧入 CO2 の挙動を推定すの発達によって、デジタルロックに対して、LBM 流ることが試みられている。また CCS だけでなく、体シミュレーションを実施できるようになってきた。

Enhanced Oil Recovery（EOR）といった幅広い分

しかし多くの LBM シミュレーションでは、計算の野で、間隙内の流体挙動の定量化に向けた研究が行制限から、低解像度かつ小さいスケールで実施されわれている。

ることが多かった。例えば、Ahrenholz et al. (2008)

格子ボルツマン法（Lattice Boltzmann Method;

は、11μm の解像度で 2003 グリッドの砂岩モデルにLBM）は、速度の分布関数の時間発展方程式を解く

対して、Capillary pressure を計算している。また近ことで、流体運動を記述する手法である(Chen and年の論文（例えば Ramstad et al. 2012)でも、6μm

Doolen, 1998)。LBM は圧力ポアソン方程式を解かの解像度で 2563 グリッドのべレア砂岩などのデジず、また局所的な計算がメインであるため、複雑なタルロックに対して、LBM を実施している。

（様式 20）成果報告書
本研究では、Graphics Processing Unit（GPU）をで管理する。並列演算可能な最少単位をスレッド、用いた並列計算を行うことによって、実際の岩石のスレッドの集合をブロックと呼ばれる。CUDA 関数

間隙モデルに対して、高解像度（2 μm）かつ大規模（カーネル）とは GPU で動くプログラムを指し、3（グリッドサイズ：1000 ）に LBM シミュレーショどれくらいのスレッド数で「カーネル関数」を実行ンを行うことを可能とした。このような巨大な間隙さ せ る か 指 定 す る こ と が で き る 。 SPMD(Singleモデルに対して、2 相の LBM シミュレーションを実Program Multiple Data)と呼ばれるプログラミング施した例は、（現時点では）世界でも初めてであるとモデルが採用されているため、実行するプログラム考えられる。このように計算領域のサイズアップをは全てのスレッドで同じである。ここで、各スレッ行うことで、ポアスケール（μm スケール）からコドがどの格子の計算を担当するのを指定する必要がアスケール（mm スケール）の数値シミュレーショある。すなわち、スレッド ID と更新される格子の空ンを可能となり、実験室で得られた結果と比較した間座標を示さなければならない。我々は、図１に示議論を行うことができるようになった。最終的には、すようにｚ方向の一列の格子点の更新を１スレッド貯留層の条件を変えた LBM シミュレーションを実で行うことにした。１ブロックが担当する計算領域施し、CO2 圧入に適した貯留層の条件などを明らかは、blockDimx, blockDimy のパラメータによって決にすることを目指す。

められる。図１のようにブロックをｘ、ｙ平面にマッピングさせる。つまり、Nx、Ny を計算領域ｘ、概要

ｙ方向のメッシュ数とすると、ｘ、ｙ方向のブロッ

GPU は、従来は PC やワークステーションにおいク数は、それぞれ Nx/blockDimx, Ny/blockDimy であて画像処理を担当する主要な部品の一つとして使わる。これによって、関数内のｘ、ｙループが必要なれてきた。しかし画像処理を行うだけでなく、GPUくなり、GPU スケジューラが自動的に複数のコアをに内蔵している大量のコア数を生かし、汎用的な並利用し、並列計算を行うことができる。スレッド、列 計 算 に 応 用 す る GPGPU （ General Purpose

ブロック ID

Computation on Graphics Processing Unit）、あるいblockID.y）と 格子座標（jx, jy）の関係は、次式のは「GPU コンピューティング」と呼ばれる技術が、通りである。

（threadID.x, threadID.y, blockID.x,

注目を集めてきた。この GPU を用いた計算技術によって、低価格かつ高演算能力のシミュレータが実jx  blockDim . x  blockIdx . x  threadIdx . x ;(1)jy  blockDim . y  blockIdx . y  threadIdx . y .

(2)

現できるようになってきた。
我々は LBM の計算コードを GPU に実装し、高効率シミュレータの開発を行ってきた。LBM では、局所的な計算がメインとなるため、大規模 GPU 並列

GPU メモリ内確保されている配列のデータレイ計算に適していると考えられる。特に、今回採用すアウトも、計算を効率的に行う上で重要である。

る RK color 2 相流 LBM モデル（Tölke et al., 2006）

GPU のメモリバンド幅を無駄なく活用するため、コは、界面の追跡にフェーズフィールド法の概念を取アレッシングというデータまとめて転送する方式をり込み、並列化が容易で、且つメモリ使用量は従来利用することが重要である。ここで LBM の分布関の方法よりも節約できる。そのため、この LBM モ数は、以下のレイアウトで保存される。

デルは、GPU 実装に適していると考えられる。
マルチ GPU 実装は、２段階に分けて行う。シンi  L z  L y  L x  jz  L y  L x  jy  L x  jx

（3）

グル GPU チップ内部の並列計算アルゴリズムは、マルチスレッド処理と呼ばれる。CUDA はソフトウ

またデータアクセス量を減らすため、衝突過程と移ェア上で、大量のコアをスレッド、ブロックの構造流過程を合併した。これにより、計算効率を 30％向（様式 20）成果報告書上することができた。

一方で異なるノード間のデータのやり取りは MPI（Message Passing Interface）規格を準じて、LANケーブルまたは Infiniband を経由して実施する。データ転送を効率するために、CPU ホスト側のメモリ保存は pinned メモリ方式で行う。このようなコードの最適化によって、二相流 LBM シミュレーションの計算効率を大幅に向上させることに成功した。
以上のように LBM の計算コードを効率化することで、10003 グリッドという大きな領域で、CO2 圧入シミュレーションを実現した。この計算によって、

図１．GPU 内部のマルチスレッド化（Jiang et al.,in review）.

これまでポアスケールで行ってきた LBM 計算を、コアスケールに拡張することができた。

結果および考察GPU を用いることで、べレア砂岩のマイクロX-CT 画像から構築した解像度 2μm、格子サイズ10003 のデジタルロック（図３）に対して、2 相 LBMシミュレーションを実施することに成功した。このデジタルロックの実際のサイズは、2 ㎜に相当し、
図２．マルチ GPU 間のデータ転送モデル（Jianget al., in review）.

実験データとの比較が可能となってきた。
図４は、Drainage 過程の CO2 圧入シミュレーションを行った例である。計算条件として、界面張力大規模計算を行う際、一枚の GPU メモリが最大を 60mN/m 、 超 臨 界 CO2 密 度 は 水 と 同 じ 値６GB の制限があるため、メモリ使用量を出来る限1000kg/m3 に 設 定 し た 。 ま た 粘 性 は 水 と 同 じり減らす必要がある。我々は、メモリを節約するた

1.0×10 6m2/s に固定した。流れを駆動するため、め、次の２つ方法を用いた。一つ目は計算する必要体積力を与え、入口と出口は周期条件に設定した。のない固体メッシュの情報を保存しない方法、二つ今回の計算では、CO2 は non-wetting な特徴を持ち、目はストリームパタンを工夫し、一つだけの配列に壁は wettable な特徴を持つと仮定した。

分布関数を保存する方法である。この二つの方法を

図４の計算では表面張力を低く設定しているため、用いることで、10 億 (1000*1000*1000) 格子の大規CO2 が小さい間隙に侵入していることが分かる。今

模シミュレーションを、20 枚の GPU（Tesla-K20）

回の計算から、Drainage 過程では、CO2 飽和度はで計算することが可能となった。

圧入時間と共に 50%程度まで増加することが分か

マルチ GPU 計算では、はじめに領域を分割するった。この LBM シミュレーションで得られた CO2必要がある。本研究では、流体格子数を均等に割り、の飽和度は、実験室で得られる値と整合的である。

ｚ方向に分割する方法を採用した。分割領域界面でこのような間隙内での CO2 の挙動から、絶対浸透は、一層の ghost レイヤーが設けられ、これは領域率や相対浸透率を計算することも可能になった。こ間のデータ転送かつカラー勾配の計算に用いる。

れらの水理特性は、貯留層シミュレーションによっGPU 間のデータ転送は、二つのパターンがある（図て、将来の圧 CO2 の分布を予測する上で、重要なパ２）。同一のマザーボード（ノード）に装着しているラメータである。

GPU ボード間のデータ転送は、CPU を介して行う。

（様式 20）成果報告書
し本研究で行ったように、GPU を駆使することで、ポアスケールとコアスケールの壁を取り除くことができた。
今後は、本研究で開発した手法（巨大間隙モデル
への LBM シミュレーション）を用いて、様々な貯留層の条件（例えば界面張力や粘性）で、CO2 の挙動を推定する。それにより、間隙内での CO2 挙動の定量化を行い、貯留層シミュレーションに組み込むことを目指す。

図３．マイクロ X-CT 画像から構築されたべレア砂岩のデジタルロック．

本研究は、地球環境産業時術研究機構（RITE）との共同研究で実施したものである。

引用文献

図４．Drainage 過程のシミュレーション結果（グリッドサイズ：10003）．白い部分が間隙内に圧入したCO2（Non-wettable phase）である．右側から左側に向かって CO2 を圧入している．

まとめ、今後の課題
GPU を用いることにより、マイクロ X-CT によって取得された巨大なデジタルロック（グリッドサイズ：10003）に対して、LBM シミュレーションを実施することが可能となった。それにより、これまでポアスケールで行うことの多かった LBM シミュレーションを、コアスケールで実施することが可能となった。スケールの違いを考慮するためには、異なった手法を用いることが多く、それらの現象を統一的に理解することは難しいと考えられている。しかAhrenholz. B, J. Tölke, P. Lehmann, A. Peters, A.
Kaestner, M. Krafczyk, W. Durner, “Prediction ofcapillary hysteresis in a porous material using
lattice-Boltzmann methods and comparison toexperimental data and a morphological pore
network model”, Advances in Water Resources, 31(9), 1151–1173, 2008.
Chen, S., and G.D. Doolen, "Lattice Boltzmann Method
for Fluid Flows", Annual Review of Fluid
Mechanics, 30, 329-364, 1998.
Jiang, F., T. Tsuji, and C. Hu, “Elucidating the role ofinterfacial tension for hydrological properties oftwo-phase flow in natural sandstone by an
improved lattice Boltzmann method”. in revision.
Ramstad. T, N. Idowu, C. Nardi, P.E. Øren, “Relative
permeability calculations from two-phase flow
simulations directly on digital images of porousrocks”, Transport in Porous Media, 94 (2),487-504, 2012.
Tölke, J., et al., “An adaptive scheme using
hierarchical grids for lattice Boltzmann
multi-phase flow simulations”, Computers &Fluids, 35, 820-830, 2006.
Tsuji, T., H. Yamabe, and T. Matsuoka, “LBM
simulation for CO2 saturation monitoring fromelastic velocity and resistivity: Migration ofsupercritical CO2 in porous media under several
PT conditions”, AGU fall meeting, GC51A-1176,
San Francisco, December 3-7, 2012.