(57)【特許請求の範囲】
【請求項１】
生成手段と解析手段とを備えたコンピュータが、
前記生成手段が、解析対象範囲に分割形成されたメッシュごとに、当該メッシュの解析対象となる有限要素方程式内の整合質量行列に対して、気泡関数要素の基底が直交する条件を満たす気泡関数を求め、前記整合質量行列に、前記気泡関数の積分値を代入することによって、前記整合質量行列を対角化した対角質量行列を生成する生成工程と、
前記解析手段が、前記解析対象となる物理量の物性を表す既知物性値および前記解析対象となる物理量の解析条件を表す境界解析物理量からなる既知解析物理量を、前記生成工程によって生成されたメッシュごとの対角質量行列に代入することによって、前記解析対象となる物理量の変化を解析する解析工程と、
を実行することを特徴とする直交基底気泡関数要素数値解析方法。
【請求項２】
前記解析手段が、前記既知物性値、前記境界解析物理量および前記解析対象となる物理量の解析開始時の値を表す初期解析物理量からなる既知解析物理量を、前記生成工程によって生成されたメッシュごとに対角質量行列に代入することによって、前記解析対象となる物理量の微少時間における変化を解析することを特徴とする請求項１に記載の直交基底気泡関数要素数値解析方法。
【請求項３】
前記生成工程は、
前記整合質量行列に対して、気泡関数要素の基底が直交する下記式（１）を満たす気泡関数を求め、前記整合質量行列に、前記気泡関数の積分値を代入することによって、前記整合質量行列を対角化した対角質量行列を生成することを特徴とする請求項１または２に記載の直交基底気泡関数要素数値解析方法。
【数１】
0004729767.tif 000114
【請求項４】
前記コンピュータは、さらに、第１の算出手段と第２の算出手段とを備え、
前記第１の算出手段が、前記生成工程によって生成された前記解析対象範囲のメッシュごとの対角質量行列の総和を求めることによって、前記解析対象範囲全域の対角質量行列を算出する第１の算出工程と、
前記第２の算出手段が、前記第１の算出工程によって算出された前記解析対象範囲全域の対角質量行列の逆行列を算出する第２の算出工程と、を実行し、
前記解析工程は、
前記解析対象範囲全域についての前記既知解析物理量と前記第２の算出工程によって算出された前記対角質量行列の逆行列とを乗算することによって、前記解析対象となる物理量の変化を解析することを特徴とする請求項１または３に記載の直交基底気泡関数要素数値解析方法。
【請求項５】
請求項１〜４のいずれか一つに記載の直交基底気泡関数要素数値解析方法を前記コンピュータに実行させることを特徴とする直交基底気泡関数要素数値解析プログラム。
【請求項６】
解析対象範囲に分割形成されたメッシュごとに、当該メッシュの解析対象となる有限要素方程式内の整合質量行列に対して、気泡関数要素の基底が直交する条件を満たす気泡関数を求め、前記整合質量行列に、前記気泡関数の積分値を代入することによって、前記整合質量行列を対角化した対角質量行列を生成する生成手段と、
前記解析対象となる物理量の物性を表す既知物性値および前記解析対象となる物理量の解析条件を表す境界解析物理量からなる既知解析物理量を、前記生成工程によって生成されたメッシュごとの対角質量行列に代入することによって、前記解析対象となる物理量の変化を解析する解析手段と、
を備えることを特徴とする直交基底気泡関数要素数値解析装置。
【請求項７】
前記解析手段は、
前記既知物性値、前記境界解析物理量および前記解析対象となる物理量の解析開始時の値を表す初期解析物理量からなる既知解析物理量を、前記生成工程によって生成されたメッシュごとに対角質量行列に代入することによって、前記解析対象となる物理量の微少時間における変化を解析することを特徴とする請求項６に記載の直交基底気泡関数要素数値解析装置。
【請求項８】
前記生成手段は、
前記整合質量行列に対して、気泡関数要素の基底が直交する下記式（２）を満たす気泡関数を求め、前記整合質量行列に、前記気泡関数の積分値を代入することによって、前記整合質量行列を対角化した対角質量行列を生成することを特徴とする請求項６または７に記載の直交基底気泡関数要素数値解析装置。
【数３】
0004729767.tif 000115
【請求項９】
前記生成手段によって生成された前記解析対象範囲のメッシュごとの対角質量行列の総和を求めることによって、前記解析対象範囲全域の対角質量行列を算出する第１の算出手段と、
前記第１の算出手段によって算出された前記解析対象範囲全域の対角質量行列の逆行列を算出する第２の算出手段と、を備え、
前記解析手段は、
前記解析対象範囲全域についての前記既知解析物理量と前記第２の算出手段によって算出された前記対角質量行列の逆行列とを乗算することによって、前記解析対象となる物理量の変化を解析することを特徴とする請求項６または８に記載の直交基底気泡関数要素数値解析装置。
【発明の詳細な説明】
【技術分野】
【０００１】
この発明は、気泡関数要素を用いた有限要素法による解析（有限要素解析）について、計算効率の良い対角項のみとなる質量行列を用いて、信頼性の高い数値シミュレーションをおこなうための直交基底気泡関数要素数値解析方法、直交基底気泡関数要素数値解析プログラム、および直交基底気泡関数要素数値解析装置に関する。
【背景技術】
【０００２】
従来の気泡関数要素について説明する。図４７は、従来の２次元の気泡関数要素を示す説明図であり、図４８は、従来の３次元の気泡関数要素を示す説明図である。図４７および図４８のように、三角形（四面体）要素を用いた気泡関数要素は、各要素において三角形（四面体）を形成する３（４）点と重心点の４（５）つの節点を用いて、アイソパラメトリック座標系[ｒ,ｓ]（｛ｒ,ｓ,ｔ｝）で下記式（１）のように表される（たとえば、下記非特許文献１、非特許文献２、非特許文献３参照）。【０００３】
【数１】
0004729767.tif 000002
【０００４】
式（１）のΦα,φBは、気泡関数要素の形状関数、ｕα,ｕBは三角形（四面体）の各頂点の値（解析物理量）、重心点の値（解析物理量）、Ｎは空間次元数を示している。ベクトル形式で記述すれば、形状関数は下記式（３）〜式（６）になる。【０００５】
【数２】
0004729767.tif 000003
【０００６】
式（２）のΨαは、２次元および３次元の一次要素を用いた形状関数であり、下記式（７）、（８）であらわされる。【０００７】
【数３】
0004729767.tif 000004
【０００８】
形状関数φBは気泡関数と呼ばれている。気泡関数は要素境界上においてその値が０となり、重心点で値が１となるように要素毎に定義される。非定常問題において、空間方向の離散化に気泡関数要素を用いた有限要素方程式は下記式（９）のようにあらわすことができる。【０００９】
【数４】
0004729767.tif 000005
【００１０】
式（９）のｕは求めるべき未知解析物理量（汚染物質濃度、温度、流量、水深、流速、圧力、変位など）であり、Ｍは質量行列、Ｆ（ｕ）は時間微分項以外をまとめた項である。式（９）の時間方向の離散化として、テイラー展開に基づいた４段解法は下記式（１０）〜式（１３）のように表される（たとえば、下記非特許文献４参照）。【００１１】
【数５】
0004729767.tif 000006
【００１２】
式（１０）〜式（１３）の上付き添え字ｎは現在時刻ｎでの既知解析物理量を表し、ｎ＋１は時刻ｎから微小時間Δｔ経過後の未知解析物理量を表している。【００１３】
【非特許文献１】D.N.Arnold, F.Brezzi and M.Fortin, “A Stable Finite Element for the Stokes Equations”, Calcolo, Vol.23, 1984, pp.337-pp.344 （ディー・エヌ・アーノルド、エフ・ブリジィ、エム・フォーティン著  「ア  ステイブル  ファイナイト  エレメント  フォー  ザ  ストークス  イクエイションズ」  カルコロ  ２３巻  １９８４年  ３３７頁−３４４頁）
【非特許文献２】J.C.Simo, F.Armero and C.A.Taylor, “Stable and Time-Dissipative Finite Element Methods for theIncompressible Navier-Stokes Equations in Advection Dominated Flows”, International Journal for Numerical Methods in Engineering,Vol.38, 1995, pp.1475-pp.1506 （ジェー・シー・シモ、エフ・アルメロ、シー・エー・テイラー著  「ステイブル  アンド  タイム−ディシペイティブ  ファイナイト  エレメント  メソーズ  フォー  ザ  インコンプレッシブル  ナビエ−ストークス  イクエイションズ  イン  アドベクション  ドミネーティド  フロウズ」  インターナショナル  ジャーナル  フォー  ニューメリカル  メソーズ  イン  エンジニアリング  ３８巻  １９９５年  １４７５頁−１５０６頁）
【非特許文献３】松本純一,「気泡関数を用いた非圧縮性粘性流れ解析のための2レベル-3レベル有限要素法」,応用力学論文集（土木学会）,７巻,２００４年８月,３３９頁−３４６頁
【非特許文献４】畑中勝守,「多段階有限要素法による非圧縮粘性流体の順・逆解析に関する計算力学的研究」,中央大学博士論文,１９９３年３月
【発明の開示】
【発明が解決しようとする課題】
【００１４】
ここで、上述した式（１０）〜式（１３）に示したように、４段解法を使用して、この既知解析物理量ｕnから未知解析物理量ｕn+1を求めるためには、質量行列の逆行列が必要になる。図４９は、従来のRotating  Cone問題の解析に用いる解析モデルを示す説明図であり、図５０は、従来のRotating  Cone問題の解析に用いる初期条件の等高線を示す説明図であり、図５１は、従来の整合質量行列によって計算されたRotating  Cone問題の解析結果の等高線を示す説明図である（非特許文献４参照）。【００１５】
図４９の解析モデル４９００を初期状態（０周時）とする。そして、質量行列の逆行列を用いて、初期状態（０周時）から、所定周、たとえば、５周回転させると、図５０に示した解析モデル４９００の初期状態における等高線モデル５０００は、図５１に示したような解析結果（等高線モデル５１００）となる。【００１６】
上述した質量行列は、時間方向の離散化に気泡関数要素を用いているため、一般的に疎な分布行列（整合質量行列）となる。したがって、この分布行列（整合質量行列）の逆行列を求めることは、数値解析上多くの記憶容量および計算時間が必要となり、装置本体のコストがかかるとともに、解析処理が遅延するという問題があった。【００１７】
この問題を解決するために、通常は、質量行列の各行の成分を足し合わせて（集中化させて）、対角項のみに成分をもたせた近似行列（集中質量行列）が使用される。集中質量行列を用いた場合には、行列の成分が対角項のみであるので、逆行列は各対角成分を逆数にした行列になり、数値解析上、近似なしの整合質量行列およびその逆行列を使用する場合に比べて、非常に少ない記憶容量、計算時間で解析を実行することができる。【００１８】
しかしながら、上述した集中質量行列を使用した場合には、集中質量行列が元の質量行列と同一の行列ではなく近似行列となる。したがって、図４９に示した解析モデル４９００の初期状態（０周時）から、所定周、たとえば、５周回転させると、図５０に示した、解析モデル４９００の初期状態における等高線モデル５０００は、図５２に示したような解析結果（等高線モデル５２００）となる。このように、解析モデル４９００の５周回転後の等高線モデル５２００は、初期状態における等高線モデル５０００や図５１に示した等高線モデル５１００に対して、大幅に変形しているため計算精度が悪く、解析結果の信頼性が低いという問題があった。【００１９】
この発明は、上述した従来技術による問題点を解消するため、簡単かつ信頼性の高い有限要素解析を実現することができる直交基底気泡関数要素数値解析方法、直交基底気泡関数要素数値解析プログラム、および直交基底気泡関数要素数値解析装置を提供することを目的とする。
【課題を解決するための手段】
【００２０】
上述した問題点を解決し、目的を達成するため、この発明の直交基底気泡関数要素数値解析方法、直交基底気泡関数要素数値解析プログラム、および直交基底気泡関数要素数値解析装置は、解析対象の各要素の整合質量行列を取得し、前記解析対象の各要素の気泡関数に基づいて、前記第２の取得工程によって取得された各要素の整合質量行列を対角化した、各要素の対角質量行列を生成し、解析対象の既知解析物理量と、生成された各要素の対角質量行列と、に基づいて、前記解析対象の挙動を解析することを特徴とする。【００２１】
また、各要素における気泡関数の積分値を、各要素の整合質量行列に代入することにより、各要素の対角質量行列を生成することとしてもよい。さらに、生成された要素の対角質量行列に基づいて、前記解析対象全域の対角質量行列を算出し、前記解析対象全域の対角質量行列の逆行列を算出し、解析対象の既知解析物理量と、前記解析対象全域の対角質量行列と、前記逆行列と、に基づいて、前記解析対象の挙動を解析することとしてもよい。【００２２】
これらの発明によれば、気泡関数要素を使用した有限要素解析において質量行列の近似を行わずに、質量行列が対角行列となる下記条件式（１４）を満たす気泡関数を用いる。
【数６】
0004729767.tif 000007
【００２３】
上記式（１５−１）〜式（１５−３）は、気泡関数要素の定式化をおこなうために必要となる積分であり、<・,・>Ωeは要素領域Ωeの積分、Ａeは要素領域Ωeの面積（体積）を示している。これにより、各要素の整合質量行列内の成分を足し合わせることによって近似された集中質量行列を用いずに、気泡関数要素の基底（形状関数）が直交する条件を導入して、各要素の高精度な対角質量行列を用いて解析することができる。また、解析対象全域の対角質量行列やその逆行列も簡単に算出することができ、解析対象の挙動を、未知の解析対象物理量を用いて高精度に解析することができる。
【発明の効果】
【００２４】
本発明にかかる直交基底気泡関数要素数値解析方法、直交基底気泡関数要素数値解析プログラム、および直交基底気泡関数要素数値解析装置によれば、簡単かつ高精度に解析処理をおこなうことにより、解析対象の挙動解析の信頼性の向上を図ることができるという効果を奏する。また、記憶容量の低減化および解析時間の短縮化など、効率的な解析手法（直交基底気泡関数要素数値解析方法）を実現することができるという効果を奏する。
【図面の簡単な説明】
【００２５】
【図１】図１は、この発明の実施の形態にかかる数値解析装置のハードウェア構成を示すブロック図である。
【図２】図２は、この発明の実施の形態にかかる数値解析装置の機能構成を示すブロック図である。
【図３】図３は、この発明の実施の形態にかかる数値解析装置のデータ格納部が格納している解析物理量情報テーブルを示す説明図である。
【図４】図４は、この発明の実施の形態にかかる数値解析装置における数値解析処理手順を示すフローチャートである。
【図５】図５は、２次元の気泡関数要素を示す説明図である。
【図６】図６は、３次元の気泡関数要素を示す説明図である。
【図７】図７は、直交基底を形成する三角形気泡関数の形状φBを示す説明図である。
【図８】図８は、直交基底を形成する三角形気泡関数の形状φB2を示す説明図である。
【図９】図９は、直交基底を形成する三角形気泡関数の形状φBを示す説明図である。
【図１０】図１０は、直交基底を形成する三角形気泡関数の形状φB2を示す説明図である。
【図１１】図１１は、２次元の気泡関数を示す説明図である。
【図１２】図１２は、３次元の気泡関数を示す説明図である。
【図１３】図１３は、Rotating Cone問題の解析における計算領域を示す説明図である。
【図１４】図１４は、Rotating Cone問題の解析における解析モデルのメッシュを示す説明図である。
【図１５】図１５は、Rotating Cone問題の解析における解析モデルの流速（流況）を示す説明図である。
【図１６】図１６は、Rotating Cone問題の解析における初期条件を示す鳥瞰図である。
【図１７】図１７は、Rotating Cone問題の解析における初期条件の等高線を示す説明図である。
【図１８】図１８は、気泡関数（式（４７））での整合質量行列を使用した計算結果を示す鳥瞰図である。
【図１９】図１９は、気泡関数（式（４７））での整合質量行列を使用した計算結果の等高線を示す説明図である。
【図２０】図２０は、気泡関数式（（４７））での集中質量行列を使用した計算結果を示す鳥瞰図である。
【図２１】図２１は、気泡関数（式（４７））での集中質量行列を使用した計算結果の等高線を示す説明図である。
【図２２】図２２は、気泡関数（式（４８））での整合質量行列を使用した計算結果を示す鳥瞰図である。
【図２３】図２３は、気泡関数（式（４８）)での整合質量行列を使用した計算結果の等高線を示す説明図である。
【図２４】図２４は、気泡関数（式（４８））での集中質量行列を使用した計算結果を示す鳥瞰図である。
【図２５】図２５は、気泡関数（式（４８）)での集中質量行列を使用した計算結果の等高線を示す説明図である。
【図２６】図２６は、気泡関数（式（４９））での対角質量行列を使用した計算結果を示す鳥瞰図である。
【図２７】図２７は、気泡関数（式（４９））での対角質量行列を使用した計算結果の等高線を示す説明図である。
【図２８】図２８は、Rotating Cone問題の解析における解析モデルのメッシュ（節点数：438413、要素数：2539200）を示す説明図である。
【図２９】図２９は、Rotating Cone問題の解析における鉛直位置＝０の断面の初期条件を示す鳥瞰図である。
【図３０】図３０は、Rotating Cone問題の解析における鉛直位置＝０の部分の初期条件の等高線を示す説明図である。
【図３１】図３１は、気泡関数（式（５０））での整合質量行列を使用した鉛直位置＝０の部分の計算結果を示す鳥瞰図である。
【図３２】図３２は、気泡関数（式（５０））での整合質量行列を使用した鉛直位置＝０の部分の計算結果の等高線を示す説明図である。
【図３３】図３３は、気泡関数（式（５０））での集中質量行列を使用した鉛直位置＝０の部分の計算結果を示す鳥瞰図である。
【図３４】図３４は、気泡関数（式（５０））での集中質量行列を使用した鉛直位置＝０の部分の計算結果の等高線を示す説明図である。
【図３５】図３５は、気泡関数（式（５１））での整合質量行列を使用した鉛直位置＝０の部分の計算結果を示す鳥瞰図である。
【図３６】図３６は、気泡関数（式（５１））での整合質量行列を使用した鉛直位置＝０の部分の計算結果の等高線を示す説明図である。
【図３７】図３７は、気泡関数（式（５１））での集中質量行列を使用した鉛直位置＝０の部分の計算結果を示す鳥瞰図である。
【図３８】図３８は、気泡関数（式（５１））での集中質量行列を使用した鉛直位置＝０の部分の計算結果の等高線を示す説明図である。
【図３９】図３９は、気泡関数（式（５２））での対角質量行列を使用した鉛直位置＝０の部分の計算結果を示す鳥瞰図である。
【図４０】図４０は、気泡関数（式（５２））での対角質量行列を使用した鉛直位置＝０の部分の計算結果の等高線を示す説明図である。
【図４１】図４１は、アイソパラメトリック座標系ｒを示す説明図（その１）である。
【図４２】図４２は、アイソパラメトリック座標系ｒを示す説明図（その２）である。
【図４３】図４３は、式（８１）、（９３）を用いた三角形気泡関数の形状φBを示す説明図である。
【図４４】図４４は、式（８１）、（９３）を用いた三角形気泡関数の形状φB2を示す説明図である。
【図４５】図４５は、式（１３２）を用いた三角形３レベル気泡関数の形状を示す説明図である。
【図４６】図４６は、式（８１）、（９３）を用いた三角形気泡関数と式（１３２）を用いた三角形３レベル気泡関数の積の形状を示す説明図である。
【図４７】図４７は、従来の２次元の気泡関数要素を示す説明図である。
【図４８】図４８は、従来の３次元の気泡関数要素を示す説明図である。
【図４９】図４９は、従来のRotating Cone問題の解析に用いる解析モデルを示す説明図である。
【図５０】図５０は、従来のRotating Cone問題の解析に用いる初期条件の等高線を示す説明図である。
【図５１】図５１は、従来の整合質量行列によって計算されたRotating Cone問題の解析結果の等高線を示す説明図である。
【図５２】図５２は、従来の集中質量行列によって計算されたRotating Cone問題の解析結果の等高線を示す説明図である。
【符号の説明】
【００２６】
２００  数値解析装置
２０１  データ格納部
２０２  第１の取得部
２０３  第２の取得部
２０４  生成部
２０５  第１の算出部
２０６  第２の算出部
２０７  解析部
【発明を実施するための最良の形態】
【００２７】
以下に添付図面を参照して、この発明にかかる直交基底気泡関数要素数値解析方法、直交基底気泡関数要素数値解析プログラムおよび直交基底気泡関数要素数値解析装置（以下、単に、「数値解析装置」という）の好適な実施の形態を詳細に説明する。【００２８】
（数値解析装置のハードウェア構成）
この実施の形態にかかる数値解析装置のハードウェア構成について説明する。図１は、この発明の実施の形態にかかる数値解析装置のハードウェア構成を示すブロック図である。数値解析装置は、ＣＰＵ１０１と、ＲＯＭ１０２と、ＲＡＭ１０３と、ＨＤＤ（ハードディスクドライブ）１０４と、ＨＤ（ハードディスク）１０５と、ＦＤＤ（フレキシブルディスクドライブ）１０６と、着脱可能な記録媒体の一例としてＦＤ（フレキシブルディスク）１０７と、ディスプレイ１０８と、Ｉ／Ｆ（インターフェース）１０９と、キーボード１１０と、マウス１１１と、スキャナ１１２と、プリンタ１１３と、を備えている。また、各構成部は、バス１００によってそれぞれ接続されている。【００２９】
ここで、ＣＰＵ１０１は、数値解析装置の全体の制御を司る。ＲＯＭ１０２は、ブートプログラムなどのプログラムを記憶している。ＲＡＭ１０３は、ＣＰＵ１０１のワークエリアとして使用される。ＨＤＤ１０４は、ＣＰＵ１０１の制御にしたがってＨＤ１０５に対するデータのリード／ライトを制御する。ＨＤ１０５は、ＨＤＤ１０４の制御で書き込まれたデータを記憶する。【００３０】
ＦＤＤ１０６は、ＣＰＵ１０１の制御にしたがってＦＤ１０７に対するデータのリード／ライトを制御する。ＦＤ１０７は、ＦＤＤ１０６の制御で書き込まれたデータを記憶したり、ＦＤ１０７に記憶されたデータを数値解析装置に読み取らせたりする。着脱可能な記録媒体として、ＦＤ１０７のほか、ＣＤ−ＲＯＭ（ＣＤ−Ｒ、ＣＤ−ＲＷ）、ＭＯ、ＤＶＤ（Digital Versatile Disk）、メモリーカードなどであってもよい。ディスプレイ１０８は、カーソル、アイコンあるいはツールボックスをはじめ、文書、画像、機能情報などのデータを表示する。このディスプレイ１０８は、たとえば、ＣＲＴ、ＴＦＴ液晶ディスプレイ、プラズマディスプレイ等を採用することができる。【００３１】
Ｉ／Ｆ１０９は、通信回線を通じてインターネットなどのネットワークに接続され、このネットワークを介して他の装置に接続される。そして、Ｉ／Ｆ１０９は、ネットワークと内部のインターフェースを司り、外部装置からのデータの入出力を制御する。Ｉ／Ｆ１０９には、たとえばモデムやＬＡＮアダプタなどを採用することができる。【００３２】
キーボード１１０は、文字、数字、各種指示などの入力のためのキーを備え、データの入力をおこなう。また、タッチパネル式の入力パッドやテンキーなどであってもよい。マウス１１１は、カーソルの移動や範囲選択、あるいはウィンドウの移動やサイズの変更などをおこなう。ポインティングデバイスとして同様に機能を備えるものであれば、トラックボールやジョイスティックなどであってもよい。【００３３】
スキャナ１１２は、画像を光学的に読み取り、数値解析装置内に画像データを取り込む。また、プリンタ１１３は、画像データや文書データを印刷する。プリンタ１１３には、たとえば、レーザプリンタやインクジェットプリンタを採用することができる。【００３４】
（数値解析装置の機能的構成）
この発明の実施の形態にかかる数値解析装置の機能的構成について説明する。図２はこの発明の実施の形態にかかる数値解析装置の機能構成を示すブロック図である。図２に示すように、数値解析装置２００は、データ格納部２０１と、第１の取得部２０２と、第２の取得部２０３と、生成部２０４と、第１の算出部２０５と、第２の算出部２０６と、解析部２０７と、から構成されている。【００３５】
データ格納部２０１は、解析対象の数値解析に用いる解析物理量情報テーブルを有する。ここで、解析物理量情報テーブルについて説明する。図３は、この発明の実施の形態にかかる解析物理量情報テーブルを示す説明図である。【００３６】
図３において、解析物理量情報テーブルには、解析対象についての、解析物理量ＩＤと、解析物理量名（汚染物質、熱、流体、構造など）と、既知解析物理量とが格納されている。また、既知解析物理量としては、既知物性値と、境界解析物理量と、初期解析物理量とが格納されている。【００３７】
解析物理量ＩＤおよび解析物理量名は、ユーザの入力操作によって指定される情報であり、解析物理量ＩＤや解析物理量名が指定されると、関連付けられている既知物性値、境界解析物理量および初期解析物理量を抽出することができる。【００３８】
既知物性値は、ユーザの入力操作によって解析物理量ＩＤまたは解析物理量名が指定された解析対象が、どのような物質であるかを決定する情報である。既知物性値としては、たとえば、拡散係数、密度、比熱、熱伝導係数、粘性係数、水底面摩擦係数、ヤング率、ポアソン比などが挙げられる。【００３９】
境界解析物理量とは、解析対象の各要素からなる解析領域（三角形や四面体に分割したメッシュ）がどのような条件下で解析がおこなわれるかという解析条件をあらわす情報である。境界解析物理量としては、たとえば、物質濃度、温度、流速、圧力、流量、水深、変位などが挙げられる。【００４０】
境界解析物理量については、解析対象がたとえば汚染物質や熱の場合、汚染物質や熱の発生（供給）源の部分で、物質濃度＝発生（供給）源値、温度＝発生（供給）源値となる。また、解析対象が流体で境界が壁ならば、流速＝０、流量＝０、解析対象が構造で境界が物体を固定ならば、変位＝０に設定することができる。定常問題（解析物理量が時間に依存せず時々刻々と変化しない問題）であれば、上述した既知物性値および境界解析物理量によって、解析することができる。【００４１】
非定常問題（解析物理量が時間に依存し時々刻々と変化する問題）の場合は、さらに初期条件となる初期解析物理量が必要になる。初期解析物理量は、解析開始時の解析物理量の値である。また、データ格納部２０１において、図示はしないが、初期設定情報として、時間方向に対して、逐次解析計算をおこなうために時間増分量、計算ステップ数などの情報も、解析物理量情報テーブルに格納されている。【００４２】
さらに、解析対象範囲に分割形成されるメッシュの分割幅、節点数、要素数、節点番号に対応した節点座標値および要素番号に対応した要素結合情報などのメッシュ情報も格納されている。このメッシュは、解析対象範囲が２次元空間である場合は、三角形の形状、解析対象範囲が３次元空間である場合は、四面体の有限な形状の要素となる。メッシュは、図１に示したスキャナ１１２により読み取ることによっても取り込むことができる。このデータ格納部２０１は、具体的には、例えば、図１に示したＲＡＭ１０３、ＨＤ１０５またはＦＤ１０７によりその機能を実現する。【００４３】
第１の取得部２０２は、解析対象の既知解析物理量を取得する。具体的には、入力装置（たとえば図１に示したキーボード１１０またはマウス１１１）を操作することにより、解析対象となる物理量の解析物理量ＩＤを選択し、データ格納部２０１から既知解析物理量を抽出する。この第１の取得部２０２は、具体的には、たとえば図１に示したＲＯＭ１０２、ＲＡＭ１０３、ＨＤ１０５、ＦＤ１０７などに格納されたプログラムをＣＰＵ１０１が実行することによって、また図１に示したＩ／Ｆ１０９によってその機能を実現する。【００４４】
第２の取得部２０３は、解析対象の各要素（メッシュ）の整合質量行列を取得する。第２の取得部２０３は、具体的には、たとえば、図１に示したＲＯＭ１０２、ＲＡＭ１０３、ＨＤ１０５、ＦＤ１０７などに格納されたプログラムをＣＰＵ１０１が実行することによって、また図１に示したＩ／Ｆ１０９によってその機能を実現する。【００４５】
生成部２０４は、解析対象の各要素の気泡関数に基づいて、第２の取得部２０３によって取得された各要素の整合質量行列を対角化した、各要素の対角質量行列を生成する。生成部２０４は、具体的には、たとえば、図１に示したＲＯＭ１０２、ＲＡＭ１０３、ＨＤ１０５、ＦＤ１０７などに格納されたプログラムをＣＰＵ１０１が実行することによって、また図１に示したＩ／Ｆ１０９によってその機能を実現する。【００４６】
第１の算出部２０５は、生成部２０４によって生成された各要素の対角質量行列に基づいて、解析対象全域の対角質量行列を算出する。第１の算出部２０５は、具体的には、たとえば、図１に示したＲＯＭ１０２、ＲＡＭ１０３、ＨＤ１０５、ＦＤ１０７などに格納されたプログラムをＣＰＵ１０１が実行することによって、また図１に示したＩ／Ｆ１０９によってその機能を実現する。【００４７】
第２の算出部２０６は、第１の算出部２０５によって算出された、解析対象全域の対角質量行列の逆行列を算出する。第２の算出部２０６は、具体的には、たとえば、図１に示したＲＯＭ１０２、ＲＡＭ１０３、ＨＤ１０５、ＦＤ１０７などに格納されたプログラムをＣＰＵ１０１が実行することによって、また図１に示したＩ／Ｆ１０９によってその機能を実現する。【００４８】
解析部２０７は、第１の取得部２０２によって取得された既知解析物理量と、生成部２０４によって生成された各要素の対角質量行列と、に基づいて、解析対象の挙動を解析する。具体的には、第１の取得部２０２によって取得された既知解析物理量と、第１の算出部２０５によって算出された解析対象全域の対角質量行列と、第２の算出部２０６によって算出された逆行列と、に基づいて、解析対象の挙動を解析する。なお、解析対象全域の対角質量行列を必要とせず、各要素の対角質量行列のみを用いて解析対象の挙動を解析する解法があるが（たとえば、O.C.Zienkiewicz,  R.L.Taylor,  S.J.Sherwin  and  J.Peiro,  "On  Discontinuous  Galerkin  Methods",  International  Journal  for  Numerical  Methods  in  Engineering,  Vol.58,  2003,  pp.1119-1148参照）、この場合には、生成部２０４と第１の算出部２０５を同一視して、第２の算出部２０６へ進むことも可能である。解析部２０７は、具体的には、たとえば、図１に示したＲＯＭ１０２、ＲＡＭ１０３、ＨＤ１０５、ＦＤ１０７などに格納されたプログラムをＣＰＵ１０１が実行することによって、また図１に示したＩ／Ｆ１０９によってその機能を実現する。【００４９】
本実施の形態にかかる数値解析装置２００の数値解析処理手順について説明する。図４はこの発明の実施の形態にかかる数値解析装置２００の数値解析処理手順を示すフローチャートである。まず、第１の取得部２０２により、解析対象の既知解析物理量を取得する（ステップＳ４０１）。つぎに、第２の取得部２０３により、各要素（要素レベル）の整合質量行列を取得する（ステップＳ４０２）。【００５０】
続いて、要素ごとに、気泡関数を積分し（ステップＳ４０３）、各要素（要素レベル）の整合質量行列に、ステップＳ４０３で積分された値を代入することによって、各要素（要素レベル）の対角質量行列を算出する（ステップＳ４０４）。つぎに、各要素（要素レベル）の対角質量行列の総和（重ね合せ）により、解析対象全域の対角質量行列を算出する（ステップＳ４０５）。【００５１】
そして、この対角質量行列の逆行列を算出し（ステップＳ４０６）、解析対象の既知解析物理量と、解析対象全域の対角質量行列と、その逆行列とに基づいて、解析対象の挙動を解析する。具体的には、たとえば、上記式（１０）〜式（１３）に代入することにより、解析することができる。【００５２】
数値解析装置２００の具体的な解析処理内容について説明する。
＜１．直交基底気泡関数要素＞
質量行列Ｍは、有限要素法による解析（有限要素解析）について、ある解析物理量を補間した関数に重み関数を掛け、対象領域で積分を行った場合に現れる。気泡関数要素を使用した場合には、解析対象となる任意領域を区分的に三角形（四面体）形状にて要素数Ｎeで分割した要素群を全体系に重ね合せる作業によって形成され、要素領域での積分を用いて下記式（１６）のようにあらわすことができる。【００５３】
【数７】
0004729767.tif 000008
【００５４】
式（１６）に対して、気泡関数要素の基底（形状関数）が直交するためには、下記式（１８）、式（２０）を満たす必要がある。なお、式（１８）の導入に際して、式（１７）を用いる。また、式（２０）の導入に際して、式（１９）を用いる。【００５５】
【数８】
0004729767.tif 000009
【００５６】
【数９】
0004729767.tif 000010
【００５７】
【数１０】
0004729767.tif 000011
【００５８】
【数１１】
0004729767.tif 000012
【００５９】
式（１８）、（２０）より以下の下記関係式（２１）を得る。【００６０】
【数１２】
0004729767.tif 000013
【００６１】
＜２．直交条件を満たす気泡関数＞
式（１８）、式（２０）を満たすことができる気泡関数として下記式（２２）で表される気泡関数を提案する。【００６２】
【数１３】
0004729767.tif 000014
【００６３】
式（１８）に式（２２）を代入することにより下記式（２３）を得る。【００６４】
【数１４】
0004729767.tif 000015
【００６５】
式（２０）に式（２２）を代入することにより下記式（２４）を得る。【００６６】
【数１５】
0004729767.tif 000016
【００６７】
式（２３）に式（２４）を代入することにより、最終的にα1を未知量とした下記式（２５）を得る。【００６８】
【数１６】
0004729767.tif 000017
【００６９】
式（２５）はα1に関する二次方程式なので、解の公式より下記式（２６）にて未知量α1を求めることができる。【００７０】
【数１７】
0004729767.tif 000018
【００７１】
＜３．直交基底を形成する具体的な気泡関数＞
直交基底を形成する具体的な気泡関数について説明する。図５は、２次元の気泡関数要素を示す説明図であり、図６は、３次元の気泡関数要素を示す説明図である。【００７２】
直交基底を形成する具体的な気泡関数を導出するため、図５および図６に示した三角形(四面体)の要素領域を、重心点を用いて３（４）つの小三角形（小四面体）ｗi,ｉ＝１・・・Ｎ＋１に分割するξ乗気泡関数（０＜ξ＜∞）を用いる。ξ乗気泡関数は小三角形(小四面体)ごとにアイソパラメトリック座標系｛ｒ,ｓ｝（｛ｒ,ｓ,ｔ｝）を用いて次の式（２７）および式（２８）のように定義される。【００７３】
【数１８】
0004729767.tif 000019
【００７４】
これらの気泡関数の要素領域における積分値は、以下の式（２９）によって求めることができる。【００７５】
【数１９】
0004729767.tif 000020
【００７６】
式（２２）を、ξ乗気泡関数を用いて下記式（３０）のように定義する。【００７７】
【数２０】
0004729767.tif 000021
【００７８】
式（３０）の各べき乗の値を下記式（３１）のように設定すると二次元（三角形）、三次元（四面体）の気泡関数に対してα1，α2は下記式（３２）〜式（３５）のようになる。【００７９】
【数２１】
0004729767.tif 000022
【００８０】
【数２２】
0004729767.tif 000023
【００８１】
式（３０）の各べき乗の値を下記式（３６）のように設定する（式（３１）とは異なった値）と二次元（三角形）、三次元（四面体）の気泡関数に対してα1, α2は下記式（３７）〜式（４０）のようになる。【００８２】
【数２３】
0004729767.tif 000024
【００８３】
【数２４】
0004729767.tif 000025
【００８４】
図７は、式（３２）のα1,α2を使用した直交基底を形成する三角形気泡関数の形状φBを示す説明図であり、図８は、式（３２）のα1,α2を使用した直交基底を形成する三角形気泡関数の形状φB2を示す説明図である。また、図９は、式（３７）のα1,α2を使用した直交基底を形成する三角形気泡関数の形状φBを示す説明図であり、図１０は、式（３７）のα1,α2を使用した直交基底を形成する三角形気泡関数の形状φB2を示す説明図である。【００８５】
式（３２）、式（３７）および図７〜図１０に示すように、式（２１）を満たす気泡関数（形状）は無数に存在するため、直交基底を形成する具体的な気泡関数やその形状には殆ど意味がなく、直交基底気泡関数要素の本質は気泡関数要素の基底が直交する条件式（２１）を発明した部分にある。【００８６】
＜４．直交基底気泡関数要素の要素レベルの質量行列＞
下記式（４１）〜式（４６）に、２次元（三角形）、３次元（四面体）における通常の気泡関数要素の要素レベルの整合質量行列、集中質量行列、直交基底気泡関数要素の要素レベルの質量行列（対角質量行列）を示す。【００８７】
【数２５】
0004729767.tif 000026
【００８８】
式（４３）の要素レベルの対角質量行列Ｍ(e)ijは、生成部２０４により、式（４１）の要素レベルの整合質量行列Ｍ(e)ijの<１，φB>Ωeおよび‖φB‖2Ωeに、気泡関数φBの積分値ＣＡeを代入することによって生成される。そして、第１の算出部２０５により、各要素の要素レベルの対角質量行列Ｍ(e)ijの総和（重ね合せ）をとることによって、解析対象全域の対角質量行列Ｍを算出することができる。また、第２の算出部２０６により、対角質量行列Ｍの逆行列Ｍ-1を算出することができる。そして、解析部２０７により、既知解析物理量ｕ、対角質量行列Ｍおよび逆行列Ｍ-1を、上述した式（１０）〜式（１３）に代入することによって、解析対象の挙動を解析することができる。【００８９】
【数２６】
0004729767.tif 000027
【００９０】
式（４６）の要素レベルの対角質量行列は、生成部２０４により、式（４４）の要素レベルの整合質量行列の<１，φB>Ωeおよび‖φB‖2Ωeに、気泡関数φBの積分値ＣＡeを代入することによって生成される。そして、第１の算出部２０５により、各要素の要素レベルの対角質量行列Ｍ(e)ijの総和（重ね合せ）をとることによって、解析対象全域の対角質量行列Ｍを算出することができる。また、第２の算出部２０６により、対角質量行列Ｍの逆行列Ｍ-1を算出することができる。そして、解析部２０７により、既知解析物理量ｕ、対角質量行列Ｍおよび逆行列Ｍ-1を、上述した式（１０）〜式（１３）に代入することによって、解析対象の挙動を解析することができる。【００９１】
＜５．気泡関数について＞
気泡関数について説明する。図１１は、２次元の気泡関数を示す説明図であり、図１２は、３次元の気泡関数を示す説明図である。気泡関数φBは、形状関数の定義（重心点が１、その他の点では０となる関数）を満たせば任意の関数を選択することができる。検証の比較計算として、最も多く使用されている以下に示す二つの気泡関数および直交基底を形成する気泡関数を用いる。図１１および図１２に示したアイソパラメトリック座標系｛ｒ,ｓ｝（｛ｒ,ｓ,ｔ｝）を用いて、式（４７）〜式（５２）のように定義される。【００９２】
【数２７】
0004729767.tif 000028
【００９３】
【数２８】
0004729767.tif 000029
【００９４】
＜６．Rotating Cone問題＞
発明手法の検証計算として、数値解析手法の性能を比較検討するための有名なベンチマーク問題であるRotating  Cone問題（物質濃度の移流問題）の解析をおこなう（非特許文献４参照）。図１３は、Rotating  Cone問題の解析における計算領域を示す説明図であり、図１４は、Rotating  Cone問題の解析における解析モデルのメッシュ（節点数：4921、要素数：9600）を示す説明図であり、図１５は、Rotating  Cone問題の解析における解析モデルの流速（流況）を示す説明図である。【００９５】
図１６は、Rotating Cone問題の解析における初期条件（初期濃度分布）を示す鳥瞰図であり、図１７は、Rotating Cone問題の解析における初期条件（初期濃度分布）の等高線を示す説明図である。このRotating Cone問題は、非定常移流方程式において図１３の計算領域に、図１４にて図１５の定常な一定流速（流況）を仮定し、図１６および図１７のようなある物理量（汚染物質濃度など）の移流状況を解析する問題である。【００９６】
実際の物理量（汚染物質濃度など）は、濃度が拡散しながら流況にそって運ばれるが、拡散の効果をなくして計算を行った場合（移流方程式を計算した場合）には、濃度が流れに乗って計算領域を周回したとき、物理量の分布が、理想的には、物理量を移流させる作用しかないので元の初期分布と一致しなければならいという問題である。この問題により数値解析を行えば、対象とした手法が持つ近似誤差を評価することができ、計算精度の検証ができる。【００９７】
図１８は、気泡関数（式（４７））での整合質量行列を使用した計算結果を示す鳥瞰図であり、図１９は、気泡関数（式（４７））での整合質量行列を使用した計算結果の等高線を示す説明図である。また、図２０は、気泡関数（式（４７））での集中質量行列を使用した計算結果を示す鳥瞰図であり、図２１は、気泡関数（式（４７））での集中質量行列を使用した計算結果の等高線を示す説明図である。【００９８】
図２２は、気泡関数（式（４８））での整合質量行列を使用した計算結果を示す鳥瞰図であり、図２３は、気泡関数（式４８）での整合質量行列を使用した計算結果の等高線を示す説明図である。また、図２４は、気泡関数（式（４８））での集中質量行列を使用した計算結果を示す鳥瞰図であり、図２５は、気泡関数（式（４８））での集中質量行列を使用した計算結果の等高線を示す説明図である。また、図２６は、気泡関数（式（４９））での対角質量行列を使用した計算結果を示す鳥瞰図であり、図２７は、気泡関数（式（４９））での対角質量行列を使用した計算結果の等高線を示す説明図である。【００９９】
解析方法としては、検証対象としたすべての気泡関数に対して、空間方向の離散化には数値解析の安定化を考慮した気泡関数要素安定化法（非特許文献３参照）を使用し、時間方向の離散化には４段解法を、微小時間Δｔはπ／４００を採用した。図１８、図１９、図２２および図２３は、気泡関数（式（４７））、（式（４８））での整合質量行列を使用した計算結果である。【０１００】
この計算結果に対して、気泡関数（式（４７）、式（４８））での集中質量行列を使用した計算結果、図２０、図２１、図２４および図２５は、物理量（円錐）の進行方向の後方に振動が発生しており、円錐の分布状況も崩れている結果となっている。これは、濃度の物理的に意味のない減衰、広がりが著しく、濃度が最大になるべき地点も大きくずれており、数値解析上、信頼性の低い計算結果となっている。【０１０１】
一方、気泡関数（式（４９））での対角質量行列を使用した計算結果は、円錐の進行方向の後方に振動が発生しておらず、その分布も崩れていなく、気泡関数（式（４７））、式４８））での整合質量行列を使用した計算結果と同等の計算精度を保っており、数値解析上、信頼性の高い計算結果が得られている。【０１０２】
＜７．Rotating Cone問題による３次元解析の検証＞
発明手法の３次元解析の検証として、２次元と同様にRotating Cone問題（物質濃度の移流問題）の解析をおこなう（非特許文献４参照）。図２８は、Rotating Cone問題の解析における解析モデルのメッシュ（節点数：438413、要素数：2539200）を示す説明図である。図２８は、図１３の２次元解析モデルに鉛直長さ＝２（−１〜１）を持たせたモデルである。解析モデルの流速（流況）および濃度分布は、鉛直方向には変化がなく２次元解析モデルの場合と同様である。【０１０３】
図２９は、Rotating Cone問題の解析における鉛直位置＝０の断面の初期条件（初期濃度分布）を示す鳥瞰図であり、図３０は、Rotating Cone問題の解析における鉛直位置＝０の部分の初期条件（初期濃度分布）の等高線を示す説明図である。【０１０４】
図３１は、気泡関数（式（５０））での整合質量行列を使用した鉛直位置＝０の部分の計算結果を示す鳥瞰図であり、図３２は、気泡関数（式（５０））での整合質量行列を使用した鉛直位置＝０の部分の計算結果の等高線を示す説明図である。また、図３３は、気泡関数（式（５０））での集中質量行列を使用した鉛直位置＝０の部分の計算結果を示す鳥瞰図であり、図３４は、気泡関数（式（５０））での集中質量行列を使用した鉛直位置＝０の部分の計算結果の等高線を示す説明図である。【０１０５】
図３５は、気泡関数（式（５１））での整合質量行列を使用した鉛直位置＝０の部分の計算結果を示す鳥瞰図であり、図３６は、気泡関数（式（５１））での整合質量行列を使用した鉛直位置＝０の部分の計算結果の等高線を示す説明図である。また、図３７は、気泡関数（式（５１））での集中質量行列を使用した鉛直位置＝０の部分の計算結果を示す鳥瞰図であり、図３８は、気泡関数（式（５１））での集中質量行列を使用した鉛直位置＝０の部分の計算結果の等高線を示す説明図である。また、図３９は、気泡関数（式（５２））での対角質量行列を使用した鉛直位置＝０の部分の計算結果を示す鳥瞰図であり、図４０は、気泡関数（式（５２））での対角質量行列を使用した鉛直位置＝０の部分の計算結果の等高線を示す説明図である。【０１０６】
解析方法としては、検証対象としたすべての気泡関数に対して、空間方向の離散化には数値解析の安定化を考慮した気泡関数要素安定化法（非特許文献３参照）を使用し、時間方向の離散化には４段解法を、微小時間Δｔはπ／６００を採用した。図３１、図３２、図３５および図３６は、気泡関数（式（５０）、式（５１））での整合質量行列を使用した計算結果である。【０１０７】
この計算結果に対して、気泡関数（式（５０）、式（５１））での集中質量行列を使用した計算結果、図３３、図３４、図３７および図３８は、物理量（円錐）の進行方向の後方に振動が発生しており、円錐の分布状況も崩れている結果となっている。これは、濃度の物理的に意味のない減衰、広がりが著しく、濃度が最大になるべき地点も大きくずれており、数値解析上、信頼性の低い計算結果となっている。【０１０８】
一方、気泡関数（式（５２））での対角質量行列を使用した計算結果は、円錐の進行方向の後方に振動が発生しておらず、その分布も崩れていなく、気泡関数（式（５０）、式（５１））での整合質量行列を使用した計算結果と同等の計算精度を保っており、数値解析上、信頼性の高い計算結果が得られている。【０１０９】
＜８．１次元の直交基底気泡関数要素と要素レベルの質量行列＞
以上、２次元、３次元の直交基底気泡関数要素を発明し、その効果を示した。有限要素法の数値計算では、実用上、問題となる解析は、殆どの場合、２次元と３次元であるが、直交基底気泡関数要素は、１次元でも定義することができる。ここでは、１次元の気泡関数要素、ξ乗気泡関数を示し、直交条件式（１４）（または式（２１））を１次元〜３次元まで空間次元数Ｎを用いて統一的に記述する。また、一次元（線気泡関数要素）の要素レベルの質量行列を示す。【０１１０】
図４１は、アイソパラメトリック座標系ｒを示す説明図（その１）である。図４１に示したアイソパラメトリック座標系ｒを、式（１）、（２）の表現で記述すると、１次元のΦT, ｕT,Ψαは下記式（５３）〜（５５）になる。【０１１１】
【数２９】
0004729767.tif 000030
【０１１２】
図４２は、アイソパラメトリック座標系ｒを示す説明図（その２）である。図４２に示したアイソパラメトリック座標系ｒにおいて、中間点を用いて、線要素領域をｗ1、ｗ2に分割し、気泡関数の要素領域での積分値が式（２９）となる１次元のξ乗気泡関数は下記式（５６）のように定義される。【０１１３】
【数３０】
0004729767.tif 000031
【０１１４】
式（５６）、（２７）、（２８）のξ（η）乗気泡関数（０＜η＜∞）は、下記式（５７）に示すように、気泡関数を空間で偏微分した積分値も統一的に記述でき、下記式（５８）、（５９）のような積分値の関係を持っている。【０１１５】
【数３１】
0004729767.tif 000032
【０１１６】
式（５７）のｘk,ｘl＝ｘ1・・・ｘNはそれぞれ空間方向を表し、式（５９）のαkdは実数を、ＫDは正の整数を示している。式（５９）は式（５８）より誘導している。Ｃ＝（Ｎ＋１）／（Ｎ＋２）とし、空間次元数Ｎを使用して直交基底気泡関数要素の条件式（１４）（または式（２１））を１次元〜３次元まで統一的に記述すると下記式（６０）になる。【０１１７】
【数３２】
0004729767.tif 000033
【０１１８】
下記式（６１）〜（６３）に、１次元（線）における通常の気泡関数要素の要素レベルの整合質量行列、集中質量行列、直交基底気泡関数要素の要素レベルの質量行列（対角質量行列）を示す。【０１１９】
【数３３】
0004729767.tif 000034
【０１２０】
＜９．拡張された直交基底気泡関数要素＞
気泡関数要素を使用した流体解析、構造解析、固有値解析などでは、一般に一次要素の形状関数Ψαの積分値に加えて、下記式（６４）〜（６６）の気泡関数の積分値が必要になる。【０１２１】
【数３４】
0004729767.tif 000035
【０１２２】
下記式（６７）と式（６４）より、下記の関係式（６８）を得る。【０１２３】
【数３５】
0004729767.tif 000036
【０１２４】
気泡関数として、下記式（６９）の関係を持っている気泡関数の場合には、必要となる気泡関数の積分値（６４）〜（６６）の式（６４）は、式（７０）に代用することができる。【０１２５】
【数３６】
0004729767.tif 000037
【０１２６】
式（６９）の関係式は、非特許文献３で現れるような特殊な気泡関数を導入しない限り、例えば式（４７）、（４８）、（５０）、（５１）などの一般的によく使用されている気泡関数で成り立つので、直交条件式（２１）の導出仮定で現れる式（１７）、（１９）についても使用しており、式（２２）は式（６９）を満たす気泡関数であると仮定している。式（２２）右辺の気泡関数として式（３０）を採用した場合には、式（５９）より式（２２）は式（６９）を満足していることがわかる。【０１２７】
直交基底気泡関数要素を使用した場合には、式（７０）、（６５）、（６６）において、式（７０）、（６５）は直交条件式（６０）よりわかるが、式（６６）は、直交条件により決定することができない。式（５６）、（２７）、（２８）のξ乗気泡関数を使用して、直交基底気泡関数を形成した場合、実際に、式（６６）は、下記式（７１）で表される。【０１２８】
【数３７】
0004729767.tif 000038
【０１２９】
式（７１）をみると、Ｄ（ξ1,ξ2,ξ3）の値が、変数ξ1、ξ2、ξ3の関数となっている。実際に、式（３１）、（３６）の値を用いて、式（７２）を計算（α1（ξ1,ξ2,ξ3）＝α1+（ξ1,ξ2,ξ3）を使用）してみると、下記式（７３）〜（７５）のように採用する変数によって計算される値が異なる。【０１３０】
【数３８】
0004729767.tif 000039
【０１３１】
本発明では、直交基底気泡関数要素の拡張として、直交条件式（６０）によって得られる積分値を利用して式（６６）の積分値を決定できる、下記式（７６）を新たに定義する。【０１３２】
【数３９】
0004729767.tif 000040
【０１３３】
基底が直交する条件式（６０）に、式（７６）の条件を追加した下記式（８０）を拡張された直交基底気泡関数要素の条件式とする。【０１３４】
【数４０】
0004729767.tif 000041
【０１３５】
＜１０．拡張された直交基底気泡関数要素を満たす気泡関数＞
式（７１）からわかるように、式（５６）、（２７）、（２８）のξ乗気泡関数を使用して直交基底気泡関数を形成した式（６６）の積分値は、変数ξ1、ξ2、ξ3の値によって決まる。本特許では、拡張された直交基底気泡関数要素の条件式（８０）を満足するような変数ξ1、ξ2、ξ3を考え、そのような変数が実際に存在することを示す。変数はξ1、ξ2、ξ3の三つあるが、その内の２変数を下記式（８１）のように固定し、ξ1の１変数の問題として、拡張された直交基底気泡関数要素の条件を満足する気泡関数を求めることを考える。【０１３６】
【数４１】
0004729767.tif 000042
【０１３７】
ξ1の１変数の問題として、下記式（８２）を定義する。【０１３８】
【数４２】
0004729767.tif 000043
【０１３９】
拡張された直交基底気泡関数要素の条件を満たす気泡関数が存在することは、実数の集合Ｒにおいて、ｆ（ξ1）＝０となる変数が存在することである。いま、下記式（８３）〜（８５）のような実数の集合Ｒの部分集合Ｓを考える。【０１４０】
【数４３】
0004729767.tif 000044
【０１４１】
ξ1∈Ｓでは、１次元〜３次元において、下記式（８６）となる。【０１４２】
【数４４】
0004729767.tif 000045
【０１４３】
従って、α1、α2は実数の範囲で値を持つ。また、ξ1∈Ｓにて１次元〜３次元では、それぞれ下記式（８７）および式（８８）が成り立つので、式（７２）は有界な値をとり、式（８２）は実変数の関数となる。【０１４４】
【数４５】
0004729767.tif 000046
【０１４５】
式（８２）と部分集合Ｓを使用すると下記式（８９）〜（９１）を得ることが出来る。【０１４６】
【数４６】
0004729767.tif 000047
【０１４７】
任意のξ1c∈Ｓに対して、１次元では、式（８９）の第１式が負、第２式が正で、第３式が正（単調増加）、２次元、３次元では、式（９０）、（９１）の第１式が正、第２式が負で、第３式が負（単調減少）であるので、１次元〜３次元において、部分集合Ｓの範囲には、それぞれｆ（ξ1）＝０の解が一意に存在する。【０１４８】
具体的に、部分集合Ｓの範囲内で、二分法を使用してｆ（ξ1）＝０の解を求めると、下記式（９２）〜（９４）の値が得られる。【０１４９】
【数４７】
0004729767.tif 000048
【０１５０】
図４３、図４４に、式（８１）、（９３）を用いた三角形気泡関数の形状φB、三角形気泡関数の形状φB2を示す。【０１５１】
または、下記式（９５）、（９６）を使用して、ｆ（ξ1）＝０の解の存在を示すことができる。【０１５２】
【数４８】
0004729767.tif 000049
【０１５３】
１次元〜３次元において、ｈを正の定数とし、任意のξ1c∈Ｓに対して下記式（９７）の結果を得る。【０１５４】
【数４９】
0004729767.tif 000050
【０１５５】
平均値の定理から、任意のξ1a, ξ1b∈Ｓに対して、下記式（９８）が得られる。【０１５６】
【数５０】
0004729767.tif 000051
【０１５７】
従って、縮小写像の原理（定理）より、式（９５）を反復的に逐次計算をさせる下記式（９９）にて得られる数列｛ξ1m｝は、部分集合Ｓの範囲内でｆ（ξ1）＝０を満足する一意の解に収束する。【０１５８】
【数５１】
0004729767.tif 000052
【０１５９】
なお、拡張された直交基底気泡関数要素の存在を示すために、式（８０）の三つの条件を満たす気泡関数として、三つの未知量をα1、α2、ξ1とし、適当な気泡関数φB1、φB2、φB3を採用しが、他にも、下記式（１００）の気泡関数φBのように、未知量をα1、α2、α3として、適当な気泡関数φB1、φB2、φB3、φB4を採用すれば、式（８０）の条件を満たす気泡関数が得られる。このように、拡張された直交基底気泡関数要素の条件式を満たす気泡関数は、唯一ではなく、幾つも存在するため、具体的な気泡関数やその形状には殆ど意味がなく、拡張された直交基底気泡関数要素で重要な意味を成すのは、式（８０）の条件となる気泡関数の積分値である。【０１６０】
【数５２】
0004729767.tif 000053
【０１６１】
＜１１．拡張された直交基底気泡関数要素の要素レベルの粘性項（拡散項）の行列＞
通常の気泡関数、直交基底気泡関数要素の数値解析において式（６６）が必要になる要素レベルの粘性項（拡散項）の行列を下記式（１０１）〜（１０３）に示す。【０１６２】
【数５３】
0004729767.tif 000054
【０１６３】
式（１０１）〜（１０３）に式（７６）を代入して整理した、拡張された直交基底気泡関数要素の要素レベルの粘性項（拡散項）の行列を下記式（１０４）〜（１０６）に示す。【０１６４】
【数５４】
0004729767.tif 000055
【０１６５】
以上、拡張された直交基底気泡関数要素の要素レベルの粘性項（拡散項）の行列（１０４）〜（１０６）を示したが、本発明に関する定式化やプログラムでは、行列（１０１）〜（１０３）に式（７６）を代入して整理を行わずに、そのまま行列（１０１）〜（１０３）と式（７６）を使用することも、もちろん可能である。【０１６６】
＜１２．直交基底気泡関数要素と正規化気泡関数＞
式（１）は、下記式（１０７）、（１０８）のように変形することができる。【０１６７】
【数５５】
0004729767.tif 000056
【０１６８】
式（１０８）は、気泡関数自体の積分値がＡe（１次元では各要素ｅの線、２次元では各要素ｅの面積、３次元では各要素ｅの体積）となるような気泡関数を使用して、式変形を行ったものであり、式（１１０）は正規化気泡関数と呼ばれている。直交基底気泡関数要素において式（１）の表現形式の他に、式（１０７）の表現形式や式（１０８）の表現形式を用いて、定式化やプログラミング、計算実行を行うことも可能である。また、式（１）、（１０７）、（１０８）の表現形式を用いて得られた有限要素方程式（変分方程式）では、静的縮約という操作によって要素ごとに重心点の自由度ｕB、ｕB’、ｕS’を消去することができるが、重心点の自由度を消去した有限要素方程式（変分方程式）を用いて、定式化やプログラミング、計算実行を行うことも可能である。直交基底気泡関数要素の気泡関数として正規化気泡関数を使用した積分値は下記式（１１１）、（１１２）のようになる。【０１６９】
【数５６】
0004729767.tif 000057
【０１７０】
さらに、拡張された直交基底気泡関数要素では、下記式（１１３）が得られる。【０１７１】
【数５７】
0004729767.tif 000058
【０１７２】
＜１３．直交基底気泡関数要素と安定化作用制御項＞
気泡関数要素の持つ数値不安定性を緩和させるために、下記式（１１４）の気泡関数の粘性項（安定化作用制御項）を用いて解析精度を向上させる方法がある。【０１７３】
【数５８】
0004729767.tif 000059
【０１７４】
式（１１４）の安定化作用制御パラメータν’（−∞＜ν’＜∞）は、粘性係数の次元（粘性係数に相当する変数を無次元化している場合は無次元）を持つパラメータであり、この値を適切に決定することによって高精度解法を実現する。この方法は、従来、気泡関数要素に通常のGalerkin法を使用して、得られた有限要素方程式に重心点のみの粘性項（安定化作用制御項）を付加していたが（松本純一,梅津剛,「気泡関数を用いた非圧縮性流体の有限要素法解析に関する一考察」,日本応用数理学会平成８年度年会,1996年,46頁−47頁、松本純一,梅津剛,川原睦人,「線形型気泡関数を用いた非圧縮性粘性流体解析と適応型有限要素法」,応用力学論文集（土木学会）,2巻,1999年8月,223頁−232頁参照）、現在では、重み関数に関してもう一つの気泡関数（３レベル気泡関数）を採用して重心点のみの粘性項（安定化作用制御項）を導入している（J.Matsumoto  and  M.Kawahara,  "Shape  Identification  for  Fluid-Structure  Interaction  Problem  Using  Improved  Bubble  Element",  International  Journal  of  Computational  Fluid  Dynamics,  Vol.15,  2001,  pp.33-45、非特許文献３参照）。３レベル気泡関数は、下記式（１１５）〜（１１７）を満たす気泡関数である。【０１７５】
【数５９】
0004729767.tif 000060
【０１７６】
式（６７）、（１１５）より下記式（１１８）を得る。【０１７７】
【数６０】
0004729767.tif 000061
【０１７８】
また、下記式（１１９）の関係を持っている３レベル気泡関数の場合には、条件式（１１５）〜（１１７）の式（１１５）は式（１２０）で代用できる。【０１７９】
【数６１】
0004729767.tif 000062
【０１８０】
式（１１９）の関係を持っている３レベル気泡関数を仮定し、上記の条件式（１２０）、（１１６）、（１１７）を満たす具体的な３レベル気泡関数として、下記式（１２１）を用いる。【０１８１】
【数６２】
0004729767.tif 000063
【０１８２】
式中のνは、粘性係数（拡散係数）である。式（１２１）を式（１２０）、（１１６）、（１１７）に代入して整理すると、下記式（１２２）を得る。【０１８３】
【数６３】
0004729767.tif 000064
【０１８４】
式（１２１）右辺の四つの気泡関数を具体的に下記式（１２３）のように定義する。【０１８５】
【数６４】
0004729767.tif 000065
【０１８６】
式（１２３）、（５９）より式（１２１）は式（１１９）を満足していることがわかる。式（１２３）の変数を下記式（１２４）のように設定する。【０１８７】
【数６５】
0004729767.tif 000066
【０１８８】
α1（ξ1,ξ2,ξ3）＝α1+（ξ1,ξ2,ξ3）として、式（３１）の変数を用いると、式（１２２）より下記式（１２５）〜（１２７）を得る。【０１８９】
【数６６】
0004729767.tif 000067
【０１９０】
α1（ξ1,ξ2,ξ3）＝α1+（ξ1,ξ2,ξ3）として、式（３６）の変数を用いると、式（１２２）より下記式（１２８）〜（１３０）を得る。【０１９１】
【数６７】
0004729767.tif 000068
【０１９２】
α1（ξ1,ξ2,ξ3）＝α1+（ξ1,ξ2,ξ3）として、式（９２）〜（９４）、式（８１）の変数を用いると、式（１２２）より下記式（１３１）〜（１３３）を得る。【０１９３】
【数６８】
0004729767.tif 000069
【０１９４】
図４５、図４６に、式（１３２）を用いた三角形３レベル気泡関数の形状（ν’＝νとした場合）、式（９３）、（８１）を用いた三角形気泡関数と式（１３２）を用いた三角形３レベル気泡関数の積の形状（ν’＝νとした場合）を示す。以上により、直交基底気泡関数要素を用いた場合にも、式（１２０）（あるいは式（１１５））、（１１６）、（１１７）を満たす３レベル気泡関数は存在し、安定化作用制御項（１１４）を導入できる。【０１９５】
なお、３レベル気泡関数の存在を示すため、式（１２０）、（１１６）、（１１７）の三つの条件を満足する三つの未知量δ1、δ2、δ3と四つの適当な気泡関数を持った３レベル気泡関数（１２１）を採用したが、この他にも、下記式（１３４）のように、安定化作用制御パラメータν’と同じ次元を持つパラメータを導入し、未知量をδ1、δ2として、式（１３４）の右辺に適当な三つの気泡関数を用いることによって、式（１２０）、（１１６）の条件を満たす３レベル気泡関数を求めることができる。そして、下記式（１３５）のように安定化作用制御パラメータν’を定義すれば、式（１１７）も導出でき、式（１２０）（あるいは式（１１５））、（１１６）、（１１７）となる３レベル気泡関数の存在を示すことが出来る。このように、３レベル気泡関数の条件式を満たす気泡関数は、唯一ではなく、幾つも存在するため、具体的な気泡関数やその形状には殆ど意味がなく、３レベル気泡関数で重要な意味を成すのは、安定化作用制御項を導くために必要な式（１２０）（あるいは式（１１５））、（１１６）、（１１７）の条件である。【０１９６】
【数６９】
0004729767.tif 000070
【０１９７】
＜１４．直交基底気泡関数要素と正規化気泡関数による安定化作用制御項＞
正規化気泡関数を使用した場合にも安定化作用制御項は導入することが出来、下記式（１３６）によって表される。【０１９８】
【数７０】
0004729767.tif 000071
【０１９９】
直交基底気泡関数要素で使用する気泡関数として正規化気泡関数を用いた場合、３レベル正規化気泡関数は下記式（１３７）と定めることができ、その積分値は下記式（１３８）〜（１４２）のように記述できる。【０２００】
【数７１】
0004729767.tif 000072
【０２０１】
上記式（１３８）〜（１４２）より、式（１１５）〜（１１７）、（１１９）、（１２０）を満足する３レベル気泡関数が存在すれば、式（１３８）〜（１４２）も自動的に満たすので３レベル正規化気泡関数も同時に存在することがわかる。【０２０２】
＜１５．直交基底気泡関数要素と気泡関数の安定化作用行列＞
式（１０７）の表現形式を用いた気泡関数要素の重心点の自由度ｕB’は、静的縮約という操作によって要素ごとに消去できる。その結果得られる有限要素方程式は、安定化有限要素法で導かれる有限要素方程式と密接な関係がある（非特許文献２参照）。静的縮約によって得られる気泡関数の安定化作用行列は、一般的に下記式（１４３）の形式によって表される。【０２０３】
【数７２】
0004729767.tif 000073
【０２０４】
Ｌは、求めるべき未知解析物理量の数である。式（１４４）は、式（６９）を満たす気泡関数を使用した場合、式（７０）、（６５）、（６６）の気泡関数の積分値から成る行列である。安定化作用制御パラメータν’は、下記式（１４５）の関係が満たされるように決定する。【０２０５】
【数７３】
0004729767.tif 000074
【０２０６】
式（１４６）は解析精度が向上することが期待できる安定化行列である。式（１４５）の関係を満たすような安定化作用制御パラメータν’を求めるために、下記式（１４８）の評価関数を最小にするような安定化作用制御パラメータν’を決定する問題を考える。【０２０７】
【数７４】
0004729767.tif 000075
【０２０８】
式（１４８）を最小にするような安定化作用制御パラメータν’を求めるために、式（１４８）の停留条件である下記関係式（１５０）を用いて、安定化作用制御パラメータν’を求める。【０２０９】
【数７５】
0004729767.tif 000076
【０２１０】
３レベル気泡関数は重み関数の中で導入しているので、求めるべき未知解析物理量の偏微分方程式ごとに異なった３レベル気泡関数を使用することができる。従って、安定化作用制御パラメータν’は、下記式（１５１）の３レベル気泡関数を用いると、求めるべき未知解析物理量の数Ｌだけ、安定化作用制御パラメータνm’を定義することができる。【０２１１】
【数７６】
0004729767.tif 000077
【０２１２】
式（１５１）の係数δ1、δ2、δ3は、式（１２２）によって決定することが出来、式（１５１）は式（１２０）、（１１６）を満たす。また、式（１１７）の代わりに下記式（１５２）の安定化作用制御項を得る。【０２１３】
【数７７】
0004729767.tif 000078
【０２１４】
式（１５２）の安定化作用制御項を用いた場合、式（１４４）は下記式（１５３）のようになる。【０２１５】
【数７８】
0004729767.tif 000079
【０２１６】
式（１４５）の関係を満たすような安定化作用制御パラメータνm’を求めるために、下記式（１５４）の評価関数を最小にするような安定化作用制御パラメータνm’を決定する問題を考える。【０２１７】
【数７９】
0004729767.tif 000080
【０２１８】
式（１５４）を最小にするような安定化作用制御パラメータνm’を求めるために、式（１５４）の停留条件である下記関係式（１５６）を用いて、安定化作用制御パラメータνm’を求める。【０２１９】
【数８０】
0004729767.tif 000081
【０２２０】
式（１５５）はベクトルであり、ベクトルの成分を下記式（１５７）のように表現する。【０２２１】
【数８１】
0004729767.tif 000082
【０２２２】
式（１５７）をみるとわかるように、安定化作用制御パラメータνm’は各ベクトル成分で独立に定義されているので、式（１４５）の関係を満たすような安定化作用制御パラメータνm’を決定する問題として、式（１５７）の成分ｍごとに、下記式（１５８）の評価関数を用いた最小化問題を考えることができる。【０２２３】
【数８２】
0004729767.tif 000083
【０２２４】
式（１５８）を最小にするような安定化作用制御パラメータνm’を求めるために、式（１５８）の停留条件である下記関係式（１５９）を用いて、安定化作用制御パラメータνm’を求める。【０２２５】
【数８３】
0004729767.tif 000084
【０２２６】
式（１５７）は、安定化作用制御パラメータνm’が各ベクトル成分で独立に定義されているので、式（１５４）、（１５８）の停留条件から導かれる安定化作用制御パラメータνm’を求める式（１５６）、（１５９）は、結果的に同じ式になる。【０２２７】
他の方法として、下記式（１６０）の関係を満たすように、安定化作用制御パラメータν’を決定する問題も考えられる。【０２２８】
【数８４】
0004729767.tif 000085
【０２２９】
式（１６０）の関係を満たすような安定化作用制御パラメータν’を求めるために、下記式（１６１）の評価関数を最小にするような安定化作用制御パラメータν’を決定する問題を考える。【０２３０】
【数８５】
0004729767.tif 000086
【０２３１】
式（１６１）を最小にするような安定化作用制御パラメータν’を求めるために、式（１６１）の停留条件である下記関係式（１６３）を用いて、安定化作用制御パラメータν’を求める。【０２３２】
【数８６】
0004729767.tif 000087
【０２３３】
安定化作用制御パラメータν’は、式（１５１）の３レベル気泡関数を用いると、式（１５２）の安定化作用制御項が導出され、求めるべき未知解析物理量の数Ｌだけ、安定化作用制御パラメータνm’を定義することができる。式（１６０）の関係を満たすような安定化作用制御パラメータνm’を求めるために、下記式（１６４）の評価関数を最小にするような安定化作用制御パラメータνm’を決定する問題を考える。【０２３４】
【数８７】
0004729767.tif 000088
【０２３５】
式（１６４）を最小にするような安定化作用制御パラメータνm’を求めるために、式（１６４）の停留条件である下記関係式（１６６）を用いて、安定化作用制御パラメータνm’を求める。【０２３６】
【数８８】
0004729767.tif 000089
【０２３７】
式（１６５）はベクトルであり、ベクトルの成分を下記式（１６７）のように表現する。【０２３８】
【数８９】
0004729767.tif 000090
【０２３９】
式（１６７）をみるとわかるように、安定化作用制御パラメータνm’は各ベクトル成分で独立に定義されているので、式（１６０）の関係を満たすような安定化作用制御パラメータνm’を決定する問題として、式（１６７）の成分ｍごとに、下記式（１６８）の評価関数を用いた最小化問題を考えることができる。【０２４０】
【数９０】
0004729767.tif 000091
【０２４１】
式（１６８）を最小にするような安定化作用制御パラメータνm’を求めるために、式（１６８）の停留条件である下記関係式（１６９）を用いて、安定化作用制御パラメータνm’を求める。【０２４２】
【数９１】
0004729767.tif 000092
【０２４３】
式（１６７）は、安定化作用制御パラメータνm’が各ベクトル成分で独立に定義されているので、式（１６４）、（１６８）の停留条件から導かれる安定化作用制御パラメータνm’を求める式（１６６）、（１６９）は、結果的に同じ式になる。【０２４４】
＜１６．直交基底気泡関数要素と正規化気泡関数の安定化作用行列＞
式（１０８）の表現形式を用いた気泡関数要素の重心点の自由度ｕS’は、静的縮約という操作によって要素ごとに消去できる。静的縮約によって得られる正規化気泡関数の安定化作用行列は、一般的に下記式（１７０）の形式によって表される。【０２４５】
【数９２】
0004729767.tif 000093
【０２４６】
式（１７０）の第二項は、式（１４１）の関係を満たす気泡関数を使用した場合、式（１１１）〜（１１３）各第一項の正規化気泡関数の積分値から成る行列であり、式（７０）、（６５）、（６６）で表した場合には第三項のように記述できる。安定化作用制御パラメータν’は、下記式（１７２）の関係が満たされるように決定する。【０２４７】
【数９３】
0004729767.tif 000094
【０２４８】
式（１７２）の関係を満たすような安定化作用制御パラメータν’を求めるために、下記式（１７４）の評価関数を最小にするような安定化作用制御パラメータν’を決定する問題を考える。【０２４９】
【数９４】
0004729767.tif 000095
【０２５０】
式（１７４）を最小にするような安定化作用制御パラメータν’を求めるために、式（１７４）の停留条件である下記関係式（１７６）を用いて、安定化作用制御パラメータν’を求める。【０２５１】
【数９５】
0004729767.tif 000096
【０２５２】
３レベル正規化気泡関数は重み関数の中で導入しているので、求めるべき未知解析物理量の偏微分方程式ごとに異なった３レベル正規化気泡関数を使用することができる。従って、安定化作用制御パラメータν’は、下記式（１７７）の３レベル正規化気泡関数を用いると、求めるべき未知解析物理量の数Ｌだけ、安定化作用制御パラメータνm’を定義することができる。【０２５３】
【数９６】
0004729767.tif 000097
【０２５４】
式（１７７）の係数δ1、δ2、δ3は、式（１２２）によって決定することが出来、式（１７７）は式（１４２）、（１３９）を満たす。また、式（１４０）の代わりに下記式（１７８）の安定化作用制御項を得る。【０２５５】
【数９７】
0004729767.tif 000098
【０２５６】
式（１７８）の安定化作用制御項を用いた場合、式（１７１）は下記式（１７９）のようになる。【０２５７】
【数９８】
0004729767.tif 000099
【０２５８】
式（１７９）の第一項は、式（１１１）〜（１１３）各第一項の正規化気泡関数の積分値から成る行列であり、式（７０）、（６５）（６６）で表した場合には第二項のように記述できる。式（１７２）の関係を満たすような安定化作用制御パラメータνm’を求めるために、下記式（１８０）の評価関数を最小にするような安定化作用制御パラメータνm’を決定する問題を考える。【０２５９】
【数９９】
0004729767.tif 000100
【０２６０】
式（１８０）を最小にするような安定化作用制御パラメータνm’を求めるために、式（１８０）の停留条件である下記関係式（１８２）を用いて、安定化作用制御パラメータνm’を求める。【０２６１】
【数１００】
0004729767.tif 000101
【０２６２】
式（１８１）はベクトルであり、ベクトルの成分を下記式（１８３）のように表現する。【０２６３】
【数１０１】
0004729767.tif 000102
【０２６４】
式（１８３）をみるとわかるように、安定化作用制御パラメータνm’は各ベクトル成分で独立に定義されているので、式（１７２）の関係を満たすような安定化作用制御パラメータνm’を決定する問題として、式（１８３）の成分ｍごとに、下記式（１８４）の評価関数を用いて最小化問題を考えることができる。【０２６５】
【数１０２】
0004729767.tif 000103
【０２６６】
式（１８４）を最小にするような安定化作用制御パラメータνm’を求めるために、式（１８４）の停留条件である下記関係式（１８５）を用いて、安定化作用制御パラメータνm’を求める。【０２６７】
【数１０３】
0004729767.tif 000104
【０２６８】
式（１８３）は、安定化作用制御パラメータνm’が各ベクトル成分で独立に定義されているので、式（１８０）、（１８４）の停留条件から導かれる安定化作用制御パラメータνm’を求める式（１８２）、（１８５）は、結果的に同じ式になる。【０２６９】
他の方法として、下記式（１８６）の関係を満たすように、安定化作用制御パラメータν’を決定する問題も考えられる。【０２７０】
【数１０４】
0004729767.tif 000105
【０２７１】
式（１８６）の関係を満たすような安定化作用制御パラメータν’を求めるために、下記式（１８７）の評価関数を最小にするような安定化作用制御パラメータν’を決定する問題を考える。【０２７２】
【数１０５】
0004729767.tif 000106
【０２７３】
式（１８７）を最小にするような安定化作用制御パラメータν’を求めるために、式（１８７）の停留条件である下記関係式（１８９）を用いて、安定化作用制御パラメータν’を求める。【０２７４】
【数１０６】
0004729767.tif 000107
【０２７５】
安定化作用制御パラメータν’は、式（１７７）の３レベル正規化気泡関数を用いると、式（１７８）の安定化作用制御項が導出され、求めるべき未知解析物理量の数Ｌだけ、安定化作用制御パラメータνm’を定義することができる。式（１８６）の関係を満たすような安定化作用制御パラメータνm’を求めるために、下記式（１９０）の評価関数を最小にするような安定化作用制御パラメータνm’を決定する問題を考える。【０２７６】
【数１０７】
0004729767.tif 000108
【０２７７】
式（１９０）を最小にするような安定化作用制御パラメータνm’を求めるために、式（１９０）の停留条件である下記関係式（１９２）を用いて、安定化作用制御パラメータνm’を求める。【０２７８】
【数１０８】
0004729767.tif 000109
【０２７９】
式（１９１）はベクトルであり、ベクトルの成分を下記式（１９３）のように表現する。【０２８０】
【数１０９】
0004729767.tif 000110
【０２８１】
式（１９３）をみるとわかるように、安定化作用制御パラメータνm’は各ベクトル成分で独立に定義されているので、式（１８６）の関係を満たすような安定化作用制御パラメータνm’を決定する問題として、式（１９３）の成分ｍごとに、下記式（１９４）の評価関数を用いた最小化問題を考えることができる。【０２８２】
【数１１０】
0004729767.tif 000111
【０２８３】
式（１９４）を最小にするような安定化作用制御パラメータνm’を求めるために、式（１９４）の停留条件である下記関係式（１９５）を用いて、安定化作用制御パラメータνm’を求める。【０２８４】
【数１１１】
0004729767.tif 000112
【０２８５】
式（１９３）は、安定化作用制御パラメータνm’が各ベクトル成分で独立に定義されているので、式（１９０）、（１９４）の停留条件から導かれる安定化作用制御パラメータνm’を求める式（１９２）、（１９５）は、結果的に同じ式になる。【０２８６】
＜１７．気泡関数と正規化気泡関数の安定化作用制御パラメータの関係＞
直交基底気泡関数要素で、式（１）、（１０７）の表現形式を採用した場合には、式（７０）、（６５）、（６６）の積分を使用して、式（１５０）、（１５６）、（１５９）、（１６３）、（１６６）、（１６９）のいずれかを採用し、式（１０８）の表現形式を採用した場合には、式（１１１）〜（１１３）各第一項の積分を使用して、式（１７６）、（１８２）、（１８５）、（１８９）、（１９２）、（１９５）のいずれかを採用すれば安定化作用制御パラメータを計算することが出来る。いま、式（１７６）、（１８２）、（１８５）、（１８９）、（１９２）、（１９５）を、式（７０）、（６５）、（６６）の積分を使用して記述すると下記式（１９６）〜（２０１）のようになる。【０２８７】
【数１１２】
0004729767.tif 000113
【０２８８】
式（７０）、（６５）、（６６）の積分を使用した式（１９６）〜（２０１）の第二項より、安定化作用制御パラメータを求めるための関係式は、式（１５０）、（１５６）、（１５９）、（１６３）、（１６６）、（１６９）の安定化作用制御パラメータを求めるための関係式と同じに（等しく）なる。従って、式（１５０）、（１５６）、（１５９）、（１６３）、（１６６）、（１６９）より式（７０）、（６５）、（６６）の気泡関数の積分を使用して求めた安定化作用制御パラメータと、式（１７６）、（１８２）、（１８５）、（１８９）、（１９２）、（１９５）より式（１１１）〜（１１３）各第一項の正規化気泡関数の積分を使用して求めた安定化作用制御パラメータは、同じ値（等しい値）になることがわかる。【０２８９】
上記より、１）気泡関数、気泡関数自体の積分値がＡeとなる正規化気泡関数、２）気泡関数の持つ数値不安定を緩和させるために導入する安定化作用制御項、および３レベル気泡関数、３レベル正規化気泡関数、３）安定化作用制御パラメータの計算で用いる安定化作用行列、の定式化やプログラミングにおいて、式（７０）、（６５）、（６６）の気泡関数の積分値、あるいは式（１１１）〜（１１３）各第一項の正規化気泡関数の積分値が必要になる。拡張された直交基底気泡関数要素は、必要となる積分値が全て事前に（複雑な積分を行うことなく）式（８０）によって与えられているので、整合質量行列と同等の解析精度を持つ対角質量行列の作成や気泡関数の持つ数値不安定性を改良する安定化作用制御項の導入で現れる一連の作業において、非常に合理的に処理を行うことが出来る。【０２９０】
直交基底気泡関数要素は、２次元、３次元の解析結果からわかるように、記憶容量、計算時間の効率の良い解法（例えば、４段解法）を採用しながら、信頼性のある結果を得ることができる。さらに、拡張された直交基底気泡関数要素の条件式（８０）により、インクジェットの液滴挙動解析、都市や河川の氾濫解析、地震等で発生する津波の解析のみならず、構造物の応力解析、自動車のエンジン固有振動解析などで必要となる気泡関数要素の積分値や気泡関数の持つ数値不安定性を改良するために必要となる処理を合理的、かつ統一的に扱うことができるので、上記のような多岐にわたる数値解析への適用が容易に行えるという効果をもたらす。【０２９１】
以上説明したように、本実施の形態によれば、気泡関数要素を使用した数値解析において、記憶容量、計算時間などの計算効率の良い解析手法（直交基底気泡関数要素数値解析方法）を実現することができる。【０２９２】
なお、本実施の形態で説明した直交基底気泡関数要素数値解析方法は、予め用意されたプログラムをパーソナル・コンピュータ、ワークステーション、ＰＣクラスタ、スーパーコンピュータ等のコンピュータで実行することにより実現できる。このプログラムは、ハードディスク、フレキシブルディスク、ＣＤ−ＲＯＭ、ＭＯ、ＤＶＤ等のコンピュータで読み取り可能な記録媒体に記録され、コンピュータにて記録媒体から読み出されることによって実行される。またこのプログラムは、インターネット等のネットワークを介して配布することが可能な伝送媒体であってもよい。
【産業上の利用可能性】
【０２９３】
以上のように、本発明にかかる直交基底気泡関数要素数値解析方法、直交基底気泡関数要素数値解析プログラムおよび直交基底気泡関数要素数値解析装置は、線、三角形および四面体形状のメッシュを使用した有限要素法に基づく数値解析に有用である。
【図１】
0004729767.tif 000116
【図２】
0004729767.tif 000117
【図３】
0004729767.tif 000118
【図４】
0004729767.tif 000119
【図５】
0004729767.tif 000120
【図６】
0004729767.tif 000121
【図７】
0004729767.tif 000122
【図８】
0004729767.tif 000123
【図９】
0004729767.tif 000124
【図１０】
0004729767.tif 000125
【図１１】
0004729767.tif 000126
【図１２】
0004729767.tif 000127
【図１３】
0004729767.tif 000128
【図１４】
0004729767.tif 000129
【図１５】
0004729767.tif 000130
【図１６】
0004729767.tif 000131
【図１７】
0004729767.tif 000132
【図１８】
0004729767.tif 000133
【図１９】
0004729767.tif 000134
【図２０】
0004729767.tif 000135
【図２１】
0004729767.tif 000136
【図２２】
0004729767.tif 000137
【図２３】
0004729767.tif 000138
【図２４】
0004729767.tif 000139
【図２５】
0004729767.tif 000140
【図２６】
0004729767.tif 000141
【図２７】
0004729767.tif 000142
【図２８】
0004729767.tif 000143
【図２９】
0004729767.tif 000144
【図３０】
0004729767.tif 000145
【図３１】
0004729767.tif 000146
【図３２】
0004729767.tif 000147
【図３３】
0004729767.tif 000148
【図３４】
0004729767.tif 000149
【図３５】
0004729767.tif 000150
【図３６】
0004729767.tif 000151
【図３７】
0004729767.tif 000152
【図３８】
0004729767.tif 000153
【図３９】
0004729767.tif 000154
【図４０】
0004729767.tif 000155
【図４１】
0004729767.tif 000156
【図４２】
0004729767.tif 000157
【図４３】
0004729767.tif 000158
【図４４】
0004729767.tif 000159
【図４５】
0004729767.tif 000160
【図４６】
0004729767.tif 000161
【図４７】
0004729767.tif 000162
【図４８】
0004729767.tif 000163
【図４９】
0004729767.tif 000164
【図５０】
0004729767.tif 000165
【図５１】
0004729767.tif 000166
【図５２】
0004729767.tif 000167