(19)【発行国】日本国特許庁(JP)
(12)【公報種別】特許公報(B2)
(11)【特許番号】
(24)【登録日】2022-06-29
(45)【発行日】2022-07-07
(54)【発明の名称】数値解析装置、数値解析方法、数値解析プログラム
(51)【国際特許分類】
G06Q 10/04 20120101AFI20220630BHJP
G06Q 50/10 20120101ALI20220630BHJP
G01V 9/02 20060101ALI20220630BHJP
【ＦＩ】
G06Q10/04
G06Q50/10
G01V9/02
【請求項の数】
3
(21)【出願番号】P 2018108019
(22)【出願日】2018-06-05
(65)【公開番号】P2019212049
(43)【公開日】2019-12-12
【審査請求日】2021-05-27
(73)【特許権者】
【識別番号】000002299
【氏名又は名称】清水建設株式会社
(74)【代理人】
【識別番号】110002147
【氏名又は名称】特許業務法人酒井国際特許事務所
(72)【発明者】
【氏名】櫻井 英行
(72)【発明者】
【氏名】山田 俊子
【審査官】速水 雄太
(56)【参考文献】
【文献】特開２００５－２７３２１８（ＪＰ，Ａ）
【文献】米国特許出願公開第２０１０／００８８０７６（ＵＳ，Ａ１）
【文献】山田俊子，櫻井英行，鈴木誠，注水／揚水孔の簡易有限要素モデルのコード検証，土木学会論文集Ｃ（地圏工学），2017年，Ｖｏｌ．７３，Ｎｏ．４，第４５０－４５９頁
(58)【調査した分野】(Int.Cl.，ＤＢ名)
Ｇ０６Ｑ １０／００－９９／００
Ｇ０１Ｖ ９／０２
(57)【特許請求の範囲】
【請求項1】
注水／揚水井の周りの地下水の浸透流を有限要素法で数値解析するための装置であって、注水／揚水井の周りの領域を有限個の任意の形状の要素に分割して、注水／揚水井を、所定の前記要素を構成する構成節点である井戸節点または所定の前記要素を構成する構成節点である井戸節点の連なりからなる節点列でモデル化した解析領域における地下水の浸透流を、コンピュータを使用した有限要素法で数値解析する解析手段と、
一つまたは二つの前記井戸節点に連結し、この前記井戸節点を構成節点として含む前記要素である連結要素において、前記井戸節点を除いた構成節点である全ての至近節点の水頭値が、前記至近節点の平均的な水頭値に等しいと仮定して、前記連結要素の透水係数を補正する補正手段と、
前記補正手段で補正した透水係数を用いて前記解析手段により得られる前記解析領域の全水頭分布に対して、前記補正手段で補正する前の透水係数を適用して前記解析領域における浸透流の流速または流量を算出する算出手段とを備えることを特徴とする数値解析装置。
【請求項2】
注水／揚水井の周りの地下水の浸透流を有限要素法で数値解析するための方法であって、注水／揚水井の周りの領域を有限個の任意の形状の要素に分割して、注水／揚水井を、所定の前記要素を構成する構成節点である井戸節点または所定の前記要素を構成する構成節点である井戸節点の連なりからなる節点列でモデル化した解析領域における地下水の浸透流を、コンピュータを使用した有限要素法で数値解析する解析ステップと、
一つまたは二つの前記井戸節点に連結し、この前記井戸節点を構成節点として含む前記要素である連結要素において、前記井戸節点を除いた構成節点である全ての至近節点の水頭値が、前記至近節点の平均的な水頭値に等しいと仮定して、前記連結要素の透水係数を補正する補正ステップと、
前記補正ステップで補正した透水係数を用いて前記解析ステップにより得られる前記解析領域の全水頭分布に対して、前記補正ステップで補正する前の透水係数を適用して前記解析領域における浸透流の流速または流量を算出する算出ステップとを備えることを特徴とする数値解析方法。
【請求項3】
請求項２に記載の数値解析方法をコンピュータに実行させることを特徴とする数値解析プログラム。
【発明の詳細な説明】
【技術分野】
【0001】
本発明は、例えば、ボーリング孔や井戸を利用した注水／揚水時の地下水流動を有限要素法の点源モデルにより数値的にシミュレーションする数値解析装置、数値解析方法、数値解析プログラムに関するものである。
【背景技術】
【0002】
従来、大規模地下構造物（例えば岩盤備蓄、放射性廃棄物処分など）の計画・設計等における地下水流動の評価（例えば孔間透水試験の結果評価）では、有限要素法や差分法などによる数値シミュレーションが行われるが、解析モデルに小さな径のボーリング孔や井戸を実寸法でモデル化しようとすると、メッシュ分割に多大な労力を要するだけでなく、非現実的規模の解析になるため、孔径を無視して、点源の連なり（点源列）としてモデル化されることが多い。
【0003】
注水／揚水孔を点源列ではなく、実寸法の空洞としてモデル化する場合にも、井戸周辺の流動場を精度よく解析するには、孔近傍を非常に細かい分割にする必要がある。
【0004】
注水／揚水孔を点源としてモデル化する際には、点源周りの要素や格子を実際の孔径と等価となるような特定の大きさにする必要があるが、その場合でも解析領域に比べて非常に小さくする必要があり、解析領域全体の分割に支障をきたすため、実際の孔径と等価にすることなく分割される場合が少なくない。
【0005】
揚水／注水孔の点源列モデルによる精度劣化の回避のため、差分法・有限体積法では多くの方法が提案されているが（例えば、非特許文献１～８を参照）、図１４（１）～（３）に示すように、構造格子（規則正しい格子）、あるいは、それに準ずる写像型メッシュに限られる。
【0006】
有限要素法でも揚水／注水孔の点源列モデルによる精度劣化の回避のための研究はあるが（例えば、非特許文献９～１３を参照）、図１４（１）～（２）に示すように、いずれも点源周りが規則正しく合同な要素で分割されていることが条件となっている。
【0007】
本発明者は、これまで有限要素法において、揚水／注水孔の点源列周りの要素の透水係数を補正することにより、解析精度の目覚ましい改善が可能になることを示したが（非特許文献１４、１５を参照）、図１４（１）～（２）に示すように、点源周りの要素形状はすべて同一（合同）の場合に限っていた。
【0008】
一方、従来の地下水の数値シミュレーション技術として、例えば特許文献１、２に示される方法が知られている。
【先行技術文献】
【非特許文献】
【0009】
【文献】Peaceman, D. W.: Interpretation of Wellblock Pressures in Numerical Reservoir Simulation, SPEJ (June), Trans., AIME, 253., pp.183-194, 1978.
【文献】Peaceman, D. W.: Interpretation of Wellblock Pressures in Numerical Reservoir Simulation With Nonsquare Gridblocksand Anisotropic Permeability, SPEJ (June), pp.531-543, 1983.
【文献】Peaceman, D. W.: Interpretation of Wellblock Pressure in Numerical Reservoir Simulation Part 3-0ff Center and Multiple Wells Within a Wellblock, SPERE (May), Trans., AIME,289. ,pp.227-232, 1990.
【文献】Peaceman, D. W.: Representation of a Horizontal Well in Numerical Reservoir Simulation, SPE Adv. Technology Series, Vol.1, No.1. ,pp.7-16, 1993.
【文献】登坂博行, 小島圭二, 三木章夫, 千野剛司：地表流と地下水流を結合した3次元陸水シミュレーション手法の開発, 地下水学会誌, Vol.38, pp.253-267, 1996.
【文献】山石毅, 小林仁, 谷藤吉郎, 岡本明夫, 登坂博行, 小島圭二：地下石油備蓄基地建設に伴う水文・水理挙動の数値シミュレーション, 地下水学会誌, Vol.40, pp.167-183, 1998.
【文献】Ding, Y.: A Generalized 3D Well Model for Reservoir Simulation, SPE Journal., December, No.SPE30724, pp.437-450, 1996.
【文献】Wolfsteiner, C., Durlofsky, L. J., and Aziz, K.: Calculation of well index for nonconventional wells on arbitrary grids, Computational Geosciences., Vol.7, pp.61-82, 2003.
【文献】Kono, I. : The equivalent radius of a source in numerical models of groundwater flow, Proc. of JSCE, No.218, pp.103-107, 1973.
【文献】上村佳司, 榊利博, 田中良弘：浸透流解析における井戸のモデル化に関する一考察, 第28回土質工学研究発表会講演概要集, E-7, No.841, pp.2245-2246, 1993.
【文献】榊利博, 上村佳司, 田中良弘：平面2次元浸透流解析における井戸のモデル化に関する一考察, 第28回土質工学研究発表会講演概要集, E-7, No.842, pp.2247-2248, 1993.
【文献】上村佳司, 榊利博, 田中良弘：3次元FEM浸透流解析における部分貫入孔のモデル化に関する一考察, 土木学会第49回年次学術講演会講演概要集, III-84, pp.162-163, 1994.
【文献】Chen, Z. and Zhang, Y. : Well Flow Models for Various Numerical Methods, Int. J. Numer. Analysis and Modeling., Vol.6, No.3, pp.375-388, 2009.
【文献】山田俊子，櫻孔英行，鈴木誠：有限要素法を用いた浸透流解析における注水・揚水孔の実用的な簡易モデル，土木学会論文集C(地圏工学)，Vol.71，No.4，pp.407-417，2015．
【文献】山田俊子，櫻孔英行，鈴木誠：注水／揚水孔の簡易有限要素モデルのコード検証，土木学会論文集C(地圏工学)，Vol.73，No.4，pp.450-459，2017．
【特許文献】
【0010】
【文献】特開２００２－２８６８６０号公報
【文献】特開２００１－３２３４７７号公報
【発明の概要】
【発明が解決しようとする課題】
【0011】
上述したように、従来の有限要素法や差分法による地下水浸透流解析において、注水／揚水井を点源としてモデル化する場合には、その周りの要素の大きさを実際の孔径と等価となるような特定の大きさにしない限り、精度の良い解析結果は得られないという問題がある。
【0012】
また、点源モデルによる注水／揚水井の精度劣化を回避するための従来の方法は、図１４（１）～（３）に示すように、構造格子やそれに準ずる写像型メッシュに限られるという問題がある。
【0013】
このため、要素分割の方法に関わらず、点源モデルによる注水／揚水井の地下水流動解析の精度を維持向上することのできる方法が求められていた。
【0014】
本発明は、上記に鑑みてなされたものであって、要素分割の方法に関わらず、点源モデルの解析精度を維持向上することができる数値解析装置、数値解析方法、数値解析プログラムを提供することを目的とする。
【課題を解決するための手段】
【0015】
上記した課題を解決し、目的を達成するために、本発明に係る数値解析装置は、物理現象を有限要素法で数値解析するための装置であって、節点列でモデル化された点源部の周りの解析領域を有限個の要素に分割して、解析領域の物理量を有限要素法を用いて数値解析する解析手段と、点源部の周りの所定の要素の物理特性を補正する補正手段と、補正手段で補正した物理特性を用いて解析手段で数値解析して得られる物理量と、補正手段で補正する前の物理特性とに基づいて、他の物理量を算出する算出手段とを備えることを特徴とする。
【0016】
また、本発明に係る他の数値解析装置は、上述した発明において、補正手段は、点源部の周りの所定の要素の節点の物理量を一定値と仮定して、点源部の周りの所定の要素の物理特性を補正することを特徴とする。
【0017】
また、本発明に係る他の数値解析装置は、上述した発明において、解析手段は、節点列でモデル化された注水／揚水井の周りの地下水流動を数値解析するものであり、補正手段は、注水／揚水井の節点列の周りの節点群の全水頭値がそれらの平均的な値に等しいと仮定して、注水／揚水井の節点列に連結する要素の透水係数を補正するものであり、算出手段は、解析結果である全水頭分布に、補正する前の透水係数を適用して流速を算出するものであることを特徴とする。
【0018】
また、本発明に係る数値解析方法は、物理現象を有限要素法で数値解析するための方法であって、節点列でモデル化された点源部の周りの解析領域を有限個の要素に分割して、解析領域の物理量を有限要素法を用いて数値解析する解析ステップと、点源部の周りの所定の要素の物理特性を補正する補正ステップと、補正ステップで補正した物理特性を用いて解析ステップで数値解析して得られる物理量と、補正ステップで補正する前の物理特性とに基づいて、他の物理量を算出する算出ステップとを備えることを特徴とする。
【0019】
また、本発明に係る他の数値解析方法は、上述した発明において、補正ステップは、点源部の周りの所定の要素の節点の物理量を一定値と仮定して、点源部の周りの所定の要素の物理特性を補正することを特徴とする。
【0020】
また、本発明に係る他の数値解析方法は、上述した発明において、解析ステップは、節点列でモデル化された注水／揚水井の周りの地下水流動を数値解析するものであり、補正ステップは、注水／揚水井の節点列の周りの節点群の全水頭値がそれらの平均的な値に等しいと仮定して、注水／揚水井の節点列に連結する要素の透水係数を補正するものであり、算出ステップは、解析結果である全水頭分布に、補正する前の透水係数を適用して流速を算出するものであることを特徴とする。
【0021】
また、本発明に係る数値解析プログラムは、上述した数値解析方法をコンピュータに実行させることを特徴とする。
【発明の効果】
【0022】
本発明に係る数値解析装置によれば、物理現象を有限要素法で数値解析するための装置であって、節点列でモデル化された点源部の周りの解析領域を有限個の要素に分割して、解析領域の物理量を有限要素法を用いて数値解析する解析手段と、点源部の周りの所定の要素の物理特性を補正する補正手段と、補正手段で補正した物理特性を用いて解析手段で数値解析して得られる物理量と、補正手段で補正する前の物理特性とに基づいて、他の物理量を算出する算出手段とを備えるので、要素分割の方法に関わらず、点源モデルの解析精度を維持向上することができるという効果を奏する。
【0023】
また、本発明に係る他の数値解析装置によれば、補正手段は、点源部の周りの所定の要素の節点の物理量を一定値と仮定して、点源部の周りの所定の要素の物理特性を補正するので、物理特性を効率的に補正することができるという効果を奏する。
【0024】
また、本発明に係る他の数値解析装置によれば、解析手段は、節点列でモデル化された注水／揚水井の周りの地下水流動を数値解析するものであり、補正手段は、注水／揚水井の節点列の周りの節点群の全水頭値がそれらの平均的な値に等しいと仮定して、注水／揚水井の節点列に連結する要素の透水係数を補正するものであり、算出手段は、解析結果である全水頭分布に、補正する前の透水係数を適用して流速を算出するものであるので、節点列でモデル化された注水／揚水井の周りの地下水の浸透流解析の精度を維持向上することができるという効果を奏する。
【0025】
また、本発明に係る数値解析方法によれば、物理現象を有限要素法で数値解析するための方法であって、節点列でモデル化された点源部の周りの解析領域を有限個の要素に分割して、解析領域の物理量を有限要素法を用いて数値解析する解析ステップと、点源部の周りの所定の要素の物理特性を補正する補正ステップと、補正ステップで補正した物理特性を用いて解析ステップで数値解析して得られる物理量と、補正ステップで補正する前の物理特性とに基づいて、他の物理量を算出する算出ステップとを備えるので、要素分割の方法に関わらず、点源モデルの解析精度を維持向上することができるという効果を奏する。
【0026】
また、本発明に係る他の数値解析方法によれば、補正ステップは、点源部の周りの所定の要素の節点の物理量を一定値と仮定して、点源部の周りの所定の要素の物理特性を補正するので、物理特性を効率的に補正することができるという効果を奏する。
【0027】
また、本発明に係る他の数値解析方法によれば、解析ステップは、節点列でモデル化された注水／揚水井の周りの地下水流動を数値解析するものであり、補正ステップは、注水／揚水井の節点列の周りの節点群の全水頭値がそれらの平均的な値に等しいと仮定して、注水／揚水井の節点列に連結する要素の透水係数を補正するものであり、算出ステップは、解析結果である全水頭分布に、補正する前の透水係数を適用して流速を算出するものであるので、節点列でモデル化された注水／揚水井の周りの地下水の浸透流解析の精度を維持向上することができるという効果を奏する。
【0028】
また、本発明に係る数値解析プログラムによれば、上述した数値解析方法をコンピュータに実行させるので、上述した数値解析方法をコンピュータを用いて効率的に実行することができるという効果を奏する。
【図面の簡単な説明】
【0029】
【図1】図１は、本発明に係る数値解析装置、数値解析方法、数値解析プログラムの実施の形態を示す図であり、（１）は井戸節点に連結する任意形状要素群（井戸軸方向の視野）、（２）は被圧帯水層中の単一孔モデル、（３）は四面体要素（Ne=4）、（４）は井戸軸に沿った要素分割の模式図である。
【発明を実施するための形態】
【0030】
以下に、本発明に係る数値解析装置、数値解析方法、数値解析プログラムの実施の形態を図面に基づいて詳細に説明する。なお、この実施の形態によりこの発明が限定されるものではない。
【0031】
本実施の形態は、有限要素法において節点列でモデル化された注水／揚水井（以降、井戸という。）の周りの地下水流動を高精度に解析するためのアルゴリズムを利用するものである。
【0032】
本実施の形態によれば、図１（１）に示すように、任意のメッシュ分割にも適用可能である。また、メッシュ要素として、一般的な二次元要素（三角形、四角形）、三次元要素（四面体、五面体（三角柱）、六面体（四角柱））などに対応可能である。
【0033】
また、点源周りの透水係数を補正するのみなので、後述するように、既存の有限要素コードに容易に組み込み可能である（図２（１）を参照）。また、後述するように、既存の有限要素コードに組み込むのではなく、入力データの透水係数を補正することも可能である（図２（２）を参照）。解析結果の全水頭分布から流速を求める際には、補正前の透水係数を使えばよい。
【0034】
透水係数を補正するための補正係数の算出に当たり、点源至近の節点群の全水頭値は、それらの平均的な値に等しいと仮定する。ただし、実際の点源至近の水頭値は有限要素解析で求めるので、等しくはならない。
【0035】
次に、図１（１）、（２）に示す被圧帯水層中の単一孔モデル問題の理論解を応用して、点源列に連結する有限要素の透水係数を補正する場合を例にとり説明する。
【0036】
【数1】
【0037】
ここに、Qtheoryは井戸の流量、rwは井戸の半径、φw^は井戸の全水頭値、kは地盤の透水係数、Dは井戸の長さ、Rは井戸中心からの距離、φR^は距離Rの位置における全水頭値である。^（ハット）は既知量を表す。
【0038】
式（１）を用いて透水係数を補正する方法は、本発明者が既に論文（非特許文献１４、１５を参照）で発表済みであるが、それは図１４（１）、（２）に示すように井戸の周りが合同な形状の要素で分割されている場合に限る。これに対し、本実施の形態は、図１（１）に示すような任意のメッシュ分割に対応するためのアルゴリズムである。任意のメッシュとは言え、通常は、解析精度の悪化を避けるため、井戸の周りを極端に大きさの異なる要素群で分割することは行われない。井戸の至近の浸透流場は、井戸の圧力・流量の影響が支配的であるので、本実施の形態では、井戸至近の節点群の水頭値は、それらの平均的な値に等しいとおいて補正方法を導出する。以下、本実施の形態の具体的な手順を示す。
【0039】
有限要素法により得られる各要素の方程式は、次のように書ける。
【数2】
【0040】
ここに、インデックスeは要素に関する量を表し、
Ce ij：要素eの係数行列で、
Ne：要素eを構成する節点の数
φe j：要素eの構成節点jの水頭値
qe i：要素eの構成節点iの流量値
Ce ij：要素eの係数行列で、
Ne：要素eを構成する節点の数
φe j：要素eの構成節点jの水頭値
qe i：要素eの構成節点iの流量値
【0041】
透水テンソルは等方性、要素内で一定とし、要素eの透水係数をkeと表すと、
【数3】
【0042】
ここで、井戸を点源としてモデル化する際に、井戸節点（井戸を構成する節点、図１（３）を参照）に連結する要素（以降、連結要素という。）の透水係数keを補正することを考える。一つの要素が異なる２本の井戸に連結することは考えないとすると、要素の頂点にのみ節点を有する要素の場合（本発明では、例えば３節点三角形要素、４節点四角形要素、４節点四面体要素、６節点三角柱要素、８節点六面体要素などがある）、要素構成節点のうち、井戸に連結する節点は一つまたは二つに限るものとする。三つ以上もありうるが、極端にゆがんだ要素となり、実際には用いられないためである。以下、連結要素の構成節点のうち、一つの節点が井戸節点と一致する場合と二つが一致する場合に分けて記述する。
【0043】
（ａ）連結要素の一つの構成節点が井戸節点に対応する場合
説明を簡単にするため、ある連結要素の構成節点のうち、第１節点が井戸節点に対応しており、水頭値φw^が与えられるとする。二次元平面解析の場合は、必ずこれに相当する。
説明を簡単にするため、ある連結要素の構成節点のうち、第１節点が井戸節点に対応しており、水頭値φw^が与えられるとする。二次元平面解析の場合は、必ずこれに相当する。
【0044】
要素のｉ＝１に関する方程式より
【数4】
【0045】
一般に、解の劣化を避けるため、有限要素はできるだけゆがまないように分割されるので、連結要素の井戸節点を除いた構成節点（以降、至近節点という。）は、井戸から概ね等距離に配置されているとみなせる。本実施の形態では、ある一つの井戸節点に連結する全要素の至近節点の水頭値には大差がないと仮定し、次のように近似する。
【0046】
【数5】
【0047】
ここに、Nwは、ある一つの井戸節点を構成節点に持つ連結要素群の至近節点の集合、φs（バー）はNwに属する節点の平均的な水頭値である。したがって、式（４）は次のように近似できる。
【0048】
【数6】
【0049】
ここで、次の係数行列の性質を利用する。
【数7】
【0050】
一般的な記述とするため、第１節点のインデックスの代わりに井戸節点を表すwを用いると、
【数8】
【0051】
（ｂ）連結要素の二つの構成節点が井戸節点に対応する場合
便宜上、ある要素の第１、第２構成節点が井戸節点に対応しており、水頭値φw^が規定されるとする。第１節点に関する方程式より
便宜上、ある要素の第１、第２構成節点が井戸節点に対応しており、水頭値φw^が規定されるとする。第１節点に関する方程式より
【数9】
【0052】
ここで、至近節点に関して上述と同様に考え、式（５）を用いる。
【数9A】
【0053】
式（７）より、
【数10】
【数11】
【0054】
ここで、一般性を持たせた記述とするため、ce
ww = ( ce
11 + ce
12 ) 、qe
w=qe
1とおくと、上記一つの節点のみの場合と同様に
【数12】
【0055】
（ｃ）節点に集まる流量の算出
ある井戸節点wを構成節点に有する有限要素（以降、至近要素という。）の集合をEwで表すと、井戸節点wに集まる流量Qwは、
ここで、式（１）の理論解を用い、φw^←φR^、Rs（バー）←Rとおいて、式（１）のQtheoryと式（１３）のQwが等しいとおくと、
ある井戸節点wを構成節点に有する有限要素（以降、至近要素という。）の集合をEwで表すと、井戸節点wに集まる流量Qwは、
【数13】
【0056】
【数14】
【数15】
【0057】
ここにDwは、井戸節点wが占有する井戸の長さであり、図１（４）に示すように、ある井戸節点に繋がる上下の井戸節点の距離の１／２である。端部の井戸節点の場合は、隣の井戸節点までの距離の１／２となる。
【0058】
ここで一旦、ある井戸節点の連結要素群の透水係数は等しいとして、式（１５）左辺のkeの代わりにke
w（～）を、右辺のkの代わりにkeを用いると、次のようになる。
【0059】
【数16】
【数17】
【数18】
【0060】
keは、入力データとして与えられる補正前の透水係数、ke
w（～）は井戸節点w至近の流動場を再現しうるように補正した連結要素の透水係数、Bwはkeの補正係数に対応する。また、式（１８）のφは、井戸節点から至近節点群までの平均的な距離に対応し、次のどちらかを用いる。
【0061】
【数19】
【数20】
【0062】
ここで、Rjは至近節点jと井戸節点wとの距離、Rs（バー）は井戸節点wにおける至近節点までの平均距離、ｍwはNwに属する至近節点の個数である（図１（１）を参照）。
【0063】
式（１９）と式（２０）は、連結要素がすべて合同な場合など、Rjがすべて等しい場合には同じ値になる。
【0064】
（ｄ）各要素の補正後の透水係数の算出
要素eが一つの井戸節点wにのみ連結している場合はｍe w=1、要素が二つの井戸節点（井戸節点wとその前後のいずれか）に連携している場合はｍe w=2として、次のいずれかで各要素の透水係数を求める。
要素eが一つの井戸節点wにのみ連結している場合はｍe w=1、要素が二つの井戸節点（井戸節点wとその前後のいずれか）に連携している場合はｍe w=2として、次のいずれかで各要素の透水係数を求める。
【0065】
【数21】
【数22】
【数23】
【0066】
（実施例１）
次に、本発明の実施例１を説明する。
上述したように、本発明は、点源周りの透水係数を補正するのみなので、既存の有限要素コード（通常のＦＥＭ解析プログラム）に容易に組み込み可能である。また、既存の有限要素コードに組み込むのではなく、入力データの透水係数を補正することも可能である。これらの例について、図２を参照しながら説明する。
次に、本発明の実施例１を説明する。
上述したように、本発明は、点源周りの透水係数を補正するのみなので、既存の有限要素コード（通常のＦＥＭ解析プログラム）に容易に組み込み可能である。また、既存の有限要素コードに組み込むのではなく、入力データの透水係数を補正することも可能である。これらの例について、図２を参照しながら説明する。
【0067】
＜既存の有限要素コードに組み込む場合の例＞
図２（１）は、上記のアルゴリズムを既存の有限要素コード（ＦＥＭ解析プログラム）に組み込む場合の一例である。この図に示すように、まず、ステップＳ１０１、Ｓ１０２において、次の入力データを準備する。
［データセット１］ 通常の有限要素解析に必要な入力データ
［データセット２］ 次の仕様の点源列からなる井戸データ
・各井戸の井戸節点の列（井戸始点から終点、または、井戸終点から始点）
・井戸の全水頭（または圧力水頭）値φw^ 。ただし、［データセット１］から引用することも可能
・井戸の半径 rw
・井戸がモデル化されている角度（対称条件などにより井戸が解析領域境界に位置しており、３６０°モデル化されていない場合、計算で求めることもできるが、与える方が効率的である）
図２（１）は、上記のアルゴリズムを既存の有限要素コード（ＦＥＭ解析プログラム）に組み込む場合の一例である。この図に示すように、まず、ステップＳ１０１、Ｓ１０２において、次の入力データを準備する。
［データセット１］ 通常の有限要素解析に必要な入力データ
［データセット２］ 次の仕様の点源列からなる井戸データ
・各井戸の井戸節点の列（井戸始点から終点、または、井戸終点から始点）
・井戸の全水頭（または圧力水頭）値φw^ 。ただし、［データセット１］から引用することも可能
・井戸の半径 rw
・井戸がモデル化されている角度（対称条件などにより井戸が解析領域境界に位置しており、３６０°モデル化されていない場合、計算で求めることもできるが、与える方が効率的である）
【0068】
次に、上記のアルゴリズムを追加したＦＥＭ解析プログラムを実行する。このプログラムは、通常のＦＥＭ解析プログラムに対応するステップＳ１０３～Ｓ１０７と、上記のアルゴリズムに対応するステップＳ１０８、Ｓ１０９とにより構成される。
【0069】
まず、上記のデータセット１、２を読み込む（ステップＳ１０３）。一方、上記のアルゴリズムで点源列の周りの要素の透水係数の補正を行っておく（ステップＳ１０８）。次に、データセット１、２と、補正した透水係数を入力値として、ＦＥＭ連立方程式を構築し（ステップＳ１０４）、解を求める（ステップＳ１０５）。一方、透水係数の補正の解除を行っておく（ステップＳ１０９）。次に、補正の解除をした透水係数（補正前の透水係数）を用いて、流速ベクトルを算出する等の後処理を行い（ステップＳ１０６）、結果を出力する（ステップＳ１０７）。このようにすれば、要素分割の方法に関わらず、解析精度の高い解析結果を得ることができる（ステップＳ１１０）。
【0070】
＜入力データの透水係数を補正する場合の例＞
図２（２）は、入力データの透水係数を補正する場合の一例である。この図に示すように、まず、ステップＳ２０１、Ｓ２０２において、上述した入力データ（データセット１、２）を準備する。
図２（２）は、入力データの透水係数を補正する場合の一例である。この図に示すように、まず、ステップＳ２０１、Ｓ２０２において、上述した入力データ（データセット１、２）を準備する。
【0071】
次に、上記のアルゴリズムに基づいて、透水係数を補正するプログラムを実行する（ステップＳ２０３）。これにより、透水係数が補正されたＦＥＭ入力データが得られる（ステップＳ２０４）。また、点源列の周りの要素の透水係数の補正データが得られる（ステップＳ２１２）。
【0072】
次に、通常のＦＥＭ解析プログラムを実行する。この場合、まず、上記のデータセット１、２、ＦＥＭ入力データを読み込む（ステップＳ２０５）。次に、ＦＥＭ連立方程式を構築し（ステップＳ２０６）、解を求める（ステップＳ２０７）。次に、流速ベクトルを算出する等の後処理を行い（ステップＳ２０８）、結果を出力する（ステップＳ２０９）。
【0073】
その後、要素の透水係数の補正データを用いて、流速ベクトル等の透水係数に関わる物理量の逆補正プログラムを実行すると（ステップＳ２１１）、逆補正後の物理量の算出結果が得られる（Ｓ２１３）。このようにすれば、要素分割の方法に関わらず、解析精度の高い解析結果を得ることができる。
【0074】
（実施例２）
次に、本発明の実施例２を説明する。
上記のアルゴリズムの具体的な手順としては、図３に示すようなコンピュータプログラム用のフローチャートを利用可能である。以下、これらの処理フローについて説明する。
次に、本発明の実施例２を説明する。
上記のアルゴリズムの具体的な手順としては、図３に示すようなコンピュータプログラム用のフローチャートを利用可能である。以下、これらの処理フローについて説明する。
【0075】
図３に示すように、まず、上記のデータセット１、２を読み込む（ステップＳ３０１）。次に、井戸番号iwを１～井戸本数nwまで以下のステップＳ３０２～Ｓ３１３の処理を繰り返す。
【0076】
ステップＳ３０２において、井戸iwを構成する各井戸節点の占有長さDwを算出する。次に、各井戸の井戸節点番号inwを１～井戸節点数nnwまで以下のステップＳ３０３～Ｓ３１１の処理を繰り返す。ここで、要素番号ieを１～全要素数neまでステップＳ３０３～Ｓ３０８の処理を繰り返す。
【0077】
ステップＳ３０３において、要素ieは井戸節点inwの連結要素か否かを判定する。連結要素の場合はステップＳ３０４に進み（Ｙｅｓ）、連結要素ではない場合は次の要素について判定する（Ｎｏ）。次に、ステップＳ３０４において、井戸節点inwに連結している要素ieを連結要素の集合Ewとして記憶する。次に、ステップＳ３０５において、要素ieは井戸節点inwの前後の井戸節点にも連結しているか否かを判定する。この結果、要素が井戸節点inwのみに連結している場合は、上記の式（８）より流量値を計算する（ステップＳ３０５）。一方、要素が井戸節点inwの前後のいずれかに連結している場合は、上記の式（１２）より流量値を計算する（ステップＳ３０７）。計算後、要素ieの構成節点のうち、井戸節点以外の節点を井戸節点inwの至近節点の集合Nwとして記憶する（ステップＳ３０８）。
【0078】
次に、Nwを利用して井戸節点inwとその至近節点群の距離を計算し（ステップＳ３０９）、上記の式（１９）または式（２０）により平均的な距離を計算する（ステップＳ３１０）。そして、上記の式（１７）により井戸節点inwでの補正後透水係数を計算し記憶する（ステップＳ３１１）。
【0079】
次に、各井戸の井戸節点番号inwを１～井戸節点数nnwまで、井戸節点inwに関するEwを利用して要素が連結する井戸節点数ｍe
wを数える処理を繰り返す（ステップＳ３１２）。続いて、要素番号ieを１～全要素数neまで、ｍe
wがゼロでない場合、上記の式（２１）、式（２２）、式（２３）のいずれかで補正後透水係数を計算する処理を繰り返す（ステップＳ３１３）。この結果、要素の補正後透水係数データが得られる（ステップＳ３１４）。
【0080】
（本発明の効果の検証）
［均等四面体要素で分割した場合］
図４に示すように、解析領域を多数の均等四面体要素で分割した検証モデルについて有限要素法による数値解析を実施して、その結果を本発明の数値解析方法（実施例）と従来のＦＥＭによる数値解析方法（比較例）とで比較した。検証モデルの井戸の半径rw =0.05m、井戸中心からの距離R=60mとした。また、r=rwにおいてφw^=100m（井戸の全水頭値）とし、r=60mにおいてφR^=0m（距離Rの位置における全水頭値）とした。図４に示すように、厳密解と実施例はよく一致しており、紙面上では重なって区別できない部分もあるが、比較例は両者から離れていることがわかる。
［均等四面体要素で分割した場合］
図４に示すように、解析領域を多数の均等四面体要素で分割した検証モデルについて有限要素法による数値解析を実施して、その結果を本発明の数値解析方法（実施例）と従来のＦＥＭによる数値解析方法（比較例）とで比較した。検証モデルの井戸の半径rw =0.05m、井戸中心からの距離R=60mとした。また、r=rwにおいてφw^=100m（井戸の全水頭値）とし、r=60mにおいてφR^=0m（距離Rの位置における全水頭値）とした。図４に示すように、厳密解と実施例はよく一致しており、紙面上では重なって区別できない部分もあるが、比較例は両者から離れていることがわかる。
【0081】
図５（１）に示すように、比較例の場合には、全節点位置での解析結果を示すプロットが厳密解から乖離しており、バラツキも大きいことがわかる。井戸の流量QFEM（数値解析結果の流量）についても、Qexact（厳密解の流量）の１．３８倍の結果となっており、厳密解から乖離していることがわかる。
【0082】
図５（２）に示すように、実施例の場合には、全節点位置での解析結果を示すプロットは厳密解とよく一致しており、バラツキも小さいことがわかる。井戸の流量QFEM（数値解析結果の流量）についても、Qexact（厳密解の流量）の０．９９倍となっており、厳密解に非常に近いことがわかる。
【0083】
［不均等四面体要素で分割した場合］
図６に示すような不均等四面体要素分割による検証モデルについて実施例と比較例とで解析結果を比較した。境界条件は図４の検証モデルと同じである。図６に示すように、厳密解と実施例はよく一致しており、紙面上では重なって区別できない部分もあるが、比較例は両者から離れていることがわかる。なお、中心部下方の荒い連結要素で等値線が奇異な分布を示しているのは、等値線を要素内で線形補間して描いているためである。
図６に示すような不均等四面体要素分割による検証モデルについて実施例と比較例とで解析結果を比較した。境界条件は図４の検証モデルと同じである。図６に示すように、厳密解と実施例はよく一致しており、紙面上では重なって区別できない部分もあるが、比較例は両者から離れていることがわかる。なお、中心部下方の荒い連結要素で等値線が奇異な分布を示しているのは、等値線を要素内で線形補間して描いているためである。
【0084】
図７（１）に示すように、比較例の場合には、全節点位置での解析結果を示すプロットが厳密解から乖離しており、バラツキも大きいことがわかる。井戸の流量の比QFEM（数値解析結果の流量）についても、Qexact（厳密解の流量）の１．６６倍の結果となっており、厳密解から乖離していることがわかる。
【0085】
図７（２）に示すように、実施例の場合には、全節点位置での解析結果を示すプロットは厳密解とよく一致しており、バラツキも小さいことがわかる。井戸の流量QFEM（数値解析結果の流量）についても、Qexact（厳密解の流量）の０．９８倍となっており、厳密解に非常に近いことがわかる。
【0086】
［均等六面体要素で分割した場合］
図８に示すような均等六面体要素分割による検証モデルについて実施例と比較例とで解析結果を比較した。境界条件は図４の検証モデルと同じである。図８に示すように、厳密解と実施例はよく一致しており、紙面上では重なって区別できない部分もあるが、比較例は両者から離れていることがわかる。
図８に示すような均等六面体要素分割による検証モデルについて実施例と比較例とで解析結果を比較した。境界条件は図４の検証モデルと同じである。図８に示すように、厳密解と実施例はよく一致しており、紙面上では重なって区別できない部分もあるが、比較例は両者から離れていることがわかる。
【0087】
図９（１）に示すように、比較例の場合には、全節点位置での解析結果を示すプロットが厳密解から乖離しており、バラツキも大きいことがわかる。井戸の流量QFEM（数値解析結果の流量）についても、Qexact（厳密解の流量）の１．３２倍の結果となっており、厳密解から乖離していることがわかる。
【0088】
図９（２）に示すように、実施例の場合には、全節点位置での解析結果を示すプロットは厳密解とよく一致しており、バラツキも小さいことがわかる。井戸の流量QFEM（数値解析結果の流量）についても、Qexact（厳密解の流量）の１．０１倍となっており、厳密解に非常に近いことがわかる。
【0089】
［不均等五面体要素で分割した場合］
図１０に示すような不均等五面体要素分割による検証モデルについて実施例と比較例とで解析結果を比較した。境界条件は図４の検証モデルと同じである。図１０に示すように、厳密解と実施例はよく一致しており、紙面上では重なって区別できない部分もあるが、比較例は両者から離れていることがわかる。
図１０に示すような不均等五面体要素分割による検証モデルについて実施例と比較例とで解析結果を比較した。境界条件は図４の検証モデルと同じである。図１０に示すように、厳密解と実施例はよく一致しており、紙面上では重なって区別できない部分もあるが、比較例は両者から離れていることがわかる。
【0090】
図１１（１）に示すように、比較例の場合には、全節点位置での解析結果を示すプロットが厳密解から乖離しており、バラツキも大きいことがわかる。井戸の流量QFEM（数値解析結果の流量）についても、Qexact（厳密解の流量）の１．５３倍の結果となっており、厳密解から乖離していることがわかる。
【0091】
図１１（２）に示すように、実施例の場合には、全節点位置での解析結果を示すプロットは厳密解とよく一致しており、バラツキも小さいことがわかる。井戸の流量QFEM（数値解析結果の流量）についても、Qexact（厳密解の流量）の０．９４倍となっており、厳密解に近いことがわかる。なお、この実施例の場合には、四面体分割の場合に比べるとバラツキが大きいが、これは故意に要素形状をゆがめた影響と考えられる。
【0092】
［不均等六面体要素で分割した場合］
図１２に示すような不均等六面体要素分割による検証モデルについて実施例と比較例とで解析結果を比較した。境界条件は図４の検証モデルと同じである。図１２に示すように、厳密解と実施例はよく一致しており、紙面上では重なって区別できない部分もあるが、比較例は両者から離れていることがわかる。
図１２に示すような不均等六面体要素分割による検証モデルについて実施例と比較例とで解析結果を比較した。境界条件は図４の検証モデルと同じである。図１２に示すように、厳密解と実施例はよく一致しており、紙面上では重なって区別できない部分もあるが、比較例は両者から離れていることがわかる。
【0093】
図１３（１）に示すように、比較例の場合には、全節点位置での解析結果を示すプロットが厳密解から乖離しており、バラツキも大きいことがわかる。井戸の流量QFEM（数値解析結果の流量）についても、Qexact（厳密解の流量）の１．４４倍の結果となっており、厳密解から乖離していることがわかる。
【0094】
図１３（２）に示すように、実施例の場合には、全節点位置での解析結果を示すプロットは厳密解とよく一致しており、バラツキも小さいことがわかる。井戸の流量QFEM（数値解析結果の流量）についても、Qexact（厳密解の流量）の０．９９倍となっており、厳密解に近いことがわかる。なお、この実施例の場合には、四面体分割の場合に比べるとバラツキが大きいが、これは故意に要素形状をゆがめた影響と考えられる。
【0095】
したがって、本実施例によれば、有限要素法による任意分割のメッシュを用いた地下水流動解析において、注水／揚水孔を点源列でモデル化しても孔周辺の流動場を精度よく解析できる。このため、孔間透水試験に基づく地盤岩盤の透水係数の同定において、本実施例を適用すれば、信頼できる結果が得られる。
【0096】
上記の実施の形態では、地下水の浸透流問題に適用する場合を例にとり説明したが、本発明はこれに限るものではなく、浸透流問題に限らず、同じ偏微分方程式で記述可能な熱伝導問題、拡散問題等にも転用可能である。
【0097】
以上説明したように、本発明に係る数値解析装置によれば、物理現象を有限要素法で数値解析するための装置であって、節点列でモデル化された点源部の周りの解析領域を有限個の要素に分割して、解析領域の物理量を有限要素法を用いて数値解析する解析手段と、点源部の周りの所定の要素の物理特性を補正する補正手段と、補正手段で補正した物理特性を用いて解析手段で数値解析して得られる物理量と、補正手段で補正する前の物理特性とに基づいて、他の物理量を算出する算出手段とを備えるので、要素分割の方法に関わらず、点源モデルの解析精度を維持向上することができる。
【0098】
また、本発明に係る他の数値解析装置によれば、補正手段は、点源部の周りの所定の要素の節点の物理量を一定値と仮定して、点源部の周りの所定の要素の物理特性を補正するので、物理特性を効率的に補正することができる。
【0099】
また、本発明に係る他の数値解析装置によれば、解析手段は、節点列でモデル化された注水／揚水井の周りの地下水流動を数値解析するものであり、補正手段は、注水／揚水井の節点列の周りの節点群の全水頭値がそれらの平均的な値に等しいと仮定して、注水／揚水井の節点列に連結する要素の透水係数を補正するものであり、算出手段は、解析結果である全水頭分布に、補正する前の透水係数を適用して流速を算出するものであるので、節点列でモデル化された注水／揚水井の周りの地下水の浸透流解析の精度を維持向上することができる。
【0100】
また、本発明に係る数値解析方法によれば、物理現象を有限要素法で数値解析するための方法であって、節点列でモデル化された点源部の周りの解析領域を有限個の要素に分割して、解析領域の物理量を有限要素法を用いて数値解析する解析ステップと、点源部の周りの所定の要素の物理特性を補正する補正ステップと、補正ステップで補正した物理特性を用いて解析ステップで数値解析して得られる物理量と、補正ステップで補正する前の物理特性とに基づいて、他の物理量を算出する算出ステップとを備えるので、要素分割の方法に関わらず、点源モデルの解析精度を維持向上することができる。
【0101】
また、本発明に係る他の数値解析方法によれば、補正ステップは、点源部の周りの所定の要素の節点の物理量を一定値と仮定して、点源部の周りの所定の要素の物理特性を補正するので、物理特性を効率的に補正することができる。
【0102】
また、本発明に係る他の数値解析方法によれば、解析ステップは、節点列でモデル化された注水／揚水井の周りの地下水流動を数値解析するものであり、補正ステップは、注水／揚水井の節点列の周りの節点群の全水頭値がそれらの平均的な値に等しいと仮定して、注水／揚水井の節点列に連結する要素の透水係数を補正するものであり、算出ステップは、解析結果である全水頭分布に、補正する前の透水係数を適用して流速を算出するものであるので、節点列でモデル化された注水／揚水井の周りの地下水の浸透流解析の精度を維持向上することができる。
【0103】
また、本発明に係る数値解析プログラムによれば、上述した数値解析方法をコンピュータに実行させるので、上述した数値解析方法をコンピュータを用いて効率的に実行することができる。
【産業上の利用可能性】
【0104】
以上のように、本発明に係る数値解析装置、数値解析方法、数値解析プログラムは、偏微分方程式で記述可能な熱伝導問題、拡散問題等を有限要素解析するのに有用であり、特に、要素分割の方法に関わらず、点源モデルの解析精度を維持向上するのに適している。
【図1】
【図2】
【図3】
【図4】
【図5】
【図6】
【図7】
【図8】
【図9】
【図10】
【図11】
【図12】
【図13】
【図14】
