(57)【要約】【課題】要素分割の方法に関わらず、点源モデルの解析精度を維持向上することができる数値解析装置、数値解析方法、数値解析プログラムを提供する。【解決手段】物理現象を有限要素法で数値解析するための装置であって、節点列でモデル化された点源部の周りの解析領域を有限個の要素に分割して、解析領域の物理量を有限要素法を用いて数値解析する解析手段と、点源部の周りの所定の要素の物理特性を補正する補正手段と、補正手段で補正した物理特性を用いて解析手段で数値解析して得られる物理量と、補正手段で補正する前の物理特性とに基づいて、他の物理量を算出する算出手段とを備えるようにする。【選択図】図１
2019212049.tif 000002
【特許請求の範囲】
【請求項１】
物理現象を有限要素法で数値解析するための装置であって、
節点列でモデル化された点源部の周りの解析領域を有限個の要素に分割して、解析領域の物理量を有限要素法を用いて数値解析する解析手段と、
点源部の周りの所定の要素の物理特性を補正する補正手段と、
補正手段で補正した物理特性を用いて解析手段で数値解析して得られる物理量と、補正手段で補正する前の物理特性とに基づいて、他の物理量を算出する算出手段とを備えることを特徴とする数値解析装置。
【請求項２】
補正手段は、点源部の周りの所定の要素の節点の物理量を一定値と仮定して、点源部の周りの所定の要素の物理特性を補正することを特徴とする請求項１に記載の数値解析装置。
【請求項３】
解析手段は、節点列でモデル化された注水／揚水井の周りの地下水流動を数値解析するものであり、
補正手段は、注水／揚水井の節点列の周りの節点群の全水頭値がそれらの平均的な値に等しいと仮定して、注水／揚水井の節点列に連結する要素の透水係数を補正するものであり、
算出手段は、解析結果である全水頭分布に、補正する前の透水係数を適用して流速を算出するものであることを特徴とする請求項１または２に記載の数値解析装置。
【請求項４】
物理現象を有限要素法で数値解析するための方法であって、
節点列でモデル化された点源部の周りの解析領域を有限個の要素に分割して、解析領域の物理量を有限要素法を用いて数値解析する解析ステップと、
点源部の周りの所定の要素の物理特性を補正する補正ステップと、
補正ステップで補正した物理特性を用いて解析ステップで数値解析して得られる物理量と、補正ステップで補正する前の物理特性とに基づいて、他の物理量を算出する算出ステップとを備えることを特徴とする数値解析方法。
【請求項５】
補正ステップは、点源部の周りの所定の要素の節点の物理量を一定値と仮定して、点源部の周りの所定の要素の物理特性を補正することを特徴とする請求項４に記載の数値解析方法。
【請求項６】
解析ステップは、節点列でモデル化された注水／揚水井の周りの地下水流動を数値解析するものであり、
補正ステップは、注水／揚水井の節点列の周りの節点群の全水頭値がそれらの平均的な値に等しいと仮定して、注水／揚水井の節点列に連結する要素の透水係数を補正するものであり、
算出ステップは、解析結果である全水頭分布に、補正する前の透水係数を適用して流速を算出するものであることを特徴とする請求項４または５に記載の数値解析方法。
【請求項７】
請求項４〜６のいずれか一つに記載の数値解析方法をコンピュータに実行させることを特徴とする数値解析プログラム。
【発明の詳細な説明】
【技術分野】
【０００１】
本発明は、例えば、ボーリング孔や井戸を利用した注水／揚水時の地下水流動を有限要素法の点源モデルにより数値的にシミュレーションする数値解析装置、数値解析方法、数値解析プログラムに関するものである。
【背景技術】
【０００２】
従来、大規模地下構造物（例えば岩盤備蓄、放射性廃棄物処分など）の計画・設計等における地下水流動の評価（例えば孔間透水試験の結果評価）では、有限要素法や差分法などによる数値シミュレーションが行われるが、解析モデルに小さな径のボーリング孔や井戸を実寸法でモデル化しようとすると、メッシュ分割に多大な労力を要するだけでなく、非現実的規模の解析になるため、孔径を無視して、点源の連なり（点源列）としてモデル化されることが多い。
【０００３】
注水／揚水孔を点源列ではなく、実寸法の空洞としてモデル化する場合にも、井戸周辺の流動場を精度よく解析するには、孔近傍を非常に細かい分割にする必要がある。
【０００４】
注水／揚水孔を点源としてモデル化する際には、点源周りの要素や格子を実際の孔径と等価となるような特定の大きさにする必要があるが、その場合でも解析領域に比べて非常に小さくする必要があり、解析領域全体の分割に支障をきたすため、実際の孔径と等価にすることなく分割される場合が少なくない。
【０００５】
揚水／注水孔の点源列モデルによる精度劣化の回避のため、差分法・有限体積法では多くの方法が提案されているが（例えば、非特許文献１〜８を参照）、図１４（１）〜（３）に示すように、構造格子（規則正しい格子）、あるいは、それに準ずる写像型メッシュに限られる。
【０００６】
有限要素法でも揚水／注水孔の点源列モデルによる精度劣化の回避のための研究はあるが（例えば、非特許文献９〜１３を参照）、図１４（１）〜（２）に示すように、いずれも点源周りが規則正しく合同な要素で分割されていることが条件となっている。
【０００７】
本発明者は、これまで有限要素法において、揚水／注水孔の点源列周りの要素の透水係数を補正することにより、解析精度の目覚ましい改善が可能になることを示したが（非特許文献１４、１５を参照）、図１４（１）〜（２）に示すように、点源周りの要素形状はすべて同一（合同）の場合に限っていた。
【０００８】
一方、従来の地下水の数値シミュレーション技術として、例えば特許文献１、２に示される方法が知られている。
【先行技術文献】
【非特許文献】
【０００９】
【非特許文献１】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．
【特許文献】
【００１０】
【特許文献１】特開２００２−２８６８６０号公報
【特許文献２】特開２００１−３２３４７７号公報
【発明の概要】
【発明が解決しようとする課題】
【００１１】
上述したように、従来の有限要素法や差分法による地下水浸透流解析において、注水／揚水井を点源としてモデル化する場合には、その周りの要素の大きさを実際の孔径と等価となるような特定の大きさにしない限り、精度の良い解析結果は得られないという問題がある。
【００１２】
また、点源モデルによる注水／揚水井の精度劣化を回避するための従来の方法は、図１４（１）〜（３）に示すように、構造格子やそれに準ずる写像型メッシュに限られるという問題がある。
【００１３】
このため、要素分割の方法に関わらず、点源モデルによる注水／揚水井の地下水流動解析の精度を維持向上することのできる方法が求められていた。
【００１４】
本発明は、上記に鑑みてなされたものであって、要素分割の方法に関わらず、点源モデルの解析精度を維持向上することができる数値解析装置、数値解析方法、数値解析プログラムを提供することを目的とする。
【課題を解決するための手段】
【００１５】
上記した課題を解決し、目的を達成するために、本発明に係る数値解析装置は、物理現象を有限要素法で数値解析するための装置であって、節点列でモデル化された点源部の周りの解析領域を有限個の要素に分割して、解析領域の物理量を有限要素法を用いて数値解析する解析手段と、点源部の周りの所定の要素の物理特性を補正する補正手段と、補正手段で補正した物理特性を用いて解析手段で数値解析して得られる物理量と、補正手段で補正する前の物理特性とに基づいて、他の物理量を算出する算出手段とを備えることを特徴とする。
【００１６】
また、本発明に係る他の数値解析装置は、上述した発明において、補正手段は、点源部の周りの所定の要素の節点の物理量を一定値と仮定して、点源部の周りの所定の要素の物理特性を補正することを特徴とする。
【００１７】
また、本発明に係る他の数値解析装置は、上述した発明において、解析手段は、節点列でモデル化された注水／揚水井の周りの地下水流動を数値解析するものであり、補正手段は、注水／揚水井の節点列の周りの節点群の全水頭値がそれらの平均的な値に等しいと仮定して、注水／揚水井の節点列に連結する要素の透水係数を補正するものであり、算出手段は、解析結果である全水頭分布に、補正する前の透水係数を適用して流速を算出するものであることを特徴とする。
【００１８】
また、本発明に係る数値解析方法は、物理現象を有限要素法で数値解析するための方法であって、節点列でモデル化された点源部の周りの解析領域を有限個の要素に分割して、解析領域の物理量を有限要素法を用いて数値解析する解析ステップと、点源部の周りの所定の要素の物理特性を補正する補正ステップと、補正ステップで補正した物理特性を用いて解析ステップで数値解析して得られる物理量と、補正ステップで補正する前の物理特性とに基づいて、他の物理量を算出する算出ステップとを備えることを特徴とする。
【００１９】
また、本発明に係る他の数値解析方法は、上述した発明において、補正ステップは、点源部の周りの所定の要素の節点の物理量を一定値と仮定して、点源部の周りの所定の要素の物理特性を補正することを特徴とする。
【００２０】
また、本発明に係る他の数値解析方法は、上述した発明において、解析ステップは、節点列でモデル化された注水／揚水井の周りの地下水流動を数値解析するものであり、補正ステップは、注水／揚水井の節点列の周りの節点群の全水頭値がそれらの平均的な値に等しいと仮定して、注水／揚水井の節点列に連結する要素の透水係数を補正するものであり、算出ステップは、解析結果である全水頭分布に、補正する前の透水係数を適用して流速を算出するものであることを特徴とする。
【００２１】
また、本発明に係る数値解析プログラムは、上述した数値解析方法をコンピュータに実行させることを特徴とする。
【発明の効果】
【００２２】
本発明に係る数値解析装置によれば、物理現象を有限要素法で数値解析するための装置であって、節点列でモデル化された点源部の周りの解析領域を有限個の要素に分割して、解析領域の物理量を有限要素法を用いて数値解析する解析手段と、点源部の周りの所定の要素の物理特性を補正する補正手段と、補正手段で補正した物理特性を用いて解析手段で数値解析して得られる物理量と、補正手段で補正する前の物理特性とに基づいて、他の物理量を算出する算出手段とを備えるので、要素分割の方法に関わらず、点源モデルの解析精度を維持向上することができるという効果を奏する。
【００２３】
また、本発明に係る他の数値解析装置によれば、補正手段は、点源部の周りの所定の要素の節点の物理量を一定値と仮定して、点源部の周りの所定の要素の物理特性を補正するので、物理特性を効率的に補正することができるという効果を奏する。
【００２４】
また、本発明に係る他の数値解析装置によれば、解析手段は、節点列でモデル化された注水／揚水井の周りの地下水流動を数値解析するものであり、補正手段は、注水／揚水井の節点列の周りの節点群の全水頭値がそれらの平均的な値に等しいと仮定して、注水／揚水井の節点列に連結する要素の透水係数を補正するものであり、算出手段は、解析結果である全水頭分布に、補正する前の透水係数を適用して流速を算出するものであるので、節点列でモデル化された注水／揚水井の周りの地下水の浸透流解析の精度を維持向上することができるという効果を奏する。
【００２５】
また、本発明に係る数値解析方法によれば、物理現象を有限要素法で数値解析するための方法であって、節点列でモデル化された点源部の周りの解析領域を有限個の要素に分割して、解析領域の物理量を有限要素法を用いて数値解析する解析ステップと、点源部の周りの所定の要素の物理特性を補正する補正ステップと、補正ステップで補正した物理特性を用いて解析ステップで数値解析して得られる物理量と、補正ステップで補正する前の物理特性とに基づいて、他の物理量を算出する算出ステップとを備えるので、要素分割の方法に関わらず、点源モデルの解析精度を維持向上することができるという効果を奏する。
【００２６】
また、本発明に係る他の数値解析方法によれば、補正ステップは、点源部の周りの所定の要素の節点の物理量を一定値と仮定して、点源部の周りの所定の要素の物理特性を補正するので、物理特性を効率的に補正することができるという効果を奏する。
【００２７】
また、本発明に係る他の数値解析方法によれば、解析ステップは、節点列でモデル化された注水／揚水井の周りの地下水流動を数値解析するものであり、補正ステップは、注水／揚水井の節点列の周りの節点群の全水頭値がそれらの平均的な値に等しいと仮定して、注水／揚水井の節点列に連結する要素の透水係数を補正するものであり、算出ステップは、解析結果である全水頭分布に、補正する前の透水係数を適用して流速を算出するものであるので、節点列でモデル化された注水／揚水井の周りの地下水の浸透流解析の精度を維持向上することができるという効果を奏する。
【００２８】
また、本発明に係る数値解析プログラムによれば、上述した数値解析方法をコンピュータに実行させるので、上述した数値解析方法をコンピュータを用いて効率的に実行することができるという効果を奏する。
【図面の簡単な説明】
【００２９】
【図１】図１は、本発明に係る数値解析装置、数値解析方法、数値解析プログラムの実施の形態を示す図であり、（１）は井戸節点に連結する任意形状要素群（井戸軸方向の視野）、（２）は被圧帯水層中の単一孔モデル、（３）は四面体要素（Ne=4）、（４）は井戸軸に沿った要素分割の模式図である。
【図２】図２は、本発明に係る数値解析装置、数値解析方法、数値解析プログラムの実施の形態を示す概略フローチャート図であり、（１）は本発明のアルゴリズムによる処理ルーチンを通常のＦＥＭ解析プログラムに組み入れた場合、（２）は本発明のアルゴリズムにより入力データを補正する場合である。
【図３】図３は、本発明のアルゴリズムによる処理ルーチンのフローチャート図である。
【図４】図４は、均等四面体要素で分割した検証対象モデルを示す図である。
【図５】図５は、図４の解析結果の比較図であり、（１）は比較例、（２）は実施例の場合である。
【図６】図６は、不均等四面体要素で分割した検証対象モデルを示す図である。
【図７】図７は、図６の解析結果の比較図であり、（１）は比較例、（２）は実施例の場合である。
【図８】図８は、均等六面体要素で分割した検証対象モデルを示す図である。
【図９】図９は、図８の解析結果の比較図であり、（１）は比較例、（２）は実施例の場合である。
【図１０】図１０は、不均等五面体要素で分割した検証対象モデルを示す図である。
【図１１】図１１は、図１０の解析結果の比較図であり、（１）は比較例、（２）は実施例の場合である。
【図１２】図１２は、不均等六面体要素で分割した検証対象モデルを示す図である。
【図１３】図１３は、図１２の解析結果の比較図であり、（１）は比較例、（２）は実施例の場合である。
【図１４】図１４は、従来の揚水／注水孔周辺の要素分割例であり、（１）は井戸節点に連結する合同な三角形要素群（井戸軸方向の視野）、（２）は井戸節点に連結する合同な四角形要素群（井戸軸方向の視野）、（３）は写像型メッシュ（出典：非特許文献８）である。
【発明を実施するための形態】
【００３０】
以下に、本発明に係る数値解析装置、数値解析方法、数値解析プログラムの実施の形態を図面に基づいて詳細に説明する。なお、この実施の形態によりこの発明が限定されるものではない。
【００３１】
本実施の形態は、有限要素法において節点列でモデル化された注水／揚水井（以降、井戸という。）の周りの地下水流動を高精度に解析するためのアルゴリズムを利用するものである。
【００３２】
本実施の形態によれば、図１（１）に示すように、任意のメッシュ分割にも適用可能である。また、メッシュ要素として、一般的な二次元要素（三角形、四角形）、三次元要素（四面体、五面体（三角柱）、六面体（四角柱））などに対応可能である。
【００３３】
また、点源周りの透水係数を補正するのみなので、後述するように、既存の有限要素コードに容易に組み込み可能である（図２（１）を参照）。また、後述するように、既存の有限要素コードに組み込むのではなく、入力データの透水係数を補正することも可能である（図２（２）を参照）。解析結果の全水頭分布から流速を求める際には、補正前の透水係数を使えばよい。
【００３４】
透水係数を補正するための補正係数の算出に当たり、点源至近の節点群の全水頭値は、それらの平均的な値に等しいと仮定する。ただし、実際の点源至近の水頭値は有限要素解析で求めるので、等しくはならない。
【００３５】
次に、図１（１）、（２）に示す被圧帯水層中の単一孔モデル問題の理論解を応用して、点源列に連結する有限要素の透水係数を補正する場合を例にとり説明する。
【００３６】
【数１】
2019212049.tif 000003
【００３７】
ここに、Qtheoryは井戸の流量、rwは井戸の半径、φw^は井戸の全水頭値、kは地盤の透水係数、Dは井戸の長さ、Rは井戸中心からの距離、φR^は距離Rの位置における全水頭値である。^（ハット）は既知量を表す。
【００３８】
式（１）を用いて透水係数を補正する方法は、本発明者が既に論文（非特許文献１４、１５を参照）で発表済みであるが、それは図１４（１）、（２）に示すように井戸の周りが合同な形状の要素で分割されている場合に限る。これに対し、本実施の形態は、図１（１）に示すような任意のメッシュ分割に対応するためのアルゴリズムである。任意のメッシュとは言え、通常は、解析精度の悪化を避けるため、井戸の周りを極端に大きさの異なる要素群で分割することは行われない。井戸の至近の浸透流場は、井戸の圧力・流量の影響が支配的であるので、本実施の形態では、井戸至近の節点群の水頭値は、それらの平均的な値に等しいとおいて補正方法を導出する。以下、本実施の形態の具体的な手順を示す。
【００３９】
有限要素法により得られる各要素の方程式は、次のように書ける。
【数２】
2019212049.tif 000004
【００４０】
ここに、インデックスeは要素に関する量を表し、
Ceij：要素eの係数行列で、
Ne：要素eを構成する節点の数
φej：要素eの構成節点jの水頭値
qei：要素eの構成節点iの流量値
【００４１】
透水テンソルは等方性、要素内で一定とし、要素eの透水係数をkeと表すと、
【数３】
2019212049.tif 000005
【００４２】
ここで、井戸を点源としてモデル化する際に、井戸節点（井戸を構成する節点、図１（３）を参照）に連結する要素（以降、連結要素という。）の透水係数keを補正することを考える。一つの要素が異なる２本の井戸に連結することは考えないとすると、要素の頂点にのみ節点を有する要素の場合（本発明では、例えば３節点三角形要素、４節点四角形要素、４節点四面体要素、６節点三角柱要素、８節点六面体要素などがある）、要素構成節点のうち、井戸に連結する節点は一つまたは二つに限るものとする。三つ以上もありうるが、極端にゆがんだ要素となり、実際には用いられないためである。以下、連結要素の構成節点のうち、一つの節点が井戸節点と一致する場合と二つが一致する場合に分けて記述する。
【００４３】
（ａ）連結要素の一つの構成節点が井戸節点に対応する場合
説明を簡単にするため、ある連結要素の構成節点のうち、第１節点が井戸節点に対応しており、水頭値φw^が与えられるとする。二次元平面解析の場合は、必ずこれに相当する。
【００４４】
要素のｉ＝１に関する方程式より
【数４】
2019212049.tif 000006
【００４５】
一般に、解の劣化を避けるため、有限要素はできるだけゆがまないように分割されるので、連結要素の井戸節点を除いた構成節点（以降、至近節点という。）は、井戸から概ね等距離に配置されているとみなせる。本実施の形態では、ある一つの井戸節点に連結する全要素の至近節点の水頭値には大差がないと仮定し、次のように近似する。
【００４６】
【数５】
2019212049.tif 000007
【００４７】
ここに、Nwは、ある一つの井戸節点を構成節点に持つ連結要素群の至近節点の集合、φs（バー）はNwに属する節点の平均的な水頭値である。したがって、式（４）は次のように近似できる。
【００４８】
【数６】
2019212049.tif 000008
【００４９】
ここで、次の係数行列の性質を利用する。
【数７】
2019212049.tif 000009
【００５０】
一般的な記述とするため、第１節点のインデックスの代わりに井戸節点を表すwを用いると、
【数８】
2019212049.tif 000010
【００５１】
（ｂ）連結要素の二つの構成節点が井戸節点に対応する場合
便宜上、ある要素の第１、第２構成節点が井戸節点に対応しており、水頭値φw^が規定されるとする。第１節点に関する方程式より
【数９】
2019212049.tif 000011
【００５２】
ここで、至近節点に関して上述と同様に考え、式（５）を用いる。
【数９Ａ】
2019212049.tif 000012
【００５３】
式（７）より、
【数１０】
2019212049.tif 000013
【数１１】
2019212049.tif 000014
【００５４】
ここで、一般性を持たせた記述とするため、ceww  =  (  ce11  +  ce12  )  、qew=qe1とおくと、上記一つの節点のみの場合と同様に
【数１２】
2019212049.tif 000015
【００５５】
（ｃ）節点に集まる流量の算出
ある井戸節点wを構成節点に有する有限要素（以降、至近要素という。）の集合をEwで表すと、井戸節点wに集まる流量Qwは、
【数１３】
2019212049.tif 000016
ここで、式（１）の理論解を用い、φw^←φR^、Rs（バー）←Rとおいて、式（１）のQtheoryと式（１３）のQwが等しいとおくと、
【００５６】
【数１４】
2019212049.tif 000017
【数１５】
2019212049.tif 000018
【００５７】
ここにDwは、井戸節点wが占有する井戸の長さであり、図１（４）に示すように、ある井戸節点に繋がる上下の井戸節点の距離の１／２である。端部の井戸節点の場合は、隣の井戸節点までの距離の１／２となる。
【００５８】
ここで一旦、ある井戸節点の連結要素群の透水係数は等しいとして、式（１５）左辺のkeの代わりにkew（〜）を、右辺のkの代わりにkeを用いると、次のようになる。
【００５９】
【数１６】
2019212049.tif 000019
【数１７】
2019212049.tif 000020
【数１８】
2019212049.tif 000021
【００６０】
keは、入力データとして与えられる補正前の透水係数、kew（〜）は井戸節点w至近の流動場を再現しうるように補正した連結要素の透水係数、Bwはkeの補正係数に対応する。また、式（１８）のφは、井戸節点から至近節点群までの平均的な距離に対応し、次のどちらかを用いる。
【００６１】
【数１９】
2019212049.tif 000022
【数２０】
2019212049.tif 000023
【００６２】
ここで、Rjは至近節点jと井戸節点wとの距離、Rs（バー）は井戸節点wにおける至近節点までの平均距離、ｍwはNwに属する至近節点の個数である（図１（１）を参照）。
【００６３】
式（１９）と式（２０）は、連結要素がすべて合同な場合など、Rjがすべて等しい場合には同じ値になる。
【００６４】
（ｄ）各要素の補正後の透水係数の算出
要素eが一つの井戸節点wにのみ連結している場合はｍew=1、要素が二つの井戸節点（井戸節点wとその前後のいずれか）に連携している場合はｍew=2として、次のいずれかで各要素の透水係数を求める。
【００６５】
【数２１】
2019212049.tif 000024
【数２２】
2019212049.tif 000025
【数２３】
2019212049.tif 000026
ここに、Newは、要素eが連結する井戸節点の集合である。
【００６６】
（実施例１）
次に、本発明の実施例１を説明する。
上述したように、本発明は、点源周りの透水係数を補正するのみなので、既存の有限要素コード（通常のＦＥＭ解析プログラム）に容易に組み込み可能である。また、既存の有限要素コードに組み込むのではなく、入力データの透水係数を補正することも可能である。これらの例について、図２を参照しながら説明する。
【００６７】
＜既存の有限要素コードに組み込む場合の例＞
図２（１）は、上記のアルゴリズムを既存の有限要素コード（ＦＥＭ解析プログラム）に組み込む場合の一例である。この図に示すように、まず、ステップＳ１０１、Ｓ１０２において、次の入力データを準備する。
［データセット１］  通常の有限要素解析に必要な入力データ
［データセット２］  次の仕様の点源列からなる井戸データ
・各井戸の井戸節点の列（井戸始点から終点、または、井戸終点から始点）
・井戸の全水頭（または圧力水頭）値φw^  。ただし、［データセット１］から引用することも可能
・井戸の半径  rw
・井戸がモデル化されている角度（対称条件などにより井戸が解析領域境界に位置しており、３６０°モデル化されていない場合、計算で求めることもできるが、与える方が効率的である）
【００６８】
次に、上記のアルゴリズムを追加したＦＥＭ解析プログラムを実行する。このプログラムは、通常のＦＥＭ解析プログラムに対応するステップＳ１０３〜Ｓ１０７と、上記のアルゴリズムに対応するステップＳ１０８、Ｓ１０９とにより構成される。
【００６９】
まず、上記のデータセット１、２を読み込む（ステップＳ１０３）。一方、上記のアルゴリズムで点源列の周りの要素の透水係数の補正を行っておく（ステップＳ１０８）。次に、データセット１、２と、補正した透水係数を入力値として、ＦＥＭ連立方程式を構築し（ステップＳ１０４）、解を求める（ステップＳ１０５）。一方、透水係数の補正の解除を行っておく（ステップＳ１０９）。次に、補正の解除をした透水係数（補正前の透水係数）を用いて、流速ベクトルを算出する等の後処理を行い（ステップＳ１０６）、結果を出力する（ステップＳ１０７）。このようにすれば、要素分割の方法に関わらず、解析精度の高い解析結果を得ることができる（ステップＳ１１０）。
【００７０】
＜入力データの透水係数を補正する場合の例＞
図２（２）は、入力データの透水係数を補正する場合の一例である。この図に示すように、まず、ステップＳ２０１、Ｓ２０２において、上述した入力データ（データセット１、２）を準備する。
【００７１】
次に、上記のアルゴリズムに基づいて、透水係数を補正するプログラムを実行する（ステップＳ２０３）。これにより、透水係数が補正されたＦＥＭ入力データが得られる（ステップＳ２０４）。また、点源列の周りの要素の透水係数の補正データが得られる（ステップＳ２１２）。
【００７２】
次に、通常のＦＥＭ解析プログラムを実行する。この場合、まず、上記のデータセット１、２、ＦＥＭ入力データを読み込む（ステップＳ２０５）。次に、ＦＥＭ連立方程式を構築し（ステップＳ２０６）、解を求める（ステップＳ２０７）。次に、流速ベクトルを算出する等の後処理を行い（ステップＳ２０８）、結果を出力する（ステップＳ２０９）。
【００７３】
その後、要素の透水係数の補正データを用いて、流速ベクトル等の透水係数に関わる物理量の逆補正プログラムを実行すると（ステップＳ２１１）、逆補正後の物理量の算出結果が得られる（Ｓ２１３）。このようにすれば、要素分割の方法に関わらず、解析精度の高い解析結果を得ることができる。
【００７４】
（実施例２）
次に、本発明の実施例２を説明する。
上記のアルゴリズムの具体的な手順としては、図３に示すようなコンピュータプログラム用のフローチャートを利用可能である。以下、これらの処理フローについて説明する。
【００７５】
図３に示すように、まず、上記のデータセット１、２を読み込む（ステップＳ３０１）。次に、井戸番号iwを１〜井戸本数nwまで以下のステップＳ３０２〜Ｓ３１３の処理を繰り返す。
【００７６】
ステップＳ３０２において、井戸iwを構成する各井戸節点の占有長さDwを算出する。次に、各井戸の井戸節点番号inwを１〜井戸節点数nnwまで以下のステップＳ３０３〜Ｓ３１１の処理を繰り返す。ここで、要素番号ieを１〜全要素数neまでステップＳ３０３〜Ｓ３０８の処理を繰り返す。
【００７７】
ステップＳ３０３において、要素ieは井戸節点inwの連結要素か否かを判定する。連結要素の場合はステップＳ３０４に進み（Ｙｅｓ）、連結要素ではない場合は次の要素について判定する（Ｎｏ）。次に、ステップＳ３０４において、井戸節点inwに連結している要素ieを連結要素の集合Ewとして記憶する。次に、ステップＳ３０５において、要素ieは井戸節点inwの前後の井戸節点にも連結しているか否かを判定する。この結果、要素が井戸節点inwのみに連結している場合は、上記の式（８）より流量値を計算する（ステップＳ３０５）。一方、要素が井戸節点inwの前後のいずれかに連結している場合は、上記の式（１２）より流量値を計算する（ステップＳ３０７）。計算後、要素ieの構成節点のうち、井戸節点以外の節点を井戸節点inwの至近節点の集合Nwとして記憶する（ステップＳ３０８）。
【００７８】
次に、Nwを利用して井戸節点inwとその至近節点群の距離を計算し（ステップＳ３０９）、上記の式（１９）または式（２０）により平均的な距離を計算する（ステップＳ３１０）。そして、上記の式（１７）により井戸節点inwでの補正後透水係数を計算し記憶する（ステップＳ３１１）。
【００７９】
次に、各井戸の井戸節点番号inwを１〜井戸節点数nnwまで、井戸節点inwに関するEwを利用して要素が連結する井戸節点数ｍewを数える処理を繰り返す（ステップＳ３１２）。続いて、要素番号ieを１〜全要素数neまで、ｍewがゼロでない場合、上記の式（２１）、式（２２）、式（２３）のいずれかで補正後透水係数を計算する処理を繰り返す（ステップＳ３１３）。この結果、要素の補正後透水係数データが得られる（ステップＳ３１４）。
【００８０】
（本発明の効果の検証）
［均等四面体要素で分割した場合］
図４に示すように、解析領域を多数の均等四面体要素で分割した検証モデルについて有限要素法による数値解析を実施して、その結果を本発明の数値解析方法（実施例）と従来のＦＥＭによる数値解析方法（比較例）とで比較した。検証モデルの井戸の半径rw  =0.05m、井戸中心からの距離R=60mとした。また、r=rwにおいてφw^=100m（井戸の全水頭値）とし、r=60mにおいてφR^=0m（距離Rの位置における全水頭値）とした。図４に示すように、厳密解と実施例はよく一致しており、紙面上では重なって区別できない部分もあるが、比較例は両者から離れていることがわかる。
【００８１】
図５（１）に示すように、比較例の場合には、全節点位置での解析結果を示すプロットが厳密解から乖離しており、バラツキも大きいことがわかる。井戸の流量QFEM（数値解析結果の流量）についても、Qexact（厳密解の流量）の１．３８倍の結果となっており、厳密解から乖離していることがわかる。
【００８２】
図５（２）に示すように、実施例の場合には、全節点位置での解析結果を示すプロットは厳密解とよく一致しており、バラツキも小さいことがわかる。井戸の流量QFEM（数値解析結果の流量）についても、Qexact（厳密解の流量）の０．９９倍となっており、厳密解に非常に近いことがわかる。
【００８３】
［不均等四面体要素で分割した場合］
図６に示すような不均等四面体要素分割による検証モデルについて実施例と比較例とで解析結果を比較した。境界条件は図４の検証モデルと同じである。図６に示すように、厳密解と実施例はよく一致しており、紙面上では重なって区別できない部分もあるが、比較例は両者から離れていることがわかる。なお、中心部下方の荒い連結要素で等値線が奇異な分布を示しているのは、等値線を要素内で線形補間して描いているためである。
【００８４】
図７（１）に示すように、比較例の場合には、全節点位置での解析結果を示すプロットが厳密解から乖離しており、バラツキも大きいことがわかる。井戸の流量の比QFEM（数値解析結果の流量）についても、Qexact（厳密解の流量）の１．６６倍の結果となっており、厳密解から乖離していることがわかる。
【００８５】
図７（２）に示すように、実施例の場合には、全節点位置での解析結果を示すプロットは厳密解とよく一致しており、バラツキも小さいことがわかる。井戸の流量QFEM（数値解析結果の流量）についても、Qexact（厳密解の流量）の０．９８倍となっており、厳密解に非常に近いことがわかる。
【００８６】
［均等六面体要素で分割した場合］
図８に示すような均等六面体要素分割による検証モデルについて実施例と比較例とで解析結果を比較した。境界条件は図４の検証モデルと同じである。図８に示すように、厳密解と実施例はよく一致しており、紙面上では重なって区別できない部分もあるが、比較例は両者から離れていることがわかる。
【００８７】
図９（１）に示すように、比較例の場合には、全節点位置での解析結果を示すプロットが厳密解から乖離しており、バラツキも大きいことがわかる。井戸の流量QFEM（数値解析結果の流量）についても、Qexact（厳密解の流量）の１．３２倍の結果となっており、厳密解から乖離していることがわかる。
【００８８】
図９（２）に示すように、実施例の場合には、全節点位置での解析結果を示すプロットは厳密解とよく一致しており、バラツキも小さいことがわかる。井戸の流量QFEM（数値解析結果の流量）についても、Qexact（厳密解の流量）の１．０１倍となっており、厳密解に非常に近いことがわかる。
【００８９】
［不均等五面体要素で分割した場合］
図１０に示すような不均等五面体要素分割による検証モデルについて実施例と比較例とで解析結果を比較した。境界条件は図４の検証モデルと同じである。図１０に示すように、厳密解と実施例はよく一致しており、紙面上では重なって区別できない部分もあるが、比較例は両者から離れていることがわかる。
【００９０】
図１１（１）に示すように、比較例の場合には、全節点位置での解析結果を示すプロットが厳密解から乖離しており、バラツキも大きいことがわかる。井戸の流量QFEM（数値解析結果の流量）についても、Qexact（厳密解の流量）の１．５３倍の結果となっており、厳密解から乖離していることがわかる。
【００９１】
図１１（２）に示すように、実施例の場合には、全節点位置での解析結果を示すプロットは厳密解とよく一致しており、バラツキも小さいことがわかる。井戸の流量QFEM（数値解析結果の流量）についても、Qexact（厳密解の流量）の０．９４倍となっており、厳密解に近いことがわかる。なお、この実施例の場合には、四面体分割の場合に比べるとバラツキが大きいが、これは故意に要素形状をゆがめた影響と考えられる。
【００９２】
［不均等六面体要素で分割した場合］
図１２に示すような不均等六面体要素分割による検証モデルについて実施例と比較例とで解析結果を比較した。境界条件は図４の検証モデルと同じである。図１２に示すように、厳密解と実施例はよく一致しており、紙面上では重なって区別できない部分もあるが、比較例は両者から離れていることがわかる。
【００９３】
図１３（１）に示すように、比較例の場合には、全節点位置での解析結果を示すプロットが厳密解から乖離しており、バラツキも大きいことがわかる。井戸の流量QFEM（数値解析結果の流量）についても、Qexact（厳密解の流量）の１．４４倍の結果となっており、厳密解から乖離していることがわかる。
【００９４】
図１３（２）に示すように、実施例の場合には、全節点位置での解析結果を示すプロットは厳密解とよく一致しており、バラツキも小さいことがわかる。井戸の流量QFEM（数値解析結果の流量）についても、Qexact（厳密解の流量）の０．９９倍となっており、厳密解に近いことがわかる。なお、この実施例の場合には、四面体分割の場合に比べるとバラツキが大きいが、これは故意に要素形状をゆがめた影響と考えられる。
【００９５】
したがって、本実施例によれば、有限要素法による任意分割のメッシュを用いた地下水流動解析において、注水／揚水孔を点源列でモデル化しても孔周辺の流動場を精度よく解析できる。このため、孔間透水試験に基づく地盤岩盤の透水係数の同定において、本実施例を適用すれば、信頼できる結果が得られる。
【００９６】
上記の実施の形態では、地下水の浸透流問題に適用する場合を例にとり説明したが、本発明はこれに限るものではなく、浸透流問題に限らず、同じ偏微分方程式で記述可能な熱伝導問題、拡散問題等にも転用可能である。
【００９７】
以上説明したように、本発明に係る数値解析装置によれば、物理現象を有限要素法で数値解析するための装置であって、節点列でモデル化された点源部の周りの解析領域を有限個の要素に分割して、解析領域の物理量を有限要素法を用いて数値解析する解析手段と、点源部の周りの所定の要素の物理特性を補正する補正手段と、補正手段で補正した物理特性を用いて解析手段で数値解析して得られる物理量と、補正手段で補正する前の物理特性とに基づいて、他の物理量を算出する算出手段とを備えるので、要素分割の方法に関わらず、点源モデルの解析精度を維持向上することができる。
【００９８】
また、本発明に係る他の数値解析装置によれば、補正手段は、点源部の周りの所定の要素の節点の物理量を一定値と仮定して、点源部の周りの所定の要素の物理特性を補正するので、物理特性を効率的に補正することができる。
【００９９】
また、本発明に係る他の数値解析装置によれば、解析手段は、節点列でモデル化された注水／揚水井の周りの地下水流動を数値解析するものであり、補正手段は、注水／揚水井の節点列の周りの節点群の全水頭値がそれらの平均的な値に等しいと仮定して、注水／揚水井の節点列に連結する要素の透水係数を補正するものであり、算出手段は、解析結果である全水頭分布に、補正する前の透水係数を適用して流速を算出するものであるので、節点列でモデル化された注水／揚水井の周りの地下水の浸透流解析の精度を維持向上することができる。
【０１００】
また、本発明に係る数値解析方法によれば、物理現象を有限要素法で数値解析するための方法であって、節点列でモデル化された点源部の周りの解析領域を有限個の要素に分割して、解析領域の物理量を有限要素法を用いて数値解析する解析ステップと、点源部の周りの所定の要素の物理特性を補正する補正ステップと、補正ステップで補正した物理特性を用いて解析ステップで数値解析して得られる物理量と、補正ステップで補正する前の物理特性とに基づいて、他の物理量を算出する算出ステップとを備えるので、要素分割の方法に関わらず、点源モデルの解析精度を維持向上することができる。
【０１０１】
また、本発明に係る他の数値解析方法によれば、補正ステップは、点源部の周りの所定の要素の節点の物理量を一定値と仮定して、点源部の周りの所定の要素の物理特性を補正するので、物理特性を効率的に補正することができる。
【０１０２】
また、本発明に係る他の数値解析方法によれば、解析ステップは、節点列でモデル化された注水／揚水井の周りの地下水流動を数値解析するものであり、補正ステップは、注水／揚水井の節点列の周りの節点群の全水頭値がそれらの平均的な値に等しいと仮定して、注水／揚水井の節点列に連結する要素の透水係数を補正するものであり、算出ステップは、解析結果である全水頭分布に、補正する前の透水係数を適用して流速を算出するものであるので、節点列でモデル化された注水／揚水井の周りの地下水の浸透流解析の精度を維持向上することができる。
【０１０３】
また、本発明に係る数値解析プログラムによれば、上述した数値解析方法をコンピュータに実行させるので、上述した数値解析方法をコンピュータを用いて効率的に実行することができる。
【産業上の利用可能性】
【０１０４】
以上のように、本発明に係る数値解析装置、数値解析方法、数値解析プログラムは、偏微分方程式で記述可能な熱伝導問題、拡散問題等を有限要素解析するのに有用であり、特に、要素分割の方法に関わらず、点源モデルの解析精度を維持向上するのに適している。
【図１】
2019212049.tif 000027
【図２】
2019212049.tif 000028
【図３】
2019212049.tif 000029
【図４】
2019212049.tif 000030
【図５】
2019212049.tif 000031
【図６】
2019212049.tif 000032
【図７】
2019212049.tif 000033
【図８】
2019212049.tif 000034
【図９】
2019212049.tif 000035
【図１０】
2019212049.tif 000036
【図１１】
2019212049.tif 000037
【図１２】
2019212049.tif 000038
【図１３】
2019212049.tif 000039
【図１４】
2019212049.tif 000040