特許 7725111 ３次元粗い離散亀裂ネットワークの構築、評価及び浸透流シミュレーションの方法

(19)【発行国】日本国特許庁(JP)

(12)【公報種別】特許公報(B2)

(11)【特許番号】

(24)【登録日】2025-08-08

(45)【発行日】2025-08-19

(54)【発明の名称】３次元粗い離散亀裂ネットワークの構築、評価及び浸透流シミュレーションの方法

(51)【国際特許分類】

   G06F 30/20 20200101AFI20250812BHJP

   G06F 30/28 20200101ALI20250812BHJP

【ＦＩ】

G06F30/20

G06F30/28

【請求項の数】 5

(21)【出願番号】P 2024207004

(22)【出願日】2024-11-28

(65)【公開番号】P2025088759

(43)【公開日】2025-06-11

【審査請求日】2024-11-28

(31)【優先権主張番号】202311617961.2

(32)【優先日】2023-11-30

(33)【優先権主張国・地域又は機関】CN

【早期審査対象出願】

(73)【特許権者】

【識別番号】523347497

【氏名又は名称】山東科技大学

(74)【代理人】

【識別番号】110000291

【氏名又は名称】弁理士法人コスモス国際特許商標事務所

(72)【発明者】

【氏名】韓 承豪

(72)【発明者】

【氏名】陳 紹杰

(72)【発明者】

【氏名】汪 鋒

(72)【発明者】

【氏名】張 継成

(72)【発明者】

【氏名】馬 宏▲発▼

(72)【発明者】

【氏名】張 偉杰

(72)【発明者】

【氏名】劉 久潭

(72)【発明者】

【氏名】趙 智超

【審査官】松浦 功

(56)【参考文献】

【文献】特開２００４－３３３２４５（ＪＰ，Ａ）

【文献】中国特許出願公開第１１５６８８４７３（ＣＮ，Ａ）

【文献】中国特許出願公開第１１５６６１３８８（ＣＮ，Ａ）

【文献】中国特許出願公開第１１３５３６７０４（ＣＮ，Ａ）

【文献】中国特許出願公開第１１１４７６９００（ＣＮ，Ａ）

【文献】中国特許出願公開第１１４６９２２３４（ＣＮ，Ａ）

【文献】中国特許出願公開第１１６３８４１７７（ＣＮ，Ａ）

【文献】中国特許出願公開第１１２１３１６４２（ＣＮ，Ａ）

(58)【調査した分野】(Int.Cl.，ＤＢ名)

Ｇ０６Ｆ ３０／００ －３０／２８

Ｇ０１Ｎ ３３／２４

Ｇ０６Ｑ ５０／０２

Ｇ０１Ｖ ９／０２

(57)【特許請求の範囲】

【請求項1】

具体的に、
野外調査データ、ボーリングデータ、物理探査データ及び他の実測データにより、亀裂ネットワークの幾何学的パラメータの確率分布モデルを取得するステップＳ１と、
Ｍｏｎｔｅ－Ｃａｒｌｏ法に基づいて亀裂円板の円心Ｐｏｉｎｔ＿ｃ、亀裂傾斜角Ｄｉｐ、亀裂方位Ｏｒｉｅｎｔａｔｉｏｎ、亀裂円板の半径Ｒａｄｉｕｓ及び亀裂開口幅Ａｐｅｒｔｕｒｅのパラメータをランダムにサンプリングし、亀裂ネットワーク内の各パラメータのデータセットを生成するステップＳ２と、
ランダム形式のワイエルシュトラス関数により、粗い亀裂の上下壁面の各座標における高さ値を生成し、それに基づいて全ての亀裂の粗い亀裂間開口幅セットを取得するステップＳ３と、
ステップＳ２で取得された亀裂ネットワーク内の各パラメータのデータセットを読み取り、バウンディングボックス法により３次元ランダム亀裂ネットワークの連結性の評価を行うステップＳ４と、
オイラー角と回転行列法により、亀裂円板の座標系の変換及び粗い亀裂開口幅値の割り当てを行い、３次元離散－粗い亀裂ネットワークの数値モデルを作成するステップＳ５と、
３次元亀裂ネットワークの浸透流のシミュレーションを行い、３次元亀裂ネットワークモデルの透水係数を計算するステップＳ６と、を含み、
ステップＳ５は、具体的に、
ステップＳ２における亀裂ネットワークの幾何学的パラメータセットを読み取り、ロードし、平面グローバル座標系の原点を亀裂円板の円心とし、半径がＲ_１の亀裂円板を作成するステップＳ５．１と、
ステップＳ３で新たに生成された上下の粗い亀裂面間の亀裂開口幅の値ｂ（ｘ_ｉ，ｙ_ｊ）を読み取り、その値を亀裂円板に割り当て、亀裂円板が粗い亀裂開口幅を有するようにするステップＳ５．２と、
グローバル座標系の原点を円板の中心点Ｏに平行移動させ、即ち、グローバル座標系の原点のｘ、ｙ、ｚの３つの方向に沿った平行移動方向をそれぞれｘ_１，ｙ_１，ｚ_１とするステップＳ５．３と、
オイラー角の定理を導入し、グローバル座標系のｘ軸が亀裂の方位方向となるように回転するように座標系をｚ軸周りに（１８０－α）°回転させ、グローバル座標系のｙ軸が亀裂の方位線と重なるように座標系をｙ軸周りにβ回転させ、デカルト座標系において亀裂円板モデルを作成するステップＳ５．４と、
ステップＳ５．１～Ｓ５．４を繰り返し、逐次反復し、３次元空間内の離散－粗い亀裂ネットワークのモデリングを完成するステップＳ５．５と、を含み、
ステップＳ６は、具体的に、
ソフトプラットフォームの支配方程式とソルバーを呼び出すステップであって、３次元亀裂ネットワークの浸透流シミュレーション中の定常流がＤａｒｃｙ方程式によって支配され、具体的に、

【数27】

（式中、ｄ_ｆは亀裂開口幅、ρは亀裂密度、μは動粘度、Ｑ_ｍは亀裂流量、ｕは速度、ｐは圧力、ｋ_ｆは透水率であり、

【数28】

は微積分の勾配演算子であり、異なる方向における微分を表す）であるステップＳ６．１と、
境界条件を設定し、ｘ方向を流体の流れ方向とし、モデルの入口境界と出口境界でそれぞれ圧力境界条件として設定するステップＳ６．２と、
モデルをメッシュ化し、流れ場及び圧力場を求解するステップＳ６．３と、
亀裂ネットワークモデルの入口での浸透流速度を計算し、入口での流速を積分することで異なる亀裂開口幅分布の下での定常状態の体積流速Ｑを取得し、具体的に、
Ｑ＝Ａｕ （１４）
（式中、Ａ＝∫ｄｈであり、流体の流れ方向に垂直な断面積、ｕは断面積を通過する流体の速度である）であるステップＳ６．４と、
ダルシー方程式に従ってモデルの亀裂ネットワークモデルの透水係数Ｋを計算し、具体的に、

【数29】

（式中、Ａは流体の流れ方向に垂直な断面積、Ｌは亀裂ネットワークの流れ方向における長さ、Δｐはモデルの入口境界と出口境界の圧力差である）であるステップＳ６．５と、を含む、ことを特徴とする３次元粗い離散亀裂ネットワークの構築、評価及び浸透流シミュレーションの方法。

【請求項2】

ステップＳ１における亀裂ネットワークの幾何学的パラメータの確率分布モデルは、亀裂円板円心モデル、亀裂傾斜角モデル、亀裂方位モデル、亀裂円板半径モデル、及び亀裂開口幅モデルを含む、ことを特徴とする請求項１に記載の３次元粗い離散亀裂ネットワークの構築、評価及び浸透流シミュレーションの方法。

【請求項3】

ステップＳ２における亀裂ネットワーク内の各データセットは、具体的に、
Ｐｏｉｎｔ＿ｃ＝｛ｘ_１，ｙ_１，ｚ_１；ｘ_２，ｙ_２，ｚ_２；ｘ_３，ｙ_３，ｚ_３；・・・；ｘ_ｎ，ｙ_ｎ，ｚ_ｎ｝
Ｄｉｐ＝｛α_１；α_２；α_３；・・・；α_ｎ｝
Ｏｒｉｅｎｔａｔｉｏｎ＝｛β_１；β_２；β_３；・・・；β_ｎ｝
Ｒａｄｉｕｓ＝｛ｒ_１；ｒ_２；ｒ_３；・・・；ｒ_ｎ｝
Ａｐｅｒｔｕｒｅ＝｛ｂ_１；ｂ_２；ｂ_３；・・・；ｂ_ｎ｝
ｎ＝Ｌ_長さ・Ｌ_幅・Ｌ_高さ・ρ
（式中、ｎは亀裂モデル内の亀裂数であり、Ｌ_長さ、Ｌ_幅及びＬ_高さはそれぞれ亀裂モデルの長さ、幅、高さを表し、ｘ_ｎ，ｙ_ｎ，ｚ_ｎはｎ番目の亀裂の中心点のデカルト座標系における座標を表し、α_ｎはｎ番目の亀裂の傾斜角であり、β_ｎはｎ番目の亀裂の方位であり、ｒ_ｎはｎ番目の亀裂円板の半径であり、ｂ_ｎはｎ番目の亀裂の開口幅であり、ρは亀裂モデルの亀裂密度である）である、ことを特徴とする請求項１に記載の３次元粗い離散亀裂ネットワークの構築、評価及び浸透流シミュレーションの方法。

【請求項4】

ステップＳ３は、具体的に、
ランダム形式のワイエルシュトラス関数により、粗い亀裂の上下壁面の各座標における高さ値を生成し、上亀裂壁面の高さが

【数30】

（式中、Ｚｕｐｐｅｒ_ｉ、ｊ（ｘ_ｉ，ｙ_ｊ）は上粗い亀裂壁面内の（ｘ_ｉ，ｙ_ｊ）における平面座標の高さ、Ｃ_Ｎは標準正規分布を持つ独立乱数、Ｎは乱数の数、Ｄとλはフラクタル変数、Ａ_ＮとＢ_Ｎは［０，２π］の一様分布を持つ独立乱数、ｉ、ｊは生成される亀裂ネットワークの空間周波数解像度である）であり、Ｃ_Ｎの乱数シードを変更して差異のある粗い亀裂の上下壁面のデータセットＺｕｐｐｅｒ（ｘ_ｉ，ｙ_ｊ）とＺｌｏｗｅｒ（ｘ_ｉ，ｙ_ｊ）を生成するステップＳ３．１と、
ステップＳ２における亀裂ネットワーク内の亀裂開口幅セットＡｐｅｒｔｕｒｅ＝｛ｂ_１；ｂ_２；ｂ_３；・・・；ｂ_ｎ｝を読み取り、Ｚｕｐｐｅｒ（ｘ_ｉ，ｙ_ｊ）’＝Ｚｕｐｐｅｒ（ｘ_ｉ，ｙ_ｊ）＋ｂ_１とし、即ち、生成される亀裂のデータセットＺｕｐｐｅｒ（ｘ_ｉ，ｙ_ｊ）’に初期開口幅ｂ_１を割り当て、Ｚｕｐｐｅｒ（ｘ_ｉ，ｙ_ｊ）は初期開口幅が割り当てられない亀裂のデータセットであるステップＳ３．２と、
地殻応力作用を考慮し、応力、法線剛性を入力することで亀裂閉合の変形量を取得し、実際の粗い亀裂開口幅の変形量及び上下壁面の接触具合を取得するステップＳ３．３と、

【数31】

（式中、Δｂ_ｆは亀裂閉合の変形量、Ｖ_ｍは初期亀裂閉合の変形量、σ_ｎは垂直応力、Ｋ_ｎｉは初期法線剛性である）
上下の亀裂が変位し、部分的に閉合した後の粗い亀裂の平均亀裂開口幅を式（３）により計算するステップＳ３．４と、

【数32】

（ただし、Ｍ及びＮはそれぞれｘ方向及びｙ方向のノードの数を表す）
全ての亀裂の粗い亀裂間開口幅セットを取得するように、ステップＳ３．１～Ｓ３．４を繰り返すステップＳ３．５と、を含む、ことを特徴とする請求項１に記載の３次元粗い離散亀裂ネットワークの構築、評価及び浸透流シミュレーションの方法。

【請求項5】

ステップＳ４は、具体的に、
バウンディングボックス法により、第１亀裂と第２亀裂の空間的関係を初期的に判断し、２つの境界球の球心間の直線距離Ｌ_ｖｂ及び球の半径を比較し、２つの亀裂の境界球の
重なりを分析し、幾何学的特徴が式（４）を満たすステップＳ４．１と、

【数33】

（式中、Ｌ_ｖｂは第１亀裂と第２亀裂の境界球の間の距離、Ｏ_１（ｘ_１，ｙ_１，ｚ_１）、Ｏ_２（ｘ_２，ｙ_２，ｚ_２）はそれぞれ第１亀裂円板、第２亀裂円板の円心のグローバル座標であり、ａ_１、ａ_２はそれぞれ第１亀裂円板、第２亀裂円板の半径を表す）
Ｌ_ｖｂ＞（ａ_１＋ａ_２）の場合、

【数34】

（ｌ、ｍ、ｎは方向ベクトル、βは方位、αは傾斜角である）で亀裂円板の法線ベクトルｎ_１＝（ｌ，ｍ，ｎ）を表すステップＳ４．２と、
ｌ_１（ｘ－ｘ_１）＋ｍ_１（ｙ－ｙ_１）＋ｎ_１（ｚ－ｚ_１）＝０ （６）
（ただし、ｌ_１、ｍ_１、ｎ_１は１番目の亀裂円板の方向ベクトルである）で１番目の亀裂円板の位置する平面を表し、

【数35】

で１番目の亀裂円板の境界輪郭線を表し、
同様に、
ｌ_２（ｘ－ｘ_２）＋ｍ_２（ｙ－ｙ_２）＋ｎ_２（ｚ－ｚ_２）＝０ （８）
（ただし、ｌ_２、ｍ_２、ｎ_２は２番目の亀裂円板の方向ベクトルである）、及び

【数36】

で２番目の亀裂円板の位置する平面及び境界輪郭線をそれぞれ表すことができるステップＳ４．３と、
２つの境界球が部分的に重なると、式（１０）により２つの亀裂円板の位置する平面の夾角θを計算するステップＳ４．４と、

【数37】

夾角θ＝０の場合、２つの亀裂円板の位置する平面が平行になり、亀裂円板の交線が存在しないことを示し、そうでない場合、２つの亀裂円板が交線と交差するか否かをそれぞれ計算するステップＳ４．５と、を含む、ことを特徴とする請求項１に記載の３次元粗い離散亀裂ネットワークの構築、評価及び浸透流シミュレーションの方法。

【発明の詳細な説明】

【技術分野】

【0001】

本発明は、亀裂岩盤の浸透流計算の技術分野に関し、具体的に、３次元粗い離散亀裂ネットワークの構築、評価及び浸透流シミュレーションの方法に関する。

【背景技術】

【0002】

天然岩盤は、長期間にわたる続成作用と複雑な地質構造の変化を経て多数の節理、亀裂を形成し、岩盤の力学的性質と浸透流特性を大きく変化させ、地下工事の施工安全性に影響を与える。そのため、野外での亀裂の産状に基づいて実の３次元離散－粗い亀裂ネットワークを作成することは、地下工事の防災・減災研究の重点となった。

【0003】

近年、コンピュータ技術と数値計算方法の発展に伴い、岩石力学、岩石工学の課題への数値計算の適用が増えてきている。岩盤の構造と形態の点から、大量の野外での岩盤の亀裂の実測データの統計に基づき、確率論と数理統計理論により離散亀裂ネットワークモデリング（ＤＦＮ）を行う方法は、多くの実際の工事において適用・検証され、２次元離散亀裂ネットワークモデルの生成、可視化及び透水係数の求解などの方法とプロセスは、広く検討されている。しかし、実世界では、亀裂ネットワークは３次元亀裂として存在するので、２次元離散亀裂ネットワークモデルには誤差がある。このため、Ｐｒｉｅｓｔ、王恩志、張国強、謝静らが３次元亀裂ネットワークのモデリングを行い、亀裂パラメータ（密度、長さ、方位、傾斜角、亀裂開口幅）を生成させることで、円板モデル、スクエアモデルにより３次元亀裂ネットワークを生成させる。また、３次元亀裂ネットワークの浸透流特性に関する数値計算研究が発展している。学者たちは、Ａｎｓｙｓ Ｆｌｕｅｎｔ、ＯｐｅｎＦＯＡＭ、Ｃｏｍｓｏｌ Ｍｕｌｔｉｐｈｙｓｉｃｓなどの、ＣＦＤ法に基づく主流となった数値シミュレーションソフトにより浸透流特性の研究を行い、このようなビジネスソフトは、高度に開発されており、理論が完全であり、基本的にほとんどのシミュレーションのニーズを満たすことができる。一方では、Ｃ＋＋、Ｊａｖａ、Ｍａｔｌａｂなどのプログラミングプラットフォームにより、利用可能な３次元浸透流シミュレーションソフトを独自に開発し、研究を行い、現在、ガレルキン法、３Ｄ Ｕｎｉｆｉｅｄ ｐｉｐｅ－ｎｅｔｗｏｒｋ Ｍｅｔｈｏｄ（ＵＰＭ）及びＶＯＦ法が広く使われている。

【0004】

しかし、現在の研究には次の３つの問題がある。（１）３次元亀裂は、滑らかな平板亀裂と仮定されることが多いが、実際の亀裂は粗いものであり、亀裂の粗さは、単一のネットワーク亀裂の透水性を大きく変化させ、現在の研究では、亀裂面の粗さについて考慮できない。（２）亀裂は３次元亀裂岩盤にランダムに分布しており、分布の変動性が大きく、数が多いという特徴を持つ傾向がある。２つの亀裂シミュレーション方法を比較すると、亀裂の数が多すぎると、ビジネスソフトは、直接モデリングできないことが多く、亀裂の属性値の割り当て、境界条件の設定に多くの困難があり、所望の二次開発が必要となる。独自に開発された模擬プログラムは、計算速度が速く、収束性が高いが、開発レベルが低く、適用範囲が狭いなどの問題があることが多い。（３）亀裂ネットワークの浸透流特性のプロセスの研究には、亀裂ネットワークの生成、連結性の計算、及び透水係数の求解が含まれ、複雑な実の亀裂ネットワークの世界に対して、一度に単一の亀裂ネットワークモデルを求解するだけでは効率が悪く、大規模な自動求解の方法に欠けている。

【0005】

そこで、現在、３次元粗い亀裂ネットワークのモデリング、評価及びシミュレーションを容易にし、実際の水関連プロジェクトに寄与するために、岩層の産状に基づく３次元離散粗い亀裂岩盤のモデリング、亀裂の連結性の評価、浸透流シミュレーションの完全な技術方法が必要である。

【発明の概要】

【0006】

本発明の主な目的は、従来技術における３次元粗い亀裂ネットワークのモデリング方法、亀裂の連結性の評価及び浸透流シミュレーションのための方法に欠けているという問題を解決するために、３次元粗い離散亀裂ネットワークの構築、評価及び浸透流シミュレーションの方法を提供することにある。

【0007】

上記の目的を実現するために、本発明は、具体的に以下のステップを含む３次元粗い離散亀裂ネットワークの構築、評価及び浸透流シミュレーションの方法を提供する。

【0008】

ステップＳ１：野外調査データ、ボーリングデータ、物理探査データ及び他の実測データにより、亀裂ネットワークの幾何学的パラメータの確率分布モデルを取得する。

【0009】

ステップＳ２：プログラミングソフトであるＭａｔｌａｂにより、Ｍｏｎｔｅ－Ｃａｒｌｏ法に基づいて亀裂円板の円心Ｐｏｉｎｔ＿ｃ、亀裂傾斜角Ｄｉｐ、亀裂方位Ｏｒｉｅｎｔａｔｉｏｎ、亀裂円板の半径Ｒａｄｉｕｓ及び亀裂開口幅Ａｐｅｒｔｕｒｅのパラメータをランダムにサンプリングし、亀裂ネットワーク内の各パラメータのデータセットを生成する。

【0010】

ステップＳ３：ランダム形式のワイエルシュトラス関数により、粗い亀裂の上下壁面の各座標における高さ値を生成し、それに基づいて全ての亀裂の粗い亀裂間開口幅セットを取得する。

【0011】

ステップＳ４：ステップＳ２で取得された亀裂ネットワーク内の各パラメータのデータセットを読み取り、バウンディングボックス法により３次元ランダム亀裂ネットワークの連結性の評価を行う。

【0012】

ステップＳ５：開発プラットフォームであるＣｏｍｓｏｌ ｗｉｔｈ Ｍａｔｌａｂに基づき、オイラー角と回転行列法により、亀裂円板の座標系の変換及び粗い亀裂開口幅値の割り当てを行い、３次元離散－粗い亀裂ネットワークの数値モデルを作成する。

【0013】

ステップＳ６：３次元亀裂ネットワークの浸透流のシミュレーションを行い、３次元亀裂ネットワークモデルの透水係数を計算する。

【0014】

さらには、ステップＳ１における亀裂ネットワークの幾何学的パラメータの確率分布モデルは、亀裂円板円心モデル、亀裂傾斜角モデル、亀裂方位モデル、亀裂円板半径モデル、及び亀裂開口幅モデルを含む。

【0015】

さらには、ステップＳ２における亀裂ネットワーク内の各データセットは、具体的に次のとおりである。
Ｐｏｉｎｔ＿ｃ＝｛ｘ_１，ｙ_１，ｚ_１；ｘ_２，ｙ_２，ｚ_２；ｘ_３，ｙ_３，ｚ_３；・・・；ｘ_ｎ，ｙ_ｎ，ｚ_ｎ｝
Ｄｉｐ＝｛α_１；α_２；α_３；・・・；α_ｎ｝
Ｏｒｉｅｎｔａｔｉｏｎ＝｛β_１；β_２；β_３；・・・；β_ｎ｝
Ｒａｄｉｕｓ＝｛ｒ_１；ｒ_２；ｒ_３；・・・；ｒ_ｎ｝
Ａｐｅｒｔｕｒｅ＝｛ｂ_１；ｂ_２；ｂ_３；・・・；ｂ_ｎ｝
ｎ＝Ｌ_長さ・Ｌ_幅・Ｌ_高さ・ρ
式中、ｎは亀裂モデル内の亀裂数であり、Ｌ_長さ、Ｌ_幅及びＬ_高さはそれぞれ亀裂モデルの長さ、幅、高さを表し、ｘ_ｎ，ｙ_ｎ，ｚ_ｎはｎ番目の亀裂の中心点のデカルト座標系における座標を表し、α_ｎはｎ番目の亀裂の傾斜角であり、β_ｎはｎ番目の亀裂の方位であり、ｒ_ｎはｎ番目の亀裂円板の半径であり、ｂ_ｎはｎ番目の亀裂の開口幅であり、ρは亀裂モデルの亀裂密度である。

【0016】

さらには、ステップＳ３は、具体的に以下のステップを含む。

【0017】

ステップＳ３．１：ランダム形式のワイエルシュトラス関数により、粗い亀裂の上下壁面の各座標における高さ値を生成し、上亀裂壁面の高さが次のとおりである。

【数1】

式中、Ｚｕｐｐｅｒ_ｉ、ｊ（ｘ_ｉ，ｙ_ｊ）は上粗い亀裂壁面内の（ｘ_ｉ，ｙ_ｊ）における平面座標の高さ、Ｃ_Ｎは標準正規分布を持つ独立乱数、Ｎは乱数の数、Ｄとλはフラクタル変数、Ａ_ＮとＢ_Ｎは［０，２π］の一様分布を持つ独立乱数、ｉ、ｊは生成される亀裂ネットワークの空間周波数解像度であり、Ｃ_Ｎの乱数シードを変更して差異のある粗い亀裂の上下壁面のデータセットＺｕｐｐｅｒ（ｘ_ｉ，ｙ_ｊ）とＺｌｏｗｅｒ（ｘ_ｉ，ｙ_ｊ）を生成する。

【0018】

ステップＳ３．２：ステップＳ２における亀裂ネットワーク内の亀裂開口幅セットＡｐｅｒｔｕｒｅ＝｛ｂ_１；ｂ_２；ｂ_３；・・・；ｂ_ｎ｝を読み取り、Ｚｕｐｐｅｒ（ｘ_ｉ，ｙ_ｊ）’＝Ｚｕｐｐｅｒ（ｘ_ｉ，ｙ_ｊ）＋ｂ_１とし、即ち、生成される亀裂のデータセットＺｕｐｐｅｒ（ｘ_ｉ，ｙ_ｊ）’に初期開口幅ｂ_１を割り当て、Ｚｕｐｐｅｒ（ｘ_ｉ，ｙ_ｊ）は初期開口幅が割り当てられない亀裂のデータセットである。

【0019】

ステップＳ３．３：地殻応力作用を考慮し、応力、法線剛性を入力することで亀裂閉合の変形量を取得し、実際の粗い亀裂開口幅の変形量及び上下壁面の接触具合を取得する。

【数2】

式中、Δｂ_ｆは亀裂閉合の変形量、Ｖ_ｍは初期亀裂閉合の変形量、σ_ｎは垂直応力、Ｋ_ｎｉは初期法線剛性である。

【0020】

ステップＳ３．４：上下の亀裂が変位し、部分的に閉合した後の粗い亀裂の平均亀裂開口幅を式（３）により計算する。

【数3】

ただし、Ｍ及びＮはそれぞれｘ方向及びｙ方向のノードの数を表す。

【0021】

ステップＳ３．５：全ての亀裂の粗い亀裂間開口幅セットを取得するように、ステップＳ３．１～Ｓ３．４を繰り返す。

【0022】

さらには、ステップＳ４は、具体的に以下のステップを含む。

【0023】

ステップＳ４．１：バウンディングボックス法により、第１亀裂と第２亀裂の空間的関係を初期的に判断し、２つの境界球の球心間の直線距離Ｌ_ｖｂ及び球の半径を比較し、２つの亀裂の境界球の重なりを分析し、幾何学的特徴が式（４）を満たす。

【数4】

式中、Ｌ_ｖｂは第１亀裂と第２亀裂の境界球の間の距離、Ｏ_１（ｘ_１，ｙ_１，ｚ_１）、Ｏ_２（ｘ_２，ｙ_２，ｚ_２）はそれぞれ第１亀裂円板、第２亀裂円板の円心のグローバル座標であり、ａ_１、ａ_２はそれぞれ第１亀裂円板、第２亀裂円板の半径を表す。

【0024】

ステップＳ４．２：Ｌ_ｖｂ＞（ａ_１＋ａ_２）の場合、亀裂円板の法線ベクトルｎ_１＝（ｌ，ｍ，ｎ）を次のように表す。

【数5】

ｌ、ｍ、ｎは方向ベクトル、βは方位、αは傾斜角である。

【0025】

ステップＳ４．３：１番目の亀裂円板の位置する平面を次のように表す。
ｌ_１（ｘ－ｘ_１）＋ｍ_１（ｙ－ｙ_１）＋ｎ_１（ｚ－ｚ_１）＝０ （６）
ただし、ｌ_１、ｍ_１、ｎ_１は１番目の亀裂円板の方向ベクトルである。

【0026】

１番目の亀裂円板の境界輪郭線を次のように表す。

【数6】

【0027】

同様に、２番目の亀裂円板の位置する平面及び境界輪郭線をそれぞれ次のように表すことができる。
ｌ_２（ｘ－ｘ_２）＋ｍ_２（ｙ－ｙ_２）＋ｎ_２（ｚ－ｚ_２）＝０ （８）
（ただし、ｌ_２、ｍ_２、ｎ_２は２番目の亀裂円板の方向ベクトルである）

【数7】

【0028】

ステップＳ４．４：２つの境界球が部分的に重なると、式（１０）により２つの亀裂円板の位置する平面の夾角θを計算する。

【数8】

【0029】

ステップＳ４．５：夾角θ＝０°の場合、２つの亀裂円板の位置する平面が平行になり、亀裂円板の交線が存在しないことを示し、そうでない場合、２つの亀裂円板が交線と交差するか否かをそれぞれ計算する。

【0030】

さらには、ステップＳ５は、具体的に以下のステップを含む。

【0031】

ステップＳ５．１：ＭａｔｌａｂでステップＳ２における亀裂ネットワークの幾何学的パラメータセットを読み取り、ロードし、平面グローバル座標系の原点を亀裂円板の円心とし、半径がＲ_１の亀裂円板を作成する。

【0032】

ステップＳ５．２：ステップＳ３で新たに生成された上下の粗い亀裂面間の亀裂開口幅の値ｂ（ｘ_ｉ，ｙ_ｊ）を読み取り、開発プラットフォームであるＣｏｍｓｏｌ ｗｉｔｈ Ｍａｔｌａｂにコマンドを入力し、その値を亀裂円板に割り当て、亀裂円板が粗い亀裂開口幅を有するようにする。

【0033】

ステップＳ５．３：グローバル座標系の原点を円板の中心点Ｏに平行移動させ、即ち、グローバル座標系の原点のｘ、ｙ、ｚの３つの方向に沿った平行移動方向をそれぞれｘ_１、ｙ_１、ｚ_１とする。

【0034】

ステップＳ５．４：オイラー角の定理を導入し、グローバル座標系のｘ軸が亀裂の方位方向となるように回転するように座標系をｚ軸周りに（１８０－α）°回転させ、グローバル座標系のｙ軸が亀裂の方位線と重なるように座標系をｙ軸周りにβ回転させ、デカルト座標系において亀裂円板モデルを作成する。

【0035】

ステップＳ５．５：ステップＳ５．１～Ｓ５．４を繰り返し、逐次反復し、３次元空間内の離散－粗い亀裂ネットワークのモデリングを完成する。

【0036】

さらには、ステップＳ６は、具体的に以下のステップを含む。

【0037】

ステップＳ６．１：ソフトプラットフォームであるＣｏｍｓｏｌ ｗｉｔｈ Ｍａｔｌａｂに基づき、ソフトプラットフォームの支配方程式とソルバーを呼び出し、３次元亀裂ネットワークの浸透流シミュレーション中の定常流がＤａｒｃｙ方程式によって支配され、具体的に次のとおりである。

【0038】

【数9】

式中、ｄ_ｆは亀裂開口幅、ρは亀裂密度、μは動粘度、Ｑ_ｍは亀裂流量、ｕは速度、ｐは圧力、ｋ_ｆは透水率であり、

【数10】

は微積分の勾配演算子であり、異なる方向における微分を表す。

【0039】

ステップＳ６．２：境界条件を設定し、ｘ方向を流体の流れ方向とし、モデルの入口境界と出口境界でそれぞれ圧力境界条件として設定する。

【0040】

ステップＳ６．３：シミュレーションソフトであるＣｏｍｓｏｌにより、モデルをメッシュ化し、流れ場及び圧力場を求解する。

【0041】

ステップＳ６．４：亀裂ネットワークモデルの入口での浸透流速度を計算し、入口での流速を積分することで異なる亀裂開口幅分布の下での定常状態の体積流速Ｑを取得し、具体的に次のとおりである。
Ｑ＝Ａｕ （１４）
式中、Ａ＝∫ｄｈであり、流体の流れ方向に垂直な断面積、ｕは断面積を通過する流体の速度である。

【0042】

ステップＳ６．５：ダルシー方程式に従ってモデルの亀裂ネットワークモデルの透水係数Ｋを計算し、具体的に次のとおりである。

【数11】

式中、Ａは流体の流れ方向に垂直な断面積、Ｌは亀裂ネットワークの流れ方向における長さ、Δｐはモデルの入口境界と出口境界の圧力差である。

【0043】

本発明は、以下の有益な効果がある。

【0044】

岩石亀裂の産状に基づいて亀裂ネットワークのパラメータを生成し、オイラー角と回転行列により、ビジネスソフトプラットフォームであるＣｏｍｓｏｌにおける亀裂円板の座標系の変換の問題を解決する。

【0045】

３次元「離散－粗い」亀裂ネットワークのモデリングと可視化を実現する。

【0046】

亀裂ネットワークのモデリングから亀裂ネットワークの連結性の評価、さらに浸透流シミュレーションまでのプロセスを一括に自動で実現し、手作業のコストを削減し、計算効率を高め、大規模な数値シミュレーションの研究が可能となる。

【0047】

本発明を実施するための形態又は従来技術における技術的解決手段をより明らかに説明するために、以下、発明を実施するための形態又は従来技術の記述に使用される必要がある図面を簡単に説明するが、明らかに、下記の図面は、本発明の一部の実施形態であり、当業者であれば、創造的な労力をせずに、これらの図面に基づいて他の図面を得ることができる。

【図面の簡単な説明】

【0048】

【図1】本発明に係る３次元粗い離散亀裂ネットワークの構築、評価及び浸透流シミュレーションの方法のフローチャートを示す。

【図2】２つの境界球が分離している状態の模式図を示す。

【図3】２つの境界球が交差している状態の模式図を示す。

【図4】亀裂が互いに含まれている状態の模式図を示す。

【図5】亀裂が交差している状態の模式図を示す。

【図6】亀裂が分離している状態の模式図を示す。

【図7】ステップＳ５．１により作成された３次元離散亀裂ネットワークモデルの図を示す。

【図8】ステップＳ５．２により得られた粗い亀裂の亀裂開口幅値の割り当ての模式図を示す。

【発明を実施するための形態】

【0049】

以下、図面を参照しながら、本発明の技術的解決手段を明確に、完全に説明し、当然ながら、説明される実施例は本発明の実施例の一部であり、全ての実施例ではない。本発明における実施例に基づき、当業者が創造的な労力を要することなく得られる他の実施例は、全て本発明の保護範囲に属するものとする。

【0050】

図１に示す３次元粗い離散亀裂ネットワークの構築、評価及び浸透流シミュレーションの方法は、具体的に以下のステップを含む。

【0051】

ステップＳ１：野外調査データ、ボーリングデータ、物理探査データ及び他の実測データにより、表１に示すような亀裂ネットワークの幾何学的パラメータの確率分布モデルを取得する。

【0052】

【表1】

【0053】

ステップＳ２：プログラミングソフトであるＭａｔｌａｂにより、Ｍｏｎｔｅ－Ｃａｒｌｏ法に基づいて亀裂円板の円心Ｐｏｉｎｔ＿ｃ、亀裂傾斜角Ｄｉｐ、亀裂方位Ｏｒｉｅｎｔａｔｉｏｎ、亀裂円板の半径Ｒａｄｉｕｓ及び亀裂開口幅Ａｐｅｒｔｕｒｅのパラメータをランダムにサンプリングし、表２に示すような亀裂ネットワーク内の各パラメータのデータセットを生成する。

【0054】

【表2】

【0055】

ステップＳ３：ランダム形式のワイエルシュトラス関数により、粗い亀裂の上下壁面の各座標における高さ値を生成し、それに基づいて全ての亀裂の粗い亀裂間開口幅セットを取得する。

【0056】

ステップＳ４：ステップＳ２で取得された亀裂ネットワーク内の各パラメータのデータセットを読み取り、バウンディングボックス法により３次元ランダム亀裂ネットワークの連結性の評価を行う。

【0057】

ステップＳ５：開発プラットフォームであるＣｏｍｓｏｌ ｗｉｔｈ Ｍａｔｌａｂに基づき、オイラー角と回転行列法により、亀裂円板の座標系の変換及び粗い亀裂開口幅値の割り当てを行い、３次元離散－粗い亀裂ネットワークの数値モデルを作成する。

【0058】

ステップＳ６：３次元亀裂ネットワークの浸透流のシミュレーションを行い、３次元亀裂ネットワークモデルの透水係数を計算する。

【0059】

具体的には、実地調査研究を行い、野外調査データ、ボーリングデータ、物理探査データなどの実測データにより、亀裂ネットワークのパラメータの確率分布モデルを取得し、各確率分布モデルのパラメータ値を決定する。本発明では、Ｂａｅｃｈｅｒ亀裂円板モデルを採用し、ステップＳ１における亀裂ネットワークの幾何学的パラメータの確率分布モデルは、亀裂円板円心モデル、亀裂傾斜角モデル、亀裂方位モデル、亀裂円板半径モデル、及び亀裂開口幅モデルを含む。

【0060】

亀裂円板円心モデルは、ポアソン分布に合致し、その式は次のとおりである。

【数12】

【0061】

亀裂傾斜角モデル、亀裂方位モデルは、Ｆｉｓｈｅｒ分布に合致し、その式は次のとおりである。

【数13】

【0062】

亀裂円板半径モデルは、べき乗分布に合致し、その式は次のとおりである。
ｆ（ｘ）＝ｘ^－α－１

【0063】

亀裂開口幅モデルは、一様分布関数に合致し、その式は次のとおりである。

【数14】

【0064】

具体的には、ステップＳ２における亀裂ネットワーク内の各データセットは、具体的に次のとおりである。
Ｐｏｉｎｔ＿ｃ＝｛ｘ_１，ｙ_１，ｚ_１；ｘ_２，ｙ_２，ｚ_２；ｘ_３，ｙ_３，ｚ_３；・・・；ｘ_ｎ，ｙ_ｎ，ｚ_ｎ｝
Ｄｉｐ＝｛α_１；α_２；α_３；・・・；α_ｎ｝
Ｏｒｉｅｎｔａｔｉｏｎ＝｛β_１；β_２；β_３；・・・；β_ｎ｝
Ｒａｄｉｕｓ＝｛ｒ_１；ｒ_２；ｒ_３；・・・；ｒ_ｎ｝
Ａｐｅｒｔｕｒｅ＝｛ｂ_１；ｂ_２；ｂ_３；・・・；ｂ_ｎ｝
＝Ｌ_長さ・Ｌ_幅・Ｌ_高さ・ρ
式中、ｎは亀裂モデル内の亀裂数であり、Ｌ_長さ、Ｌ_幅及びＬ_高さはそれぞれ亀裂モデルの長さ、幅、高さを表し、ｘ_ｎ，ｙ_ｎ，ｚ_ｎはｎ番目の亀裂の中心点のデカルト座標系における座標を表し、α_ｎはｎ番目の亀裂の傾斜角であり、β_ｎはｎ番目の亀裂の方位であり、ｒ_ｎはｎ番目の亀裂円板の半径であり、ｂ_ｎはｎ番目の亀裂の開口幅であり、ρは亀裂モデルの亀裂密度である。

【0065】

具体的には、実世界では、亀裂は粗いものであり、亀裂の上下の粗い壁面間の亀裂空間（即ち、亀裂開口幅）は亀裂内の浸透流特性に影響を与える。現在、粗い亀裂の研究のほとんどは単一の亀裂に対して行われ、３次元離散亀裂ネットワーク内の粗い亀裂のモデリングと浸透流特徴に関する研究は比較的に少ない。３次元離散粗い亀裂ネットワークを生成するには、まず、亀裂ごとにそれぞれの粗い亀裂間開口幅データセットを生成する必要がある。単一の粗い亀裂の生成プロセスを例として、ステップＳ３は、具体的に以下のステップを含む。

【0066】

ステップＳ３．１：生成された亀裂開口幅をより実状に合わせるために、まず単一の亀裂の上下の粗い亀裂壁面を生成する必要がある。亀裂面の粗さは、フラクタル次元で表されることができ、ランダム形式のワイエルシュトラス関数により、粗い亀裂の上下壁面の各座標における高さ値を生成し、上亀裂壁面の高さが次のとおりである。

【0067】

【数15】

式中、Ｚｕｐｐｅｒ_ｉ、ｊ（ｘ_ｉ，ｙ_ｊ）は上粗い亀裂壁面内の（ｘ_ｉ，ｙ_ｊ）における平面座標の高さ、Ｃ_Ｎは標準正規分布を持つ独立乱数、Ｎは乱数の数、Ｄとλはフラクタル変数、Ａ_ＮとＢ_Ｎは［０，２π］の一様分布を持つ独立乱数、ｉ、ｊは生成される亀裂ネットワークの空間周波数解像度であり、通常、６４、１２８、２５６、５１２、１０２４などとされる。

【0068】

粗い亀裂の上下壁面の生成プロセスにおいて、他のパラメータを一定にしながら、Ｃ_Ｎの乱数シードを変更して差異のある粗い亀裂の上下壁面のデータセットＺｕｐｐｅｒ（ｘ_ｉ，ｙ_ｊ）とＺｌｏｗｅｒ（ｘ_ｉ，ｙ_ｊ）を生成することができる。

【0069】

ステップＳ３．２：ステップＳ２における亀裂ネットワーク内の亀裂開口幅セットＡｐｅｒｔｕｒｅ＝｛ｂ_１；ｂ_２；ｂ_３；・・・；ｂ_ｎ｝を読み取り、Ｚｕｐｐｅｒ（ｘ_ｉ，ｙ_ｊ）’＝Ｚｕｐｐｅｒ（ｘ_ｉ，ｙ_ｊ）＋ｂ_１とし、即ち、生成される亀裂のデータセットＺｕｐｐｅｒ（ｘ_ｉ，ｙ_ｊ）’に初期開口幅ｂ_１を割り当て、Ｚｕｐｐｅｒ（ｘ_ｉ，ｙ_ｊ）は初期開口幅が割り当てられない亀裂のデータセットである。

【0070】

ステップＳ３．３：地殻応力作用を考慮し、応力、法線剛性を入力することで亀裂閉合の変形量を取得し、実際の粗い亀裂開口幅の変形量及び上下壁面の接触具合を取得する。

【数16】

式中、Δｂ_ｆは亀裂閉合の変形量、Ｖ_ｍは初期亀裂閉合の変形量、σ_ｎは垂直応力、Ｋ_ｎｉは初期法線剛性である。

【0071】

ステップＳ３．４：Ｚｕｐｐｅｒ（ｘ_ｉ，ｙ_ｊ）をΔｂ_ｆだけ下げる（下に移動させる）。即ち、Ｚｕｐｐｅｒ（ｘ_ｉ，ｙ_ｊ）＝Ｚｕｐｐｅｒ（ｘ_ｉ，ｙ_ｊ）－Δｂ_ｆである。ｂ（ｘ_ｉ，ｙ_ｊ）を改めて計算し、即ち、ｂ（ｘ_ｉ，ｙ_ｊ）＝Ｚｕｐｐｅｒ（ｘ_ｉ，ｙ_ｊ）－Ｚｌｏｗｅｒ（ｘ_ｉ，ｙ_ｊ）であり、ｂ（ｘ_ｉ，ｙ_ｊ）＝Ｚｕｐｐｅｒ（ｘ_ｉ，ｙ_ｊ）－Ｚｕｐｐｅｒ（ｘ_ｉ，ｙ_ｊ）≦０の場合、このとき上下の亀裂壁面が接触し、亀裂がこの位置で閉合し、流体が通過できないことを示し、ｂ（ｘ_ｉ，ｙ_ｊ）＝０とする。データを更新し、単一の粗い亀裂間開口幅データセットを得ることができる。上下の亀裂が変位し、部分的に閉合した後の粗い亀裂の平均亀裂開口幅を式（３）により計算する。

【数17】

ただし、Ｍ及びＮはそれぞれｘ方向及びｙ方向のノードの数を表す。

【0072】

ステップＳ３．５：全ての亀裂の粗い亀裂間開口幅セットを取得するように、ステップＳ３．１～Ｓ３．４を繰り返す。

【0073】

具体的には、各亀裂ネットワークのパラメータセットを読み取り、３次元ランダム亀裂ネットワークの連結性の評価を行う。亀裂岩盤の浸透流プロセスにおいて、孤立した亀裂が流体の移動過程に実質的な影響を与えることができないため、亀裂ネットワークにおける流体の流れは、主に岩盤亀裂ネットワークの連結性に依存する。亀裂の連結性は、亀裂の交差点と線の数で表され、連結性指数ＣＩは次のように定義できる。

【数18】

式中、研究領域が３次元オブジェクトである場合、Ｑは領域内の亀裂円板間の交差する線分の数、Ｐは領域内の亀裂円板の数である。

【0074】

ステップＳ４は、具体的に以下のステップを含む。

【0075】

ステップＳ４．１：亀裂円板の交差頻度を計算するために、バウンディングボックス法により、第１亀裂と第２亀裂の空間的関係を初期的に判断し、２つの境界球の球心間の直線距離Ｌ_ｖｂ及び球の半径を比較し、２つの亀裂の境界球の重なりを分析し、図２及び図３に示すように、幾何学的特徴が式（４）を満たす。

【数19】

式中、Ｌ_ｖｂは第１亀裂と第２亀裂の境界球の間の距離、Ｏ_１（ｘ_１，ｙ_１，ｚ_１）、Ｏ_２（ｘ_２，ｙ_２，ｚ_２）はそれぞれ第１亀裂円板、第２亀裂円板の円心のグローバル座標であり、ａ_１、ａ_２はそれぞれ第１亀裂円板、第２亀裂円板の半径を表す。

【0076】

ステップＳ４．２：Ｌ_ｖｂ＞（ａ_１＋ａ_２）の場合、２つの亀裂の位置関係をさらに計算し、判断する必要がある。亀裂円板の産状（傾斜角αと方位β、傾斜角のラジアン値、方位のラジアン値に応じて変換できる）に応じて、亀裂円板の法線ベクトルｎ_１＝（ｌ，ｍ，ｎ）を次のように表す。

【数20】

ｌ、ｍ、ｎは方向ベクトル、βは方位、αは傾斜角である。

【0077】

ステップＳ４．３：１番目の亀裂円板の位置する平面を次のように表す。
ｌ_１（ｘ－ｘ_１）＋ｍ_１（ｙ－ｙ_１）＋ｎ_１（ｚ－ｚ_１）＝０ （６）
ただし、ｌ_１、ｍ_１、ｎ_１は１番目の亀裂円板の方向ベクトルである。

【0078】

１番目の亀裂円板の境界輪郭線を次のように表す。

【数21】

【0079】

同様に、２番目の亀裂円板の位置する平面及び境界輪郭線をそれぞれ次のように表すことができる。
ｌ_２（ｘ－ｘ_２）＋ｍ_２（ｙ－ｙ_２）＋ｎ_２（ｚ－ｚ_２）＝０ （８）
（ただし、ｌ_２、ｍ_２、ｎ_２は２番目の亀裂円板の方向ベクトルである）

【数22】

【0080】

ステップＳ４．４：２つの境界球が部分的に重なると、式（１０）により２つの亀裂円板の位置する平面の夾角θを計算する。

【数23】

【0081】

ステップＳ４．５：夾角θ＝０の場合、２つの亀裂円板の位置する平面が平行になり、亀裂円板の交線が存在しないことを示し、そうでない場合、式（５）～式（１０）により２つの亀裂円板が交線と交差するか否かをそれぞれ計算する。２つの亀裂円板がいずれも交線と交差する場合、２つの亀裂円板と交線との４つの空間交点により亀裂の位置関係を分析し、判断する必要がある。図４及び図５に示すように、２番目の亀裂が１番目の亀裂に含まれるか又は１番目の亀裂と交差することは、２つの亀裂円板が交差関係にあることを意味する。それと同時に、４つの座標を用いて交線の空間方程式と交差長さを求めることができる。図６には、２つの亀裂が分離しており、２つの亀裂円板が同一の交線を有するが、交差していないことを示す。ここで、Ａ_１１とＡ_１２は、それぞれ亀裂Ａと交線との１番目、２番目の交点であり、Ｂ_１１とＢ_１２は、それぞれ亀裂Ｂと交線との１番目、２番目の交点である。

【0082】

具体的には、１番目の亀裂円板の生成プロセスを例として、１番目の亀裂円板のパラメータをそれぞれＯ（ｘ_１，ｙ_１，ｚ_１）、傾斜角α_１と方位β_１、円板半径Ｒ_１として仮定する。

【0083】

ステップＳ５は、具体的に以下のステップを含む。

【0084】

ステップＳ５．１：ＭａｔｌａｂでステップＳ２における亀裂ネットワークの幾何学的パラメータセットを読み取り、ロードし、平面グローバル座標系の原点を亀裂円板の円心とし、半径がＲ_１の亀裂円板を作成し、図７に示すように、３次元離散亀裂ネットワークモデルの可視化を実現する。

【0085】

ステップＳ５．２：図８に示すように、ステップＳ３で新たに生成された上下の粗い亀裂面間の亀裂開口幅の値ｂ（ｘ_ｉ，ｙ_ｊ）を読み取り、開発プラットフォームであるＣｏｍｓｏｌ ｗｉｔｈ Ｍａｔｌａｂにコマンドを入力し、その値を亀裂円板に割り当て、亀裂円板が粗い亀裂開口幅を有するようにし、粗い亀裂の亀裂開口幅値の割り当てと可視化を実現できる。

【0086】

ステップＳ５．３：グローバル座標系の原点を円板の中心点Ｏに平行移動させ、即ち、グローバル座標系の原点のｘ、ｙ、ｚの３つの方向に沿った平行移動方向をそれぞれｘ_１、ｙ_１、ｚ_１とする。

【0087】

ステップＳ５．４：オイラー角の定理を導入し、グローバル座標系のｘ軸が亀裂の方位方向となるように回転するように座標系をｚ軸周りに（１８０－α）°回転させ、グローバル座標系のｙ軸が亀裂の方位線と重なるように座標系をｙ軸周りにβ回転させ、デカルト座標系において亀裂円板モデルを作成する。

【0088】

ステップＳ５．５：ステップＳ５．１～Ｓ５．４を繰り返し、逐次反復し、３次元空間内の離散－粗い亀裂ネットワークのモデリングを完成する。

【0089】

具体的には、ステップＳ６は、具体的に以下のステップを含む。
ステップＳ６．１：ソフトプラットフォームであるＣｏｍｓｏｌ ｗｉｔｈ Ｍａｔｌａｂに基づき、ソフトプラットフォームの支配方程式とソルバーを呼び出し、３次元亀裂ネットワークの浸透流シミュレーション中の定常流がＤａｒｃｙ方程式によって支配され、具体的に次のとおりである。

【0090】

【数24】

式中、ｄ_ｆは亀裂開口幅、ρは亀裂密度、μは動粘度、Ｑ_ｍは亀裂流量、ｕは速度、ｐは圧力、ｋ_ｆは透水率であり、

【数25】

は微積分の勾配演算子であり、異なる方向における微分を表す。

【0091】

ステップＳ６．２：境界条件を設定し、ｘ方向を流体の流れ方向とし、モデルの入口境界と出口境界でそれぞれ圧力境界条件として設定する。

【0092】

ステップＳ６．３：シミュレーションソフトであるＣｏｍｓｏｌにより、モデルをメッシュ化し、流れ場及び圧力場を求解する。

【0093】

ステップＳ６．４：亀裂ネットワークモデルの入口での浸透流速度を計算し、入口での流速を積分することで異なる亀裂開口幅分布の下での定常状態の体積流速Ｑを取得し、具体的に次のとおりである。
Ｑ＝Ａｕ （１４）
式中、Ａ＝∫ｄｈであり、流体の流れ方向に垂直な断面積、ｕは断面積を通過する流体の速度である。

【0094】

ステップＳ６．５：ダルシー方程式に従ってモデルの亀裂ネットワークモデルの透水係数Ｋを計算し、具体的に次のとおりである。

【数26】

式中、Ａは流体の流れ方向に垂直な断面積、Ｌは亀裂ネットワークの流れ方向における長さ、Δｐはモデルの入口境界と出口境界の圧力差である。

【0095】

上記は、単一の３次元粗い亀裂ネットワークのモデリング、連結性の評価及び浸透流の求解プロセスである。スクリプトプログラムを作成し、試験方案を引き続き読み取ることにより、亀裂ネットワークのパラメータセットを生成し、プログラムループを行う。

【0096】

当然ながら、上記の説明は、本発明を制限するものではなく、本発明は、上記の例に限定されず、本発明の実質的な範囲で当業者が行う変形、修正、追加又は置換も、本発明の保護範囲に属するものとすべきである。

【図1】

【図2】

【図3】

【図4】

【図5】

【図6】

【図7】

【図8】