JP4729767B2 - 直交基底気泡関数要素数値解析方法、直交基底気泡関数要素数値解析プログラムおよび直交基底気泡関数要素数値解析装置 - Google Patents

Description

この発明は、気泡関数要素を用いた有限要素法による解析（有限要素解析）について、計算効率の良い対角項のみとなる質量行列を用いて、信頼性の高い数値シミュレーションをおこなうための直交基底気泡関数要素数値解析方法、直交基底気泡関数要素数値解析プログラム、および直交基底気泡関数要素数値解析装置に関する。

従来の気泡関数要素について説明する。図４７は、従来の２次元の気泡関数要素を示す説明図であり、図４８は、従来の３次元の気泡関数要素を示す説明図である。図４７および図４８のように、三角形（四面体）要素を用いた気泡関数要素は、各要素において三角形（四面体）を形成する３（４）点と重心点の４（５）つの節点を用いて、アイソパラメトリック座標系[ｒ,ｓ]（｛ｒ,ｓ,ｔ｝）で下記式（１）のように表される（たとえば、下記非特許文献１、非特許文献２、非特許文献３参照）。

式（１）のΦα,φ_Bは、気泡関数要素の形状関数、ｕα,ｕ_Bは三角形（四面体）の各頂点の値（解析物理量）、重心点の値（解析物理量）、Ｎは空間次元数を示している。ベクトル形式で記述すれば、形状関数は下記式（３）〜式（６）になる。

式（２）のΨαは、２次元および３次元の一次要素を用いた形状関数であり、下記式（７）、（８）であらわされる。

形状関数φ_Bは気泡関数と呼ばれている。気泡関数は要素境界上においてその値が０となり、重心点で値が１となるように要素毎に定義される。非定常問題において、空間方向の離散化に気泡関数要素を用いた有限要素方程式は下記式（９）のようにあらわすことができる。

式（９）のｕは求めるべき未知解析物理量（汚染物質濃度、温度、流量、水深、流速、圧力、変位など）であり、Ｍは質量行列、Ｆ（ｕ）は時間微分項以外をまとめた項である。式（９）の時間方向の離散化として、テイラー展開に基づいた４段解法は下記式（１０）〜式（１３）のように表される（たとえば、下記非特許文献４参照）。

式（１０）〜式（１３）の上付き添え字ｎは現在時刻ｎでの既知解析物理量を表し、ｎ＋１は時刻ｎから微小時間Δｔ経過後の未知解析物理量を表している。

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レベル有限要素法」,応用力学論文集（土木学会）,７巻,２００４年８月,３３９頁−３４６頁 畑中勝守,「多段階有限要素法による非圧縮粘性流体の順・逆解析に関する計算力学的研究」,中央大学博士論文,１９９３年３月

ここで、上述した式（１０）〜式（１３）に示したように、４段解法を使用して、この既知解析物理量ｕⁿから未知解析物理量ｕⁿ⁺¹を求めるためには、質量行列の逆行列が必要になる。図４９は、従来のRotating Cone問題の解析に用いる解析モデルを示す説明図であり、図５０は、従来のRotating Cone問題の解析に用いる初期条件の等高線を示す説明図であり、図５１は、従来の整合質量行列によって計算されたRotating Cone問題の解析結果の等高線を示す説明図である（非特許文献４参照）。

図４９の解析モデル４９００を初期状態（０周時）とする。そして、質量行列の逆行列を用いて、初期状態（０周時）から、所定周、たとえば、５周回転させると、図５０に示した解析モデル４９００の初期状態における等高線モデル５０００は、図５１に示したような解析結果（等高線モデル５１００）となる。

上述した質量行列は、時間方向の離散化に気泡関数要素を用いているため、一般的に疎な分布行列（整合質量行列）となる。したがって、この分布行列（整合質量行列）の逆行列を求めることは、数値解析上多くの記憶容量および計算時間が必要となり、装置本体のコストがかかるとともに、解析処理が遅延するという問題があった。

この問題を解決するために、通常は、質量行列の各行の成分を足し合わせて（集中化させて）、対角項のみに成分をもたせた近似行列（集中質量行列）が使用される。集中質量行列を用いた場合には、行列の成分が対角項のみであるので、逆行列は各対角成分を逆数にした行列になり、数値解析上、近似なしの整合質量行列およびその逆行列を使用する場合に比べて、非常に少ない記憶容量、計算時間で解析を実行することができる。

しかしながら、上述した集中質量行列を使用した場合には、集中質量行列が元の質量行列と同一の行列ではなく近似行列となる。したがって、図４９に示した解析モデル４９００の初期状態（０周時）から、所定周、たとえば、５周回転させると、図５０に示した、解析モデル４９００の初期状態における等高線モデル５０００は、図５２に示したような解析結果（等高線モデル５２００）となる。このように、解析モデル４９００の５周回転後の等高線モデル５２００は、初期状態における等高線モデル５０００や図５１に示した等高線モデル５１００に対して、大幅に変形しているため計算精度が悪く、解析結果の信頼性が低いという問題があった。

この発明は、上述した従来技術による問題点を解消するため、簡単かつ信頼性の高い有限要素解析を実現することができる直交基底気泡関数要素数値解析方法、直交基底気泡関数要素数値解析プログラム、および直交基底気泡関数要素数値解析装置を提供することを目的とする。

上述した問題点を解決し、目的を達成するため、この発明の直交基底気泡関数要素数値解析方法、直交基底気泡関数要素数値解析プログラム、および直交基底気泡関数要素数値解析装置は、解析対象の各要素の整合質量行列を取得し、前記解析対象の各要素の気泡関数に基づいて、前記第２の取得工程によって取得された各要素の整合質量行列を対角化した、各要素の対角質量行列を生成し、解析対象の既知解析物理量と、生成された各要素の対角質量行列と、に基づいて、前記解析対象の挙動を解析することを特徴とする。

また、各要素における気泡関数の積分値を、各要素の整合質量行列に代入することにより、各要素の対角質量行列を生成することとしてもよい。さらに、生成された要素の対角質量行列に基づいて、前記解析対象全域の対角質量行列を算出し、前記解析対象全域の対角質量行列の逆行列を算出し、解析対象の既知解析物理量と、前記解析対象全域の対角質量行列と、前記逆行列と、に基づいて、前記解析対象の挙動を解析することとしてもよい。

これらの発明によれば、気泡関数要素を使用した有限要素解析において質量行列の近似を行わずに、質量行列が対角行列となる下記条件式（１４）を満たす気泡関数を用いる。

上記式（１５−１）〜式（１５−３）は、気泡関数要素の定式化をおこなうために必要となる積分であり、<・,・>Ω_eは要素領域Ωeの積分、Ａ_eは要素領域Ωeの面積（体積）を示している。これにより、各要素の整合質量行列内の成分を足し合わせることによって近似された集中質量行列を用いずに、気泡関数要素の基底（形状関数）が直交する条件を導入して、各要素の高精度な対角質量行列を用いて解析することができる。また、解析対象全域の対角質量行列やその逆行列も簡単に算出することができ、解析対象の挙動を、未知の解析対象物理量を用いて高精度に解析することができる。

本発明にかかる直交基底気泡関数要素数値解析方法、直交基底気泡関数要素数値解析プログラム、および直交基底気泡関数要素数値解析装置によれば、簡単かつ高精度に解析処理をおこなうことにより、解析対象の挙動解析の信頼性の向上を図ることができるという効果を奏する。また、記憶容量の低減化および解析時間の短縮化など、効率的な解析手法（直交基底気泡関数要素数値解析方法）を実現することができるという効果を奏する。

図１は、この発明の実施の形態にかかる数値解析装置のハードウェア構成を示すブロック図である。 図２は、この発明の実施の形態にかかる数値解析装置の機能構成を示すブロック図である。 図３は、この発明の実施の形態にかかる数値解析装置のデータ格納部が格納している解析物理量情報テーブルを示す説明図である。 図４は、この発明の実施の形態にかかる数値解析装置における数値解析処理手順を示すフローチャートである。 図５は、２次元の気泡関数要素を示す説明図である。 図６は、３次元の気泡関数要素を示す説明図である。 図７は、直交基底を形成する三角形気泡関数の形状φ_Bを示す説明図である。 図８は、直交基底を形成する三角形気泡関数の形状φ_B ²を示す説明図である。 図９は、直交基底を形成する三角形気泡関数の形状φ_Bを示す説明図である。 図１０は、直交基底を形成する三角形気泡関数の形状φ_B ²を示す説明図である。 図１１は、２次元の気泡関数を示す説明図である。 図１２は、３次元の気泡関数を示す説明図である。 図１３は、Rotating Cone問題の解析における計算領域を示す説明図である。 図１４は、Rotating Cone問題の解析における解析モデルのメッシュを示す説明図である。 図１５は、Rotating Cone問題の解析における解析モデルの流速（流況）を示す説明図である。 図１６は、Rotating Cone問題の解析における初期条件を示す鳥瞰図である。 図１７は、Rotating Cone問題の解析における初期条件の等高線を示す説明図である。 図１８は、気泡関数（式（４７））での整合質量行列を使用した計算結果を示す鳥瞰図である。 図１９は、気泡関数（式（４７））での整合質量行列を使用した計算結果の等高線を示す説明図である。 図２０は、気泡関数式（（４７））での集中質量行列を使用した計算結果を示す鳥瞰図である。 図２１は、気泡関数（式（４７））での集中質量行列を使用した計算結果の等高線を示す説明図である。 図２２は、気泡関数（式（４８））での整合質量行列を使用した計算結果を示す鳥瞰図である。 図２３は、気泡関数（式（４８）)での整合質量行列を使用した計算結果の等高線を示す説明図である。 図２４は、気泡関数（式（４８））での集中質量行列を使用した計算結果を示す鳥瞰図である。 図２５は、気泡関数（式（４８）)での集中質量行列を使用した計算結果の等高線を示す説明図である。 図２６は、気泡関数（式（４９））での対角質量行列を使用した計算結果を示す鳥瞰図である。 図２７は、気泡関数（式（４９））での対角質量行列を使用した計算結果の等高線を示す説明図である。 図２８は、Rotating Cone問題の解析における解析モデルのメッシュ（節点数：438413、要素数：2539200）を示す説明図である。 図２９は、Rotating Cone問題の解析における鉛直位置＝０の断面の初期条件を示す鳥瞰図である。 図３０は、Rotating Cone問題の解析における鉛直位置＝０の部分の初期条件の等高線を示す説明図である。 図３１は、気泡関数（式（５０））での整合質量行列を使用した鉛直位置＝０の部分の計算結果を示す鳥瞰図である。 図３２は、気泡関数（式（５０））での整合質量行列を使用した鉛直位置＝０の部分の計算結果の等高線を示す説明図である。 図３３は、気泡関数（式（５０））での集中質量行列を使用した鉛直位置＝０の部分の計算結果を示す鳥瞰図である。 図３４は、気泡関数（式（５０））での集中質量行列を使用した鉛直位置＝０の部分の計算結果の等高線を示す説明図である。 図３５は、気泡関数（式（５１））での整合質量行列を使用した鉛直位置＝０の部分の計算結果を示す鳥瞰図である。 図３６は、気泡関数（式（５１））での整合質量行列を使用した鉛直位置＝０の部分の計算結果の等高線を示す説明図である。 図３７は、気泡関数（式（５１））での集中質量行列を使用した鉛直位置＝０の部分の計算結果を示す鳥瞰図である。 図３８は、気泡関数（式（５１））での集中質量行列を使用した鉛直位置＝０の部分の計算結果の等高線を示す説明図である。 図３９は、気泡関数（式（５２））での対角質量行列を使用した鉛直位置＝０の部分の計算結果を示す鳥瞰図である。 図４０は、気泡関数（式（５２））での対角質量行列を使用した鉛直位置＝０の部分の計算結果の等高線を示す説明図である。 図４１は、アイソパラメトリック座標系ｒを示す説明図（その１）である。 図４２は、アイソパラメトリック座標系ｒを示す説明図（その２）である。 図４３は、式（８１）、（９３）を用いた三角形気泡関数の形状φ_Bを示す説明図である。 図４４は、式（８１）、（９３）を用いた三角形気泡関数の形状φ_B ²を示す説明図である。 図４５は、式（１３２）を用いた三角形３レベル気泡関数の形状を示す説明図である。 図４６は、式（８１）、（９３）を用いた三角形気泡関数と式（１３２）を用いた三角形３レベル気泡関数の積の形状を示す説明図である。 図４７は、従来の２次元の気泡関数要素を示す説明図である。 図４８は、従来の３次元の気泡関数要素を示す説明図である。 図４９は、従来の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で分割した要素群を全体系に重ね合せる作業によって形成され、要素領域での積分を用いて下記式（１６）のようにあらわすことができる。

式（１６）に対して、気泡関数要素の基底（形状関数）が直交するためには、下記式（１８）、式（２０）を満たす必要がある。なお、式（１８）の導入に際して、式（１７）を用いる。また、式（２０）の導入に際して、式（１９）を用いる。

式（１８）、（２０）より以下の下記関係式（２１）を得る。

＜２．直交条件を満たす気泡関数＞
式（１８）、式（２０）を満たすことができる気泡関数として下記式（２２）で表される気泡関数を提案する。

式（１８）に式（２２）を代入することにより下記式（２３）を得る。

式（２０）に式（２２）を代入することにより下記式（２４）を得る。

式（２３）に式（２４）を代入することにより、最終的にα₁を未知量とした下記式（２５）を得る。

式（２５）はα₁に関する二次方程式なので、解の公式より下記式（２６）にて未知量α₁を求めることができる。

＜３．直交基底を形成する具体的な気泡関数＞
直交基底を形成する具体的な気泡関数について説明する。図５は、２次元の気泡関数要素を示す説明図であり、図６は、３次元の気泡関数要素を示す説明図である。

直交基底を形成する具体的な気泡関数を導出するため、図５および図６に示した三角形(四面体)の要素領域を、重心点を用いて３（４）つの小三角形（小四面体）ｗ_i,ｉ＝１・・・Ｎ＋１に分割するξ乗気泡関数（０＜ξ＜∞）を用いる。ξ乗気泡関数は小三角形(小四面体)ごとにアイソパラメトリック座標系｛ｒ,ｓ｝（｛ｒ,ｓ,ｔ｝）を用いて次の式（２７）および式（２８）のように定義される。

これらの気泡関数の要素領域における積分値は、以下の式（２９）によって求めることができる。

式（２２）を、ξ乗気泡関数を用いて下記式（３０）のように定義する。

式（３０）の各べき乗の値を下記式（３１）のように設定すると二次元（三角形）、三次元（四面体）の気泡関数に対してα₁，α₂は下記式（３２）〜式（３５）のようになる。

式（３０）の各べき乗の値を下記式（３６）のように設定する（式（３１）とは異なった値）と二次元（三角形）、三次元（四面体）の気泡関数に対してα₁, α₂は下記式（３７）〜式（４０）のようになる。

図７は、式（３２）のα₁,α₂を使用した直交基底を形成する三角形気泡関数の形状φ_Bを示す説明図であり、図８は、式（３２）のα₁,α₂を使用した直交基底を形成する三角形気泡関数の形状φ_B ²を示す説明図である。また、図９は、式（３７）のα₁,α₂を使用した直交基底を形成する三角形気泡関数の形状φ_Bを示す説明図であり、図１０は、式（３７）のα₁,α₂を使用した直交基底を形成する三角形気泡関数の形状φ_B ²を示す説明図である。

式（３２）、式（３７）および図７〜図１０に示すように、式（２１）を満たす気泡関数（形状）は無数に存在するため、直交基底を形成する具体的な気泡関数やその形状には殆ど意味がなく、直交基底気泡関数要素の本質は気泡関数要素の基底が直交する条件式（２１）を発明した部分にある。

＜４．直交基底気泡関数要素の要素レベルの質量行列＞
下記式（４１）〜式（４６）に、２次元（三角形）、３次元（四面体）における通常の気泡関数要素の要素レベルの整合質量行列、集中質量行列、直交基底気泡関数要素の要素レベルの質量行列（対角質量行列）を示す。

式（４３）の要素レベルの対角質量行列Ｍ^(e) _ijは、生成部２０４により、式（４１）の要素レベルの整合質量行列Ｍ^(e) _ijの<１，φ_B>Ω_eおよび‖φ_B‖²Ω_eに、気泡関数φ_Bの積分値ＣＡ_eを代入することによって生成される。そして、第１の算出部２０５により、各要素の要素レベルの対角質量行列Ｍ^(e) _ijの総和（重ね合せ）をとることによって、解析対象全域の対角質量行列Ｍを算出することができる。また、第２の算出部２０６により、対角質量行列Ｍの逆行列Ｍ^-1を算出することができる。そして、解析部２０７により、既知解析物理量ｕ、対角質量行列Ｍおよび逆行列Ｍ^-1を、上述した式（１０）〜式（１３）に代入することによって、解析対象の挙動を解析することができる。

式（４６）の要素レベルの対角質量行列は、生成部２０４により、式（４４）の要素レベルの整合質量行列の<１，φ_B>Ω_eおよび‖φ_B‖²Ω_eに、気泡関数φ_Bの積分値ＣＡ_eを代入することによって生成される。そして、第１の算出部２０５により、各要素の要素レベルの対角質量行列Ｍ^(e) _ijの総和（重ね合せ）をとることによって、解析対象全域の対角質量行列Ｍを算出することができる。また、第２の算出部２０６により、対角質量行列Ｍの逆行列Ｍ^-1を算出することができる。そして、解析部２０７により、既知解析物理量ｕ、対角質量行列Ｍおよび逆行列Ｍ^-1を、上述した式（１０）〜式（１３）に代入することによって、解析対象の挙動を解析することができる。

＜５．気泡関数について＞
気泡関数について説明する。図１１は、２次元の気泡関数を示す説明図であり、図１２は、３次元の気泡関数を示す説明図である。気泡関数φ_Bは、形状関数の定義（重心点が１、その他の点では０となる関数）を満たせば任意の関数を選択することができる。検証の比較計算として、最も多く使用されている以下に示す二つの気泡関数および直交基底を形成する気泡関数を用いる。図１１および図１２に示したアイソパラメトリック座標系｛ｒ,ｓ｝（｛ｒ,ｓ,ｔ｝）を用いて、式（４７）〜式（５２）のように定義される。

＜６．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,Ψαは下記式（５３）〜（５５）になる。

図４２は、アイソパラメトリック座標系ｒを示す説明図（その２）である。図４２に示したアイソパラメトリック座標系ｒにおいて、中間点を用いて、線要素領域をｗ₁、ｗ₂に分割し、気泡関数の要素領域での積分値が式（２９）となる１次元のξ乗気泡関数は下記式（５６）のように定義される。

式（５６）、（２７）、（２８）のξ（η）乗気泡関数（０＜η＜∞）は、下記式（５７）に示すように、気泡関数を空間で偏微分した積分値も統一的に記述でき、下記式（５８）、（５９）のような積分値の関係を持っている。

式（５７）のｘ_k,ｘ_l＝ｘ₁・・・ｘ_Nはそれぞれ空間方向を表し、式（５９）のα_kdは実数を、Ｋ_Dは正の整数を示している。式（５９）は式（５８）より誘導している。Ｃ＝（Ｎ＋１）／（Ｎ＋２）とし、空間次元数Ｎを使用して直交基底気泡関数要素の条件式（１４）（または式（２１））を１次元〜３次元まで統一的に記述すると下記式（６０）になる。

下記式（６１）〜（６３）に、１次元（線）における通常の気泡関数要素の要素レベルの整合質量行列、集中質量行列、直交基底気泡関数要素の要素レベルの質量行列（対角質量行列）を示す。

＜９．拡張された直交基底気泡関数要素＞
気泡関数要素を使用した流体解析、構造解析、固有値解析などでは、一般に一次要素の形状関数Ψαの積分値に加えて、下記式（６４）〜（６６）の気泡関数の積分値が必要になる。

下記式（６７）と式（６４）より、下記の関係式（６８）を得る。

気泡関数として、下記式（６９）の関係を持っている気泡関数の場合には、必要となる気泡関数の積分値（６４）〜（６６）の式（６４）は、式（７０）に代用することができる。

式（６９）の関係式は、非特許文献３で現れるような特殊な気泡関数を導入しない限り、例えば式（４７）、（４８）、（５０）、（５１）などの一般的によく使用されている気泡関数で成り立つので、直交条件式（２１）の導出仮定で現れる式（１７）、（１９）についても使用しており、式（２２）は式（６９）を満たす気泡関数であると仮定している。式（２２）右辺の気泡関数として式（３０）を採用した場合には、式（５９）より式（２２）は式（６９）を満足していることがわかる。

直交基底気泡関数要素を使用した場合には、式（７０）、（６５）、（６６）において、式（７０）、（６５）は直交条件式（６０）よりわかるが、式（６６）は、直交条件により決定することができない。式（５６）、（２７）、（２８）のξ乗気泡関数を使用して、直交基底気泡関数を形成した場合、実際に、式（６６）は、下記式（７１）で表される。

式（７１）をみると、Ｄ（ξ₁,ξ₂,ξ₃）の値が、変数ξ₁、ξ₂、ξ₃の関数となっている。実際に、式（３１）、（３６）の値を用いて、式（７２）を計算（α₁（ξ₁,ξ₂,ξ₃）＝α₁₊（ξ₁,ξ₂,ξ₃）を使用）してみると、下記式（７３）〜（７５）のように採用する変数によって計算される値が異なる。

本発明では、直交基底気泡関数要素の拡張として、直交条件式（６０）によって得られる積分値を利用して式（６６）の積分値を決定できる、下記式（７６）を新たに定義する。

基底が直交する条件式（６０）に、式（７６）の条件を追加した下記式（８０）を拡張された直交基底気泡関数要素の条件式とする。

＜１０．拡張された直交基底気泡関数要素を満たす気泡関数＞
式（７１）からわかるように、式（５６）、（２７）、（２８）のξ乗気泡関数を使用して直交基底気泡関数を形成した式（６６）の積分値は、変数ξ₁、ξ₂、ξ₃の値によって決まる。本特許では、拡張された直交基底気泡関数要素の条件式（８０）を満足するような変数ξ₁、ξ₂、ξ₃を考え、そのような変数が実際に存在することを示す。変数はξ₁、ξ₂、ξ₃の三つあるが、その内の２変数を下記式（８１）のように固定し、ξ₁の１変数の問題として、拡張された直交基底気泡関数要素の条件を満足する気泡関数を求めることを考える。

ξ₁の１変数の問題として、下記式（８２）を定義する。

拡張された直交基底気泡関数要素の条件を満たす気泡関数が存在することは、実数の集合Ｒにおいて、ｆ（ξ₁）＝０となる変数が存在することである。いま、下記式（８３）〜（８５）のような実数の集合Ｒの部分集合Ｓを考える。

ξ₁∈Ｓでは、１次元〜３次元において、下記式（８６）となる。

従って、α₁、α₂は実数の範囲で値を持つ。また、ξ₁∈Ｓにて１次元〜３次元では、それぞれ下記式（８７）および式（８８）が成り立つので、式（７２）は有界な値をとり、式（８２）は実変数の関数となる。

式（８２）と部分集合Ｓを使用すると下記式（８９）〜（９１）を得ることが出来る。

任意のξ_1c∈Ｓに対して、１次元では、式（８９）の第１式が負、第２式が正で、第３式が正（単調増加）、２次元、３次元では、式（９０）、（９１）の第１式が正、第２式が負で、第３式が負（単調減少）であるので、１次元〜３次元において、部分集合Ｓの範囲には、それぞれｆ（ξ₁）＝０の解が一意に存在する。

具体的に、部分集合Ｓの範囲内で、二分法を使用してｆ（ξ₁）＝０の解を求めると、下記式（９２）〜（９４）の値が得られる。

図４３、図４４に、式（８１）、（９３）を用いた三角形気泡関数の形状φ_B、三角形気泡関数の形状φ_B ²を示す。

または、下記式（９５）、（９６）を使用して、ｆ（ξ₁）＝０の解の存在を示すことができる。

１次元〜３次元において、ｈを正の定数とし、任意のξ_1c∈Ｓに対して下記式（９７）の結果を得る。

平均値の定理から、任意のξ_1a, ξ_1b∈Ｓに対して、下記式（９８）が得られる。

従って、縮小写像の原理（定理）より、式（９５）を反復的に逐次計算をさせる下記式（９９）にて得られる数列｛ξ₁ ^m｝は、部分集合Ｓの範囲内でｆ（ξ₁）＝０を満足する一意の解に収束する。

なお、拡張された直交基底気泡関数要素の存在を示すために、式（８０）の三つの条件を満たす気泡関数として、三つの未知量をα₁、α₂、ξ₁とし、適当な気泡関数φ_B1、φ_B2、φ_B3を採用しが、他にも、下記式（１００）の気泡関数φ_Bのように、未知量をα₁、α₂、α₃として、適当な気泡関数φ_B1、φ_B2、φ_B3、φ_B4を採用すれば、式（８０）の条件を満たす気泡関数が得られる。このように、拡張された直交基底気泡関数要素の条件式を満たす気泡関数は、唯一ではなく、幾つも存在するため、具体的な気泡関数やその形状には殆ど意味がなく、拡張された直交基底気泡関数要素で重要な意味を成すのは、式（８０）の条件となる気泡関数の積分値である。

＜１１．拡張された直交基底気泡関数要素の要素レベルの粘性項（拡散項）の行列＞
通常の気泡関数、直交基底気泡関数要素の数値解析において式（６６）が必要になる要素レベルの粘性項（拡散項）の行列を下記式（１０１）〜（１０３）に示す。

式（１０１）〜（１０３）に式（７６）を代入して整理した、拡張された直交基底気泡関数要素の要素レベルの粘性項（拡散項）の行列を下記式（１０４）〜（１０６）に示す。

以上、拡張された直交基底気泡関数要素の要素レベルの粘性項（拡散項）の行列（１０４）〜（１０６）を示したが、本発明に関する定式化やプログラムでは、行列（１０１）〜（１０３）に式（７６）を代入して整理を行わずに、そのまま行列（１０１）〜（１０３）と式（７６）を使用することも、もちろん可能である。

＜１２．直交基底気泡関数要素と正規化気泡関数＞
式（１）は、下記式（１０７）、（１０８）のように変形することができる。

式（１０８）は、気泡関数自体の積分値がＡ_e（１次元では各要素ｅの線、２次元では各要素ｅの面積、３次元では各要素ｅの体積）となるような気泡関数を使用して、式変形を行ったものであり、式（１１０）は正規化気泡関数と呼ばれている。直交基底気泡関数要素において式（１）の表現形式の他に、式（１０７）の表現形式や式（１０８）の表現形式を用いて、定式化やプログラミング、計算実行を行うことも可能である。また、式（１）、（１０７）、（１０８）の表現形式を用いて得られた有限要素方程式（変分方程式）では、静的縮約という操作によって要素ごとに重心点の自由度ｕ_B、ｕ_B’、ｕ_S’を消去することができるが、重心点の自由度を消去した有限要素方程式（変分方程式）を用いて、定式化やプログラミング、計算実行を行うことも可能である。直交基底気泡関数要素の気泡関数として正規化気泡関数を使用した積分値は下記式（１１１）、（１１２）のようになる。

さらに、拡張された直交基底気泡関数要素では、下記式（１１３）が得られる。

＜１３．直交基底気泡関数要素と安定化作用制御項＞
気泡関数要素の持つ数値不安定性を緩和させるために、下記式（１１４）の気泡関数の粘性項（安定化作用制御項）を用いて解析精度を向上させる方法がある。

式（１１４）の安定化作用制御パラメータν’（−∞＜ν’＜∞）は、粘性係数の次元（粘性係数に相当する変数を無次元化している場合は無次元）を持つパラメータであり、この値を適切に決定することによって高精度解法を実現する。この方法は、従来、気泡関数要素に通常の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、非特許文献３参照）。３レベル気泡関数は、下記式（１１５）〜（１１７）を満たす気泡関数である。

式（６７）、（１１５）より下記式（１１８）を得る。

また、下記式（１１９）の関係を持っている３レベル気泡関数の場合には、条件式（１１５）〜（１１７）の式（１１５）は式（１２０）で代用できる。

式（１１９）の関係を持っている３レベル気泡関数を仮定し、上記の条件式（１２０）、（１１６）、（１１７）を満たす具体的な３レベル気泡関数として、下記式（１２１）を用いる。

式中のνは、粘性係数（拡散係数）である。式（１２１）を式（１２０）、（１１６）、（１１７）に代入して整理すると、下記式（１２２）を得る。

式（１２１）右辺の四つの気泡関数を具体的に下記式（１２３）のように定義する。

式（１２３）、（５９）より式（１２１）は式（１１９）を満足していることがわかる。式（１２３）の変数を下記式（１２４）のように設定する。

α₁（ξ₁,ξ₂,ξ₃）＝α₁₊（ξ₁,ξ₂,ξ₃）として、式（３１）の変数を用いると、式（１２２）より下記式（１２５）〜（１２７）を得る。

α₁（ξ₁,ξ₂,ξ₃）＝α₁₊（ξ₁,ξ₂,ξ₃）として、式（３６）の変数を用いると、式（１２２）より下記式（１２８）〜（１３０）を得る。

α₁（ξ₁,ξ₂,ξ₃）＝α₁₊（ξ₁,ξ₂,ξ₃）として、式（９２）〜（９４）、式（８１）の変数を用いると、式（１２２）より下記式（１３１）〜（１３３）を得る。

図４５、図４６に、式（１３２）を用いた三角形３レベル気泡関数の形状（ν’＝νとした場合）、式（９３）、（８１）を用いた三角形気泡関数と式（１３２）を用いた三角形３レベル気泡関数の積の形状（ν’＝νとした場合）を示す。以上により、直交基底気泡関数要素を用いた場合にも、式（１２０）（あるいは式（１１５））、（１１６）、（１１７）を満たす３レベル気泡関数は存在し、安定化作用制御項（１１４）を導入できる。

なお、３レベル気泡関数の存在を示すため、式（１２０）、（１１６）、（１１７）の三つの条件を満足する三つの未知量δ₁、δ₂、δ₃と四つの適当な気泡関数を持った３レベル気泡関数（１２１）を採用したが、この他にも、下記式（１３４）のように、安定化作用制御パラメータν’と同じ次元を持つパラメータを導入し、未知量をδ₁、δ₂として、式（１３４）の右辺に適当な三つの気泡関数を用いることによって、式（１２０）、（１１６）の条件を満たす３レベル気泡関数を求めることができる。そして、下記式（１３５）のように安定化作用制御パラメータν’を定義すれば、式（１１７）も導出でき、式（１２０）（あるいは式（１１５））、（１１６）、（１１７）となる３レベル気泡関数の存在を示すことが出来る。このように、３レベル気泡関数の条件式を満たす気泡関数は、唯一ではなく、幾つも存在するため、具体的な気泡関数やその形状には殆ど意味がなく、３レベル気泡関数で重要な意味を成すのは、安定化作用制御項を導くために必要な式（１２０）（あるいは式（１１５））、（１１６）、（１１７）の条件である。

＜１４．直交基底気泡関数要素と正規化気泡関数による安定化作用制御項＞
正規化気泡関数を使用した場合にも安定化作用制御項は導入することが出来、下記式（１３６）によって表される。

直交基底気泡関数要素で使用する気泡関数として正規化気泡関数を用いた場合、３レベル正規化気泡関数は下記式（１３７）と定めることができ、その積分値は下記式（１３８）〜（１４２）のように記述できる。

上記式（１３８）〜（１４２）より、式（１１５）〜（１１７）、（１１９）、（１２０）を満足する３レベル気泡関数が存在すれば、式（１３８）〜（１４２）も自動的に満たすので３レベル正規化気泡関数も同時に存在することがわかる。

＜１５．直交基底気泡関数要素と気泡関数の安定化作用行列＞
式（１０７）の表現形式を用いた気泡関数要素の重心点の自由度ｕ_B’は、静的縮約という操作によって要素ごとに消去できる。その結果得られる有限要素方程式は、安定化有限要素法で導かれる有限要素方程式と密接な関係がある（非特許文献２参照）。静的縮約によって得られる気泡関数の安定化作用行列は、一般的に下記式（１４３）の形式によって表される。

Ｌは、求めるべき未知解析物理量の数である。式（１４４）は、式（６９）を満たす気泡関数を使用した場合、式（７０）、（６５）、（６６）の気泡関数の積分値から成る行列である。安定化作用制御パラメータν’は、下記式（１４５）の関係が満たされるように決定する。

式（１４６）は解析精度が向上することが期待できる安定化行列である。式（１４５）の関係を満たすような安定化作用制御パラメータν’を求めるために、下記式（１４８）の評価関数を最小にするような安定化作用制御パラメータν’を決定する問題を考える。

式（１４８）を最小にするような安定化作用制御パラメータν’を求めるために、式（１４８）の停留条件である下記関係式（１５０）を用いて、安定化作用制御パラメータν’を求める。

３レベル気泡関数は重み関数の中で導入しているので、求めるべき未知解析物理量の偏微分方程式ごとに異なった３レベル気泡関数を使用することができる。従って、安定化作用制御パラメータν’は、下記式（１５１）の３レベル気泡関数を用いると、求めるべき未知解析物理量の数Ｌだけ、安定化作用制御パラメータν_m’を定義することができる。

式（１５１）の係数δ₁、δ₂、δ₃は、式（１２２）によって決定することが出来、式（１５１）は式（１２０）、（１１６）を満たす。また、式（１１７）の代わりに下記式（１５２）の安定化作用制御項を得る。

式（１５２）の安定化作用制御項を用いた場合、式（１４４）は下記式（１５３）のようになる。

式（１４５）の関係を満たすような安定化作用制御パラメータν_m’を求めるために、下記式（１５４）の評価関数を最小にするような安定化作用制御パラメータν_m’を決定する問題を考える。

式（１５４）を最小にするような安定化作用制御パラメータν_m’を求めるために、式（１５４）の停留条件である下記関係式（１５６）を用いて、安定化作用制御パラメータν_m’を求める。

式（１５５）はベクトルであり、ベクトルの成分を下記式（１５７）のように表現する。

式（１５７）をみるとわかるように、安定化作用制御パラメータν_m’は各ベクトル成分で独立に定義されているので、式（１４５）の関係を満たすような安定化作用制御パラメータν_m’を決定する問題として、式（１５７）の成分ｍごとに、下記式（１５８）の評価関数を用いた最小化問題を考えることができる。

式（１５８）を最小にするような安定化作用制御パラメータν_m’を求めるために、式（１５８）の停留条件である下記関係式（１５９）を用いて、安定化作用制御パラメータν_m’を求める。

式（１５７）は、安定化作用制御パラメータν_m’が各ベクトル成分で独立に定義されているので、式（１５４）、（１５８）の停留条件から導かれる安定化作用制御パラメータν_m’を求める式（１５６）、（１５９）は、結果的に同じ式になる。

他の方法として、下記式（１６０）の関係を満たすように、安定化作用制御パラメータν’を決定する問題も考えられる。

式（１６０）の関係を満たすような安定化作用制御パラメータν’を求めるために、下記式（１６１）の評価関数を最小にするような安定化作用制御パラメータν’を決定する問題を考える。

式（１６１）を最小にするような安定化作用制御パラメータν’を求めるために、式（１６１）の停留条件である下記関係式（１６３）を用いて、安定化作用制御パラメータν’を求める。

安定化作用制御パラメータν’は、式（１５１）の３レベル気泡関数を用いると、式（１５２）の安定化作用制御項が導出され、求めるべき未知解析物理量の数Ｌだけ、安定化作用制御パラメータν_m’を定義することができる。式（１６０）の関係を満たすような安定化作用制御パラメータν_m’を求めるために、下記式（１６４）の評価関数を最小にするような安定化作用制御パラメータν_m’を決定する問題を考える。

式（１６４）を最小にするような安定化作用制御パラメータν_m’を求めるために、式（１６４）の停留条件である下記関係式（１６６）を用いて、安定化作用制御パラメータν_m’を求める。

式（１６５）はベクトルであり、ベクトルの成分を下記式（１６７）のように表現する。

式（１６７）をみるとわかるように、安定化作用制御パラメータν_m’は各ベクトル成分で独立に定義されているので、式（１６０）の関係を満たすような安定化作用制御パラメータν_m’を決定する問題として、式（１６７）の成分ｍごとに、下記式（１６８）の評価関数を用いた最小化問題を考えることができる。

式（１６８）を最小にするような安定化作用制御パラメータν_m’を求めるために、式（１６８）の停留条件である下記関係式（１６９）を用いて、安定化作用制御パラメータν_m’を求める。

式（１６７）は、安定化作用制御パラメータν_m’が各ベクトル成分で独立に定義されているので、式（１６４）、（１６８）の停留条件から導かれる安定化作用制御パラメータν_m’を求める式（１６６）、（１６９）は、結果的に同じ式になる。

＜１６．直交基底気泡関数要素と正規化気泡関数の安定化作用行列＞
式（１０８）の表現形式を用いた気泡関数要素の重心点の自由度ｕ_S’は、静的縮約という操作によって要素ごとに消去できる。静的縮約によって得られる正規化気泡関数の安定化作用行列は、一般的に下記式（１７０）の形式によって表される。

式（１７０）の第二項は、式（１４１）の関係を満たす気泡関数を使用した場合、式（１１１）〜（１１３）各第一項の正規化気泡関数の積分値から成る行列であり、式（７０）、（６５）、（６６）で表した場合には第三項のように記述できる。安定化作用制御パラメータν’は、下記式（１７２）の関係が満たされるように決定する。

式（１７２）の関係を満たすような安定化作用制御パラメータν’を求めるために、下記式（１７４）の評価関数を最小にするような安定化作用制御パラメータν’を決定する問題を考える。

式（１７４）を最小にするような安定化作用制御パラメータν’を求めるために、式（１７４）の停留条件である下記関係式（１７６）を用いて、安定化作用制御パラメータν’を求める。

３レベル正規化気泡関数は重み関数の中で導入しているので、求めるべき未知解析物理量の偏微分方程式ごとに異なった３レベル正規化気泡関数を使用することができる。従って、安定化作用制御パラメータν’は、下記式（１７７）の３レベル正規化気泡関数を用いると、求めるべき未知解析物理量の数Ｌだけ、安定化作用制御パラメータν_m’を定義することができる。

式（１７７）の係数δ₁、δ₂、δ₃は、式（１２２）によって決定することが出来、式（１７７）は式（１４２）、（１３９）を満たす。また、式（１４０）の代わりに下記式（１７８）の安定化作用制御項を得る。

式（１７８）の安定化作用制御項を用いた場合、式（１７１）は下記式（１７９）のようになる。

式（１７９）の第一項は、式（１１１）〜（１１３）各第一項の正規化気泡関数の積分値から成る行列であり、式（７０）、（６５）（６６）で表した場合には第二項のように記述できる。式（１７２）の関係を満たすような安定化作用制御パラメータν_m’を求めるために、下記式（１８０）の評価関数を最小にするような安定化作用制御パラメータν_m’を決定する問題を考える。

式（１８０）を最小にするような安定化作用制御パラメータν_m’を求めるために、式（１８０）の停留条件である下記関係式（１８２）を用いて、安定化作用制御パラメータν_m’を求める。

式（１８１）はベクトルであり、ベクトルの成分を下記式（１８３）のように表現する。

式（１８３）をみるとわかるように、安定化作用制御パラメータν_m’は各ベクトル成分で独立に定義されているので、式（１７２）の関係を満たすような安定化作用制御パラメータν_m’を決定する問題として、式（１８３）の成分ｍごとに、下記式（１８４）の評価関数を用いて最小化問題を考えることができる。

式（１８４）を最小にするような安定化作用制御パラメータν_m’を求めるために、式（１８４）の停留条件である下記関係式（１８５）を用いて、安定化作用制御パラメータν_m’を求める。

式（１８３）は、安定化作用制御パラメータν_m’が各ベクトル成分で独立に定義されているので、式（１８０）、（１８４）の停留条件から導かれる安定化作用制御パラメータν_m’を求める式（１８２）、（１８５）は、結果的に同じ式になる。

他の方法として、下記式（１８６）の関係を満たすように、安定化作用制御パラメータν’を決定する問題も考えられる。

式（１８６）の関係を満たすような安定化作用制御パラメータν’を求めるために、下記式（１８７）の評価関数を最小にするような安定化作用制御パラメータν’を決定する問題を考える。

式（１８７）を最小にするような安定化作用制御パラメータν’を求めるために、式（１８７）の停留条件である下記関係式（１８９）を用いて、安定化作用制御パラメータν’を求める。

安定化作用制御パラメータν’は、式（１７７）の３レベル正規化気泡関数を用いると、式（１７８）の安定化作用制御項が導出され、求めるべき未知解析物理量の数Ｌだけ、安定化作用制御パラメータν_m’を定義することができる。式（１８６）の関係を満たすような安定化作用制御パラメータν_m’を求めるために、下記式（１９０）の評価関数を最小にするような安定化作用制御パラメータν_m’を決定する問題を考える。

式（１９０）を最小にするような安定化作用制御パラメータν_m’を求めるために、式（１９０）の停留条件である下記関係式（１９２）を用いて、安定化作用制御パラメータν_m’を求める。

式（１９１）はベクトルであり、ベクトルの成分を下記式（１９３）のように表現する。

式（１９３）をみるとわかるように、安定化作用制御パラメータν_m’は各ベクトル成分で独立に定義されているので、式（１８６）の関係を満たすような安定化作用制御パラメータν_m’を決定する問題として、式（１９３）の成分ｍごとに、下記式（１９４）の評価関数を用いた最小化問題を考えることができる。

式（１９４）を最小にするような安定化作用制御パラメータν_m’を求めるために、式（１９４）の停留条件である下記関係式（１９５）を用いて、安定化作用制御パラメータν_m’を求める。

式（１９３）は、安定化作用制御パラメータν_m’が各ベクトル成分で独立に定義されているので、式（１９０）、（１９４）の停留条件から導かれる安定化作用制御パラメータν_m’を求める式（１９２）、（１９５）は、結果的に同じ式になる。

＜１７．気泡関数と正規化気泡関数の安定化作用制御パラメータの関係＞
直交基底気泡関数要素で、式（１）、（１０７）の表現形式を採用した場合には、式（７０）、（６５）、（６６）の積分を使用して、式（１５０）、（１５６）、（１５９）、（１６３）、（１６６）、（１６９）のいずれかを採用し、式（１０８）の表現形式を採用した場合には、式（１１１）〜（１１３）各第一項の積分を使用して、式（１７６）、（１８２）、（１８５）、（１８９）、（１９２）、（１９５）のいずれかを採用すれば安定化作用制御パラメータを計算することが出来る。いま、式（１７６）、（１８２）、（１８５）、（１８９）、（１９２）、（１９５）を、式（７０）、（６５）、（６６）の積分を使用して記述すると下記式（１９６）〜（２０１）のようになる。

式（７０）、（６５）、（６６）の積分を使用した式（１９６）〜（２０１）の第二項より、安定化作用制御パラメータを求めるための関係式は、式（１５０）、（１５６）、（１５９）、（１６３）、（１６６）、（１６９）の安定化作用制御パラメータを求めるための関係式と同じに（等しく）なる。従って、式（１５０）、（１５６）、（１５９）、（１６３）、（１６６）、（１６９）より式（７０）、（６５）、（６６）の気泡関数の積分を使用して求めた安定化作用制御パラメータと、式（１７６）、（１８２）、（１８５）、（１８９）、（１９２）、（１９５）より式（１１１）〜（１１３）各第一項の正規化気泡関数の積分を使用して求めた安定化作用制御パラメータは、同じ値（等しい値）になることがわかる。

上記より、１）気泡関数、気泡関数自体の積分値がＡ_eとなる正規化気泡関数、２）気泡関数の持つ数値不安定を緩和させるために導入する安定化作用制御項、および３レベル気泡関数、３レベル正規化気泡関数、３）安定化作用制御パラメータの計算で用いる安定化作用行列、の定式化やプログラミングにおいて、式（７０）、（６５）、（６６）の気泡関数の積分値、あるいは式（１１１）〜（１１３）各第一項の正規化気泡関数の積分値が必要になる。拡張された直交基底気泡関数要素は、必要となる積分値が全て事前に（複雑な積分を行うことなく）式（８０）によって与えられているので、整合質量行列と同等の解析精度を持つ対角質量行列の作成や気泡関数の持つ数値不安定性を改良する安定化作用制御項の導入で現れる一連の作業において、非常に合理的に処理を行うことが出来る。

直交基底気泡関数要素は、２次元、３次元の解析結果からわかるように、記憶容量、計算時間の効率の良い解法（例えば、４段解法）を採用しながら、信頼性のある結果を得ることができる。さらに、拡張された直交基底気泡関数要素の条件式（８０）により、インクジェットの液滴挙動解析、都市や河川の氾濫解析、地震等で発生する津波の解析のみならず、構造物の応力解析、自動車のエンジン固有振動解析などで必要となる気泡関数要素の積分値や気泡関数の持つ数値不安定性を改良するために必要となる処理を合理的、かつ統一的に扱うことができるので、上記のような多岐にわたる数値解析への適用が容易に行えるという効果をもたらす。

以上説明したように、本実施の形態によれば、気泡関数要素を使用した数値解析において、記憶容量、計算時間などの計算効率の良い解析手法（直交基底気泡関数要素数値解析方法）を実現することができる。

なお、本実施の形態で説明した直交基底気泡関数要素数値解析方法は、予め用意されたプログラムをパーソナル・コンピュータ、ワークステーション、ＰＣクラスタ、スーパーコンピュータ等のコンピュータで実行することにより実現できる。このプログラムは、ハードディスク、フレキシブルディスク、ＣＤ−ＲＯＭ、ＭＯ、ＤＶＤ等のコンピュータで読み取り可能な記録媒体に記録され、コンピュータにて記録媒体から読み出されることによって実行される。またこのプログラムは、インターネット等のネットワークを介して配布することが可能な伝送媒体であってもよい。

以上のように、本発明にかかる直交基底気泡関数要素数値解析方法、直交基底気泡関数要素数値解析プログラムおよび直交基底気泡関数要素数値解析装置は、線、三角形および四面体形状のメッシュを使用した有限要素法に基づく数値解析に有用である。

Claims (9)

1. 生成手段と解析手段とを備えたコンピュータが、
   前記生成手段が、解析対象範囲に分割形成されたメッシュごとに、当該メッシュの解析対象となる有限要素方程式内の整合質量行列に対して、気泡関数要素の基底が直交する条件を満たす気泡関数を求め、前記整合質量行列に、前記気泡関数の積分値を代入することによって、前記整合質量行列を対角化した対角質量行列を生成する生成工程と、
   前記解析手段が、前記解析対象となる物理量の物性を表す既知物性値および前記解析対象となる物理量の解析条件を表す境界解析物理量からなる既知解析物理量を、前記生成工程によって生成されたメッシュごとの対角質量行列に代入することによって、前記解析対象となる物理量の変化を解析する解析工程と、
   を実行することを特徴とする直交基底気泡関数要素数値解析方法。
2. 前記解析手段が、前記既知物性値、前記境界解析物理量および前記解析対象となる物理量の解析開始時の値を表す初期解析物理量からなる既知解析物理量を、前記生成工程によって生成されたメッシュごとに対角質量行列に代入することによって、前記解析対象となる物理量の微少時間における変化を解析することを特徴とする請求項１に記載の直交基底気泡関数要素数値解析方法。
3. 前記生成工程は、
   前記整合質量行列に対して、気泡関数要素の基底が直交する下記式（１）を満たす気泡関数を求め、前記整合質量行列に、前記気泡関数の積分値を代入することによって、前記整合質量行列を対角化した対角質量行列を生成することを特徴とする請求項１または２に記載の直交基底気泡関数要素数値解析方法。
   

   
4. 前記コンピュータは、さらに、第１の算出手段と第２の算出手段とを備え、
   前記第１の算出手段が、前記生成工程によって生成された前記解析対象範囲のメッシュごとの対角質量行列の総和を求めることによって、前記解析対象範囲全域の対角質量行列を算出する第１の算出工程と、
   前記第２の算出手段が、前記第１の算出工程によって算出された前記解析対象範囲全域の対角質量行列の逆行列を算出する第２の算出工程と、を実行し、
   前記解析工程は、
   前記解析対象範囲全域についての前記既知解析物理量と前記第２の算出工程によって算出された前記対角質量行列の逆行列とを乗算することによって、前記解析対象となる物理量の変化を解析することを特徴とする請求項１または３に記載の直交基底気泡関数要素数値解析方法。
5. 請求項１〜４のいずれか一つに記載の直交基底気泡関数要素数値解析方法を前記コンピュータに実行させることを特徴とする直交基底気泡関数要素数値解析プログラム。
6. 解析対象範囲に分割形成されたメッシュごとに、当該メッシュの解析対象となる有限要素方程式内の整合質量行列に対して、気泡関数要素の基底が直交する条件を満たす気泡関数を求め、前記整合質量行列に、前記気泡関数の積分値を代入することによって、前記整合質量行列を対角化した対角質量行列を生成する生成手段と、
   前記解析対象となる物理量の物性を表す既知物性値および前記解析対象となる物理量の解析条件を表す境界解析物理量からなる既知解析物理量を、前記生成工程によって生成されたメッシュごとの対角質量行列に代入することによって、前記解析対象となる物理量の変化を解析する解析手段と、
   を備えることを特徴とする直交基底気泡関数要素数値解析装置。
7. 前記解析手段は、
   前記既知物性値、前記境界解析物理量および前記解析対象となる物理量の解析開始時の値を表す初期解析物理量からなる既知解析物理量を、前記生成工程によって生成されたメッシュごとに対角質量行列に代入することによって、前記解析対象となる物理量の微少時間における変化を解析することを特徴とする請求項６に記載の直交基底気泡関数要素数値解析装置。
8. 前記生成手段は、
   前記整合質量行列に対して、気泡関数要素の基底が直交する下記式（２）を満たす気泡関数を求め、前記整合質量行列に、前記気泡関数の積分値を代入することによって、前記整合質量行列を対角化した対角質量行列を生成することを特徴とする請求項６または７に記載の直交基底気泡関数要素数値解析装置。
   

   
9. 前記生成手段によって生成された前記解析対象範囲のメッシュごとの対角質量行列の総和を求めることによって、前記解析対象範囲全域の対角質量行列を算出する第１の算出手段と、
   前記第１の算出手段によって算出された前記解析対象範囲全域の対角質量行列の逆行列を算出する第２の算出手段と、を備え、
   前記解析手段は、
   前記解析対象範囲全域についての前記既知解析物理量と前記第２の算出手段によって算出された前記対角質量行列の逆行列とを乗算することによって、前記解析対象となる物理量の変化を解析することを特徴とする請求項６または８に記載の直交基底気泡関数要素数値解析装置。

Images