析出アスファルテンの孔隙閉塞数値モデルの構築
村田澄彦*1
1.

研究の目的

原油中のアスファルテンが油層内や坑井内で析出すると、油層の孔隙やチュービングパイプが閉塞し、重大な生産障害となる。この問題の対策を検討するため、アスファルテン析出の重要な指標の 1 つである AOP（Asphaltene Onset Pressure）の評価およびアスファルテン析出挙動の観察が FLASS（Flow Assurance System)を用いてなされている 1)。しかし、様々な因子の影響を受けていると考えられる複雑なアスファルテンの析出挙動を閉塞に至るまで FLASS で検討するには困難がある。
そこで、本研究では、アスファルテン粒子同士の相互作用を考慮して油層孔隙内のアスファルテン流動場におけるアスファルテン粒子の凝集とそれによる孔隙の閉塞過程をシミュレートできる計算手法について検討し、その手法を用いてアスファルテン析出に関する各因子の影響を個別に考察できるシミュレーションモデルを構築することを目的とする。これにより、アスファルテン析出による生産障害に対する対策検討に資することができると考える。
2.

研究の方法

2.1 シミュレーション手法の検討
アスファルテン粒子同士の相互作用を考慮して、油層孔隙内のアスファルテン流動場におけるアスファルテン粒子の凝集とそれによる孔隙の閉塞過程をシミュレーションするには、① アスファルテンの流動場、② アスファルテンの凝集、③ アスファルテンの固体壁面への吸着・固着、④ アスファルテン以外の分子との相互作用を考慮する必要がある。
油層の孔隙のような複雑な境界面を有する空間内の流体流動計算には、空間を細かく分割した立方体格子（2 次元では正方形格子）上における仮想粒子の単純な並進演算と衝突演算の繰り返しでNavier-Stokes 方程式を解くことができる格子ボルツマン法（LBM: Lattice Boltzmann Method）が良く用いられている。また、LBM に対して流体と固体の相互作用を考慮して流動場における固体の流動を計算するための衝突演算項が提案され 2)、LBM と個別要素法（DEM: Distinct Element Method）との連成解析手法も提案されている 3)。しかしながら、この LBM-DEM 連成解析手法を本研究に適用しようとすると、流動場に対して十分に小さいアスファルテン粒子よりさらに細かい LBM 格子を用いる必要があり、計算コストが非常に高くなって非現実的であると考えられる。
これに対して、確率的回転動力学法（SRD: Stochastic Rotation Dynamics）4)は計算格子を用いない方法であり、後に Ihle と Kroll5)が指摘した任意の状態でガリレイ不変性を保証するためのグリッドシフト手順を含めると溶媒粒子の単純な並進演算と衝突演算により Navier-Stokes 方程式に従う流体流動を正しく計算できる他、溶媒粒子と微小固体粒子の相互作用ポテンシャルおよび微小固体粒子間の相互作用ポテンシャルを考慮することで、微小固体粒子のブラウン運動も正しく計算することができる。Boek ら 6)は SRD を用いて毛細管内のアスファルテン粒子の流動と凝集による毛細管閉塞のシミュレーションを行い、顕微鏡下の観察結果とよく一致することを報告している。そこで、本研究では Boek らに従って SRD を用いて目的とするアスファルテン析出に関する各因子の影響を個別に考察できるシミュレーションモデルを構築することとする。その初期段階として SRD を用いてアスファルテン粒子を含まない状態の平行平板内の単純な 2 次元流体流動計算を実施した。以下に SRD による流体流動計算の概要を示し、その結果について報告する。
2.2 SRD による流体流動計算の概要
SRD は格子を用いない計算手法ではあるが、流体の流速、密度、温度、粘度などのマクロな物理量を局所的に計算するため、計算領域をセルに分割する。なお、各セルに含まれる溶媒粒子の数に制限はない。流体流動は、溶媒粒子の単純な並進演算と衝突演算によって計算され、それぞれ次式で表される。

*

京都大学・大学院工学研究科・准教授

𝐟ex 2
𝛿𝑡 ,2𝑚𝑖並進演算：

𝐫𝑖 (𝑡 + 𝛿𝑡) = 𝐫𝑖 (𝑡) + 𝐯𝑖 (𝑡)𝛿𝑡 +衝突演算：

𝐯𝑖 (𝑡 + 𝛿𝑡) = 𝐮′(𝑡 + 𝛿𝑡) + 𝐑(𝐯′𝑖 (𝑡 + 𝛿𝑡) − 𝐮′(𝑡 + 𝛿𝑡))𝐯′𝑖 (𝑡 + 𝛿𝑡) = 𝐯𝑖 (𝑡) +𝐟ex𝛿𝑡𝑚𝑖(1)(2)

ここで、𝐫𝑖 (𝑡)と𝐯𝑖 (𝑡)はそれぞれ溶媒粒子𝑖の時刻𝑡における位置ベクトルと速度ベクトルであり、𝛿𝑡は時間ステップ間隔、𝑚𝑖 は溶媒粒子𝑖の質量、𝐟ex は全溶媒粒子に作用する外力ベクトルである。また、𝐯′𝑖 (𝑡 + 𝛿𝑡)は衝突演算前の溶媒粒子𝑖の速度、𝐮′ (𝑡 + 𝛿𝑡)は𝐯′𝑖 (𝑡 + 𝛿𝑡)に対する各セルの重心速度、
𝐑はタイムステップ毎に各セルに設定するランダムな方向の回転軸周りの角度𝛼の回転行列である。
なお、溶媒粒子の平均自由行程がセルサイズ𝑎0 より小さくなるとガリレイ不変性が保証されなくなる。この問題を解決するため、衝突演算前に±𝑎0 /2以下の成分を持つランダムベクトルで溶媒粒子全体をシフトし、衝突演算後に同じ量だけシフトバックするグリッドシフトを行う。また、固体
壁面における境界条件にはいくつかあるが、本研究では固体壁面をサーモスタットとして機能させる確率論的反射境界条件 7)を採用し、流体の流入流出口に対しては周期境界条件を設定する。
3.

得られた成果9.0

9.08.08.07.0

7.0

6.0

6.0

5.00

4.020ky = -0.051x2 + 1.2763x + 0.3073

3.0100k

2.0多項式 (20k)1.0

多項式 (100k)

R² = 0.9951

3.0

1.00.0
15Y-position

図1謝

100k

2.0

y = -0.052x2 + 1.3032x + 0.1739
R² = 0.9965

-1.0
1020k

4.0-1.0
505.0

0.00

4.

Y-velocity

X-velocity

平行平板間の 2 次元断面モデルとして、𝑥方向に 150、𝑦方向に 25 の正方形セルからなる計算領域を設定し、静止した溶媒粒子を各セル 5 個ずつランダムな配置で含む状態から+𝑥方向に一様な物体力を外力として作用させ+𝑥方向に流体を流動させた。𝑦 = 0と𝑦 = 25の位置が固体壁面となっている。このときの、時間ステップが 0、20,000（20k）、100,000（100k）の𝑥方向と𝑦方向それぞれの流速分布を図 1 に示す。図 1 に示したとおり、固体壁面の流速がほぼ 0 であり、𝑥方向の流速分布は放物線で近似でき、𝑦方向の流速はほぼ 0 となっていることから、作成した SRD のシミュレーションモデルで流体流動が正しく計算できているものと考えられる。今後、アスファルテン粒子の流動と凝集をこのシミュレーションモデルに組み込む予定である。

20250510152025Y-position

平行平板間の 2 次元断面モデルにおける流速分布

辞

本研究は、国際石油開発帝石(株)により委託されたものである。関係各位の協力と助言に厚く感謝申し上げます。
参 考 文 献1) 宮川善洋, 渡辺拓己, 米林英治, 石油技術協会誌, 第 83 巻 第 3 号, pp. 201-213, 20182) D. Noble and J. Torczynski, Int. J. Modern Phys. C, 9(08), pp. 1189-1201, 1988
3) S.A. Galindo-Torres, Comput. Methods Appl. Mech. Engrg. 265, pp. 107–119, 2013
4) A. Malevanets and R. Kapral, J. Chem. Phys. 110(17), pp. 8605-8613, 1999
5) T. Ihle and D. M. Kroll, Physical Review E, 63(2), pp. 020201-1-4, 2001
6) E. S. Boek, H. K. Ladva, J. P. Crawshaw and J. T. Padding, Energy & Fuels, 22(2), pp. 805-813, 2008
7) Y. Inoue, Y. Chen and H. Ohashi, Journal of Statistical Physics, 107(1/2), pp. 85-100, 2002