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

Abstract

　まず、第１の取得部（２０２）により、解析対象の既知解析物理量を取得する（Ｓ４０１）。つぎに、第２の取得部（２０３）により、各要素の要素レベルの整合質量行列を取得する（Ｓ４０２）。そして、要素ごとに、気泡関数を積分し（Ｓ４０３）、各要素の要素レベルの整合質量行列に、ステップＳ４０３で積分された値を代入することによって、各要素の要素レベルの対角質量行列を算出する（Ｓ４０４）。つぎに、各要素の要素レベルの対角質量行列の総和（重ね合せ）により、解析対象全域の対角質量行列を算出する（Ｓ４０５）。そして、この対角質量行列の逆行列を算出し（Ｓ４０６）、解析対象の既知解析物理量と、解析対象全域の対角質量行列と、その逆行列とに基づいて、解析対象の挙動を解析する（Ｓ４０７）。

Description

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

技術分野

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

背景技術

[0002] 従来の気泡関数要素について説明する。図 47は、従来の 2次元の気泡関数要素 を示す説明図であり、図 48は、従来の 3次元の気泡関数要素を示す説明図である。 図 47および図 48のように、三角形（四面体)要素を用いた気泡関数要素は、各要素 において三角形（四面体)を形成する 3 (4)点と重心点の 4 (5)つの節点を用いて、ァ イソパラメトリック座標系 [r,_S] ({r,s,t})で下記式（1)のように表される（たとえば、下記 非特許文献 1、非特許文献 2、非特許文献 3参照)。

[0003] [数 1]

N+1

«=ι · · ·式 （ 1 )

*α = *α - 7 ~； < , α = 1 - - - Ν + 1

W · · ·式 （2 )

[0004] 式（1)の Φ α , φ は、気泡関数要素の形状関数、 u a ,uは三角形（四面体)の各頂

B B

点の値 (解析物理量）、重心点の値 (解析物理量）、 Nは空間次元数を示している。 ベクトル形式で記述すれば、形状関数は下記式（3)〜式 (6)になる。

[0005] [数 2] Φ^Γ = [Φχ _{2 3} ΦΒ] · · .式 （ ₃ )

«^Τ = [ui u₂ u₃ U_B] . . .式 ）

3次元

Φ^Γ = [ ι Φ₂ *3 ^ Φβ] · . .式 （ ₅ )

U^T = [Ui U₂ U_S U4 u_B] . ■ .式 （ ₆ )

[0006] 式（2)の Ψ aは、 2次元および 3次元の一次要素を用いた形状関数であり、下記式

(7)、（8)であらわされる。

[0007] [数 3]

2次元

Φΐ = 1 - Γ— S, ^2 = Τ, $3 = « , · . ¹■ 7 )

3次元

Φι = 1 - r - s - 1, Φ₂ = r, *3 = -s, Φ₄ = ί . . .式 （ 8 )

[0008] 形状関数 φ

Βは気泡関数と呼ばれている。気泡関数は要素境界上においてその値 力^となり、重心点で値が 1となるように要素毎に定義される。非定常問題において、 空間方向の離散化に気泡関数要素を用いた有限要素方程式は下記式 (9)のように あらわすことができる。

[0009] [数 4]

M u + F(u) = . . . （_{9 )}

[0010] 式（9)の uは求めるべき未知解析物理量 (汚染物質濃度、温度、流量、水深、流速 、圧力、変位など)であり、 Mは質量行列、 F (u)は時間微分項以外をまとめた項であ る。式 (9)の時間方向の離散化として、ティラー展開に基づいた 4段解法は下記式（ 10)〜式（13)のように表される（たとえば、下記非特許文献 4参照)。

[0011] [数 5] < 1st step > 式 （1 0 )

< 2 式 （1 1 )

< 3

式 ( 1 2 ) く 4th step >

uⁿ⁺¹ = u ^_1Δί (uⁿ⁺³/⁴) 式 （1 3 )

[0012] 式（10)〜式（13)の上付き添え字 nは現在時刻 nでの既知解析物理量を表し、 n+ 1は時刻 nから微小時間 Δ t経過後の未知解析物理量を表して ヽる。

[0013] 非特許文献 1： D.N.Arnold, F.Brezzi and M.Fortin, "A Stable Finite Element for the Stokes Equations", Calcolo, Vol.23, 1984, pp.337— pp.344 (ディ^ ~ ·ェヌ 'アーノルド 、エフ.ブリジィ、ェム.フォーティン著 「ァ スティブル フアイナイト エレメント フォ 一 ザ スト一タス イクエイシヨンズ」 カルコ口 23卷 1984年 337頁 344頁） 特干文献 2 :J.し.; mo, F.Armero andし. A.Tayior, Stable and Time— Dissipative Fi nite Element Methods for thelncompressiole Navier— Stokes Equations in Advection Dominated Flows", International Journal for Numerical Methods in Engineering'Vol. 38, 1995, pp.1475- pp.1506 (ジェ^ ~ ·シ' ~ ·シモ、エフ'アルメロ、シ' ~ ·エ^ ~ ·ティラー 著 「スティブル アンド タイムーディシペイティブ フアイナイト エレメント メソーズ フォー ザ インコンプレツシブノレ ナビエーストークス イクエイシヨンズ イン アド ベクシヨン ドミネーテイド フロウズ」 インターナショナノレ ジャーナノレ フォー ニュ 一メリカル メソーズ イン エンジニアリング 38卷 1995年 1475頁— 1506頁） 非特許文献 3 :松本純一,「気泡関数を用いた非圧縮性粘性流れ解析のための 2レべ ル -3レベル有限要素法」,応用力学論文集（土木学会）， 7卷, 2004年 8月， 339頁— 3 46頁

非特許文献 4：畑中勝守,「多段階有限要素法による非圧縮粘性流体の順 ·逆解析に 関する計算力学的研究」,中央大学博士論文, 1993年 3月

発明の開示 発明が解決しょうとする課題

[0014] ここで、上述した式（10)〜式（13)に示したように、 4段解法を使用して、この既知 解析物理量 uⁿカゝら未知解析物理量 uⁿ⁺¹を求めるためには、質量行列の逆行列が必 要になる。図 49は、従来の Rotating Cone問題の解析に用いる解析モデルを示す説 明図であり、図 50は、従来の Rotating Cone問題の解析に用いる初期条件の等高線 を示す説明図であり、図 51は、従来の整合質量行列によって計算された Rotating C one問題の解析結果の等高線を示す説明図である (非特許文献 4参照)。

[0015] 図 49の解析モデル 4900を初期状態 (0周時）とする。そして、質量行列の逆行列 を用いて、初期状態 (0周時)から、所定周、たとえば、 5周回転させると、図 50に示し た解析モデル 4900の初期状態における等高線モデル 5000は、図 51に示したよう な解析結果 (等高線モデル 5100)となる。

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

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

[0018] しかしながら、上述した集中質量行列を使用した場合には、集中質量行列が元の 質量行列と同一の行列ではなく近似行列となる。したがって、図 49に示した解析モ デル 4900の初期状態 (0周時)から、所定周、たとえば、 5周回転させると、図 50に 示した、解析モデル 4900の初期状態における等高線モデル 5000は、図 52に示し たような解析結果 (等高線モデル 5200)となる。このように、解析モデル 4900の 5周 回転後の等高線モデル 5200は、初期状態における等高線モデル 5000や図 51に 示した等高線モデル 5100に対して、大幅に変形しているため計算精度が悪ぐ解析 結果の信頼性が低 、と 、う問題があった。

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

課題を解決するための手段

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

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

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

cA_e 式 （1 4 )

2次元

G =

3次元 4 4 3 -

式 （1 5— 1 )

式 （1 5 - 2 ) ゾ Ω, 式 （1 5 3 )

[0023] 上記式（15— 1)〜式（15— 3)は、気泡関数要素の定式ィ匕をおこなうために必要と なる積分であり、く ·， ·〉 Ωは要素領域 Ω eの積分、 Aは要素領域 Ω eの面積 (体積)を

e e

示している。これにより、各要素の整合質量行列内の成分を足し合わせることによつ て近似された集中質量行列を用いずに、気泡関数要素の基底 (形状関数)が直交す る条件を導入して、各要素の高精度な対角質量行列を用いて解析することができる。 また、解析対象全域の対角質量行列やその逆行列も簡単に算出することができ、解 析対象の挙動を、未知の解析対象物理量を用いて高精度に解析することができる。 発明の効果

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

図面の簡単な説明

[0025] [図 1]図 1は、この発明の実施の形態に力かる数値解析装置のハードウェア構成を示 すブロック図である。

[図 2]図 2は、この発明の実施の形態にカゝかる数値解析装置の機能構成を示すブロッ ク図である。 [図 3]図 3は、この発明の実施の形態に力かる数値解析装置のデータ格納部が格納 して 、る解析物理量情報テーブルを示す説明図である。

[図 4]図 4は、この発明の実施の形態にカゝかる数値解析装置における数値解析処理 手順を示すフローチャートである。

[図 5]図 5は、 2次元の気泡関数要素を示す説明図である。

[図 6]図 6は、 3次元の気泡関数要素を示す説明図である。

[図 7]図 7は、直交基底を形成する三角形気泡関数の形状 φ を示す説明図である。

B

[図 8]図 8は、直交基底を形成する三角形気泡関数の形状 φ ²を示す説明図である。

B

[図 9]図 9は、直交基底を形成する三角形気泡関数の形状 φ を示す説明図である。

B

[図 10]図 10は、直交基底を形成する三角形気泡関数の形状 φ ²を示す説明図であ

B

る。

[図 11]図 11は、 2次元の気泡関数を示す説明図である。

[図 12]図 12は、 3次元の気泡関数を示す説明図である。

[図 13]図 13は、 Rotating Cone問題の解析における計算領域を示す説明図である。

[図 14]図 14は、 Rotating Cone問題の解析における解析モデルのメッシュを示す説明 図である。

[図 15]図 15は、 Rotating Cone問題の解析における解析モデルの流速（流況）を示す 説明図である。

[図 16]図 16は、 Rotating Cone問題の解析における初期条件を示す鳥瞰図である。

[図 17]図 17は、 Rotating Cone問題の解析における初期条件の等高線を示す説明図 である。

圆 18]図 18は、気泡関数 (式 (47) )での整合質量行列を使用した計算結果を示す 鳥瞰図である。

圆 19]図 19は、気泡関数 (式 (47) )での整合質量行列を使用した計算結果の等高 線を示す説明図である。

[図 20]図 20は、気泡関数式（(47) )での集中質量行列を使用した計算結果を示す 鳥瞰図である。

圆 21]図 21は、気泡関数 (式 (47) )での集中質量行列を使用した計算結果の等高 線を示す説明図である。

[図 22]図 22は、気泡関数 (式 (48) )での整合質量行列を使用した計算結果を示す 鳥瞰図である。

[図 23]図 23は、気泡関数 (式 (48) )での整合質量行列を使用した計算結果の等高 線を示す説明図である。

[図 24]図 24は、気泡関数 (式 (48) )での集中質量行列を使用した計算結果を示す 鳥瞰図である。

[図 25]図 25は、気泡関数 (式 (48) )での集中質量行列を使用した計算結果の等高 線を示す説明図である。

[図 26]図 26は、気泡関数 (式 (49) )での対角質量行列を使用した計算結果を示す 鳥瞰図である。

[図 27]図 27は、気泡関数 (式 (49) )での対角質量行列を使用した計算結果の等高 線を示す説明図である。

[図 28]図 28は、 Rotating Cone問題の解析における解析モデルのメッシュ（節点数： 4 38413、要素数： 2539200)を示す説明図である。

[図 29]図 29は、 Rotating Cone問題の解析における鉛直位置 =0の断面の初期条件 を示す鳥瞰図である。

[図 30]図 30は、 Rotating Cone問題の解析における鉛直位置 =0の部分の初期条件 の等高線を示す説明図である。

[図 31]図 31は、気泡関数 (式 (50) )での整合質量行列を使用した鉛直位置 =0の部 分の計算結果を示す鳥瞰図である。

[図 32]図 32は、気泡関数 (式 (50) )での整合質量行列を使用した鉛直位置 =0の部 分の計算結果の等高線を示す説明図である。

[図 33]図 33は、気泡関数 (式 (50) )での集中質量行列を使用した鉛直位置 =0の部 分の計算結果を示す鳥瞰図である。

[図 34]図 34は、気泡関数 (式 (50) )での集中質量行列を使用した鉛直位置 =0の部 分の計算結果の等高線を示す説明図である。

[図 35]図 35は、気泡関数 (式 (51) )での整合質量行列を使用した鉛直位置 =0の部 分の計算結果を示す鳥瞰図である。

[図 36]図 36は、気泡関数 (式 (51) )での整合質量行列を使用した鉛直位置 =0の部 分の計算結果の等高線を示す説明図である。

[図 37]図 37は、気泡関数 (式 (51) )での集中質量行列を使用した鉛直位置 =0の部 分の計算結果を示す鳥瞰図である。

[図 38]図 38は、気泡関数 (式 (51) )での集中質量行列を使用した鉛直位置 =0の部 分の計算結果の等高線を示す説明図である。

[図 39]図 39は、気泡関数 (式 (52) )での対角質量行列を使用した鉛直位置 =0の部 分の計算結果を示す鳥瞰図である。

[図 40]図 40は、気泡関数 (式 (52) )での対角質量行列を使用した鉛直位置 =0の部 分の計算結果の等高線を示す説明図である。

[図 41]図 41は、アイソパラメトリック座標系 rを示す説明図（その 1)である。

[図 42]図 42は、アイソパラメトリック座標系 rを示す説明図（その 2)である。

[図 43]図 43は、式 (81)、（93)を用いた三角形気泡関数の形状 φ を示す説明図で

B

ある。

圆 44]図 44は、式 (81)、（93)を用いた三角形気泡関数の形状 φ ²を示す説明図で

B

ある。

[図 45]図 45は、式（132)を用いた三角形 3レベル気泡関数の形状を示す説明図で ある。

[図 46]図 46は、式 (81)、（93)を用いた三角形気泡関数と式（132)を用いた三角形

3レベル気泡関数の積の形状を示す説明図である。

[図 47]図 47は、従来の 2次元の気泡関数要素を示す説明図である。

[図 48]図 48は、従来の 3次元の気泡関数要素を示す説明図である。

[図 49]図 49は、従来の Rotating Cone問題の解析に用いる解析モデルを示す説明図 である。

[図 50]図 50は、従来の Rotating Cone問題の解析に用いる初期条件の等高線を示す 説明図である。

[図 51]図 51は、従来の整合質量行列によって計算された Rotating Cone問題の解析 結果の等高線を示す説明図である。

[図 52]図 52は、従来の集中質量行列によって計算された Rotating Cone問題の解析 結果の等高線を示す説明図である。

符号の説明

[0026] 200 数値解析装置

201 データ格納部

202 第 1の取得部

203 第 2の取得部

204 生成部

205 第 1の算出部

206 第 2の算出部

207 解析部

発明を実施するための最良の形態

[0027] 以下に添付図面を参照して、この発明にかかる直交基底気泡関数要素数値解析 方法、直交基底気泡関数要素数値解析プログラムおよび直交基底気泡関数要素数 値解析装置 (以下、単に、「数値解析装置」という）の好適な実施の形態を詳細に説 明する。

[0028] (数値解析装置のハードウェア構成）

この実施の形態に力かる数値解析装置のハードウェア構成について説明する。図 1 は、この発明の実施の形態に力かる数値解析装置のハードウェア構成を示すブロッ ク図である。数値解析装置は、 CPU101と、 ROM102と、 RAM103と、 HDD (ハー ドディスクドライブ） 104と、 HD (ノヽードディスク） 105と、 FDD (フレキシブルディスクド ライブ） 106と、着脱可能な記録媒体の一例として FD (フレキシブルディスク） 107と、 ディスプレイ 108と、 IZF (インターフェース） 109と、キーボード 110と、マウス 111と 、スキャナ 112と、プリンタ 113と、を備えている。また、各構成部は、バス 100によつ てそれぞれ接続されて ヽる。

[0029] ここで、 CPU101は、数値解析装置の全体の制御を司る。 ROM102は、ブートプ ログラムなどのプログラムを記 '慮している。 RAM103は、 CPU101のワークエリアとし て使用される。 HDD104は、 CPU101の制御にしたがって HD105に対するデータ のリード Zライトを制御する。 HD105は、 HDD104の制御で書き込まれたデータを feす。。

[0030] FDD106は、 CPU101の制御にしたがって FD107に対するデータのリード Zライ トを制御する。 FD107は、 FDD106の制御で書き込まれたデータを記憶したり、 FD 107に記憶されたデータを数値解析装置に読み取らせたりする。着脱可能な記録媒 体として、 FD107のほ力、、 CD— ROM (CD— R、 CD— RW)、 MO、 DVD (Digital Versatile Disk)、メモリーカードなどであってもよい。ディスプレイ 108は、カーソル、 アイコンあるいはツールボックスをはじめ、文書、画像、機能情報などのデータを表示 する。このディスプレイ 108は、たとえば、 CRT, TFT液晶ディスプレイ、プラズマディ スプレイ等を採用することができる。

[0031] IZF109は、通信回線を通じてインターネットなどのネットワークに接続され、このネ ットワークを介して他の装置に接続される。そして、 IZF109は、ネットワークと内部の インターフェースを司り、外部装置からのデータの入出力を制御する。 IZF109には 、たとえばモデムや LANアダプタなどを採用することができる。

[0032] キーボード 110は、文字、数字、各種指示などの入力のためのキーを備え、データ の入力をおこなう。また、タツチパネル式の入力パッドやテンキーなどであってもよい 。マウス 111は、カーソルの移動や範囲選択、あるいはウィンドウの移動やサイズの変 更などをおこなう。ポインティングデバイスとして同様に機能を備えるものであれば、ト ラックボールやジョイスティックなどであってもよい。

[0033] スキャナ 112は、画像を光学的に読み取り、数値解析装置内に画像データを取り 込む。また、プリンタ 113は、画像データや文書データを印刷する。プリンタ 113には 、たとえば、レーザプリンタやインクジェットプリンタを採用することができる。

[0034] (数値解析装置の機能的構成）

この発明の実施の形態に力かる数値解析装置の機能的構成について説明する。 図 2はこの発明の実施の形態に力かる数値解析装置の機能構成を示すブロック図で ある。図 2に示すように、数値解析装置 200は、データ格納部 201と、第 1の取得部 2 02と、第 2の取得部 203と、生成部 204と、第 1の算出部 205と、第 2の算出部 206と 、解析部 207と、力 構成されている。

[0035] データ格納部 201は、解析対象の数値解析に用いる解析物理量情報テーブルを 有する。ここで、解析物理量情報テーブルについて説明する。図 3は、この発明の実 施の形態に力かる解析物理量情報テーブルを示す説明図である。

[0036] 図 3において、解析物理量情報テーブルには、解析対象についての、解析物理量 I Dと、解析物理量名（汚染物質、熱、流体、構造など)と、既知解析物理量とが格納さ れている。また、既知解析物理量としては、既知物性値と、境界解析物理量と、初期 解析物理量とが格納されて 、る。

[0037] 解析物理量 IDおよび解析物理量名は、ユーザの入力操作によって指定される情 報であり、解析物理量 IDや解析物理量名が指定されると、関連付けられている既知 物性値、境界解析物理量および初期解析物理量を抽出することができる。

[0038] 既知物性値は、ユーザの入力操作によって解析物理量 IDまたは解析物理量名が 指定された解析対象が、どのような物質であるかを決定する情報である。既知物性値 としては、たとえば、拡散係数、密度、比熱、熱伝導係数、粘性係数、水底面摩擦係 数、ヤング率、ポアソン比などが挙げられる。

[0039] 境界解析物理量とは、解析対象の各要素からなる解析領域 (三角形や四面体に分 割したメッシュ）がどのような条件下で解析がおこなわれるかと!/、う解析条件をあらわ す情報である。境界解析物理量としては、たとえば、物質濃度、温度、流速、圧力、 流量、水深、変位などが挙げられる。

[0040] 境界解析物理量については、解析対象がたとえば汚染物質や熱の場合、汚染物 質や熱の発生 (供給)源の部分で、物質濃度 =発生 (供給)源値、温度 =発生 (供給 )源値となる。また、解析対象が流体で境界が壁ならば、流速 =0、流量 =0、解析対 象が構造で境界が物体を固定ならば、変位 =0に設定することができる。定常問題（ 解析物理量が時間に依存せず時々刻々と変化しない問題)であれば、上述した既知 物性値および境界解析物理量によって、解析することができる。

[0041] 非定常問題 (解析物理量が時間に依存し時々刻々と変化する問題)の場合は、さら に初期条件となる初期解析物理量が必要になる。初期解析物理量は、解析開始時 の解析物理量の値である。また、データ格納部 201において、図示はしないが、初期 設定情報として、時間方向に対して、逐次解析計算をおこなうために時間増分量、計 算ステップ数などの情報も、解析物理量情報テーブルに格納されて 、る。

[0042] さらに、解析対象範囲に分割形成されるメッシュの分割幅、節点数、要素数、節点 番号に対応した節点座標値および要素番号に対応した要素結合情報などのメッシュ 情報も格納されている。このメッシュは、解析対象範囲が 2次元空間である場合は、 三角形の形状、解析対象範囲が 3次元空間である場合は、四面体の有限な形状の 要素となる。メッシュは、図 1に示したスキャナ 112により読み取ることによつても取り込 むことができる。このデータ格納部 201は、具体的には、例えば、図 1に示した RAM 103、HD105または FD107によりその機能を実現する。

[0043] 第 1の取得部 202は、解析対象の既知解析物理量を取得する。具体的には、入力 装置 (たとえば図 1に示したキーボード 110またはマウス 111)を操作することにより、 解析対象となる物理量の解析物理量 IDを選択し、データ格納部 201から既知解析 物理量を抽出する。この第 1の取得部 202は、具体的には、たとえば図 1に示した R OM102、 RAM103、 HD105、 FD107などに格納されたプログラムを CPU101が 実行することによって、また図 1に示した IZF109によってその機能を実現する。

[0044] 第 2の取得部 203は、解析対象の各要素 (メッシュ）の整合質量行列を取得する。

第 2の取得部 203は、具体的には、たとえば、図 1に示した ROM102、 RAM103、 HD105、FD107などに格納されたプログラムを CPU101が実行することによって、 また図 1に示した IZF109によってその機能を実現する。

[0045] 生成部 204は、解析対象の各要素の気泡関数に基づいて、第 2の取得部 203によ つて取得された各要素の整合質量行列を対角化した、各要素の対角質量行列を生 成する。生成部 204は、具体的には、たとえば、図 1に示した ROM102、 RAM103 、 HD105、 FD107などに格納されたプログラムを CPU101が実行することによって 、また図 1に示した IZF109によってその機能を実現する。

[0046] 第 1の算出部 205は、生成部 204によって生成された各要素の対角質量行列に基 づいて、解析対象全域の対角質量行列を算出する。第 1の算出部 205は、具体的に は、たとえば、図 1に示した ROM102、 RAM103、 HD105、 FD107などに格納さ れたプログラムを CPU101が実行することによって、また図 1に示した IZF109によ つてその機能を実現する。

[0047] 第 2の算出部 206は、第 1の算出部 205によって算出された、解析対象全域の対角 質量行列の逆行列を算出する。第 2の算出部 206は、具体的には、たとえば、図 1に 示した ROM102、 RAM103、 HD105、 FD107などに格納されたプログラムを CP U101が実行することによって、また図 1に示した IZF109によってその機能を実現 する。

[0048] 解析部 207は、第 1の取得部 202によって取得された既知解析物理量と、生成部 2 04によって生成された各要素の対角質量行列と、に基づいて、解析対象の挙動を 解析する。具体的には、第 1の取得部 202によって取得された既知解析物理量と、 第 1の算出部 205によって算出された解析対象全域の対角質量行列と、第 2の算出 部 206によって算出された逆行列と、に基づいて、解析対象の挙動を解析する。な お、解析対象全域の対角質量行列を必要とせず、各要素の対角質量行列のみを用 いて解析対象の挙動を解析する解法があるが（たとえば、 O.C.Zienkiewicz, R丄. Tay lor, b.J.Sherwin and J.Peiro, On Discontinuous Galerkin Methods , Internati onal Journal for Numerical Methods in Engineering, Vol.58, 2003, pp.1119— 1 148参照）、この場合には、生成部 204と第 1の算出部 205を同一視して、第 2の算出 部 206へ進むことも可能である。解析部 207は、具体的には、たとえば、図 1に示した ROM102、 RAM103、 HD105、 FD107などに格納されたプログラムを CPU101 が実行することによって、また図 1に示した IZF109によってその機能を実現する。

[0049] 本実施の形態にカゝかる数値解析装置 200の数値解析処理手順にっ ヽて説明する 。図 4はこの発明の実施の形態にカゝかる数値解析装置 200の数値解析処理手順を 示すフローチャートである。まず、第 1の取得部 202により、解析対象の既知解析物 理量を取得する (ステップ S401)。つぎに、第 2の取得部 203により、各要素（要素レ ベル）の整合質量行列を取得する（ステップ S402)。

[0050] 続いて、要素ごとに、気泡関数を積分し (ステップ S403)、各要素（要素レベル)の 整合質量行列に、ステップ S403で積分された値を代入することによって、各要素（要 素レベル)の対角質量行列を算出する (ステップ S404)。つぎに、各要素（要素レべ ル)の対角質量行列の総和 (重ね合せ）により、解析対象全域の対角質量行列を算 出する (ステップ S405)。

[0051] そして、この対角質量行列の逆行列を算出し (ステップ S406)、解析対象の既知解 析物理量と、解析対象全域の対角質量行列と、その逆行列とに基づいて、解析対象 の挙動を解析する。具体的には、たとえば、上記式（10)〜式（13)に代入することに より、解析することができる。

[0052] 数値解析装置 200の具体的な解析処理内容にっ 、て説明する。

<1.直交基底気泡関数要素 >

質量行列 Mは、有限要素法による解析 (有限要素解析）について、ある解析物理 量を補間した関数に重み関数を掛け、対象領域で積分を行った場合に現れる。気泡 関数要素を使用した場合には、解析対象となる任意領域を区分的に三角形 (四面体 )形状にて要素数 Nで分割した要素群を全体系に重ね合せる作業によって形成され e

、要素領域での積分を用いて下記式（16)のようにあらわすことができる。

[0053] [数 7]

M = {*,$^T)n_e=∑^^e) , i = l---N + 2, = l---N + 2

e=l e=l · · ·式 （ 1 6 )

[0054] 式（16)に対して、気泡関数要素の基底 (形状関数)が直交するためには、下記式（ 18)、式（20)を満たす必要がある。なお、式（18)の導入に際して、式（17)を用いる 。また、式（20)の導入に際して、式（19)を用いる。

[0055] [数 8]

〈 。'^〉 Ω,

, a = l-''N + l · [0056] [数 9] 式 （1 7) より

(1,ΦΒ =

'式 ( 1 8)

[0057] [数 10] :(1,^Β>Ω. = 0 , α≠β , β = 1·--Ν+1

(N+l):

式 （1 9)

[0058] [数 11] 式 （1 9) より

式 （20)

[0059]

[0060] 式 （2 1)

3次元

G-

[0061] < 2.直交条件を満たす気泡関数 >

式（ 18)、式（20)を満たすことができる気泡関数として下記式（22)で表される気泡 関数を提案する。

[0062] [数 13]

ΟίΐψΒ + Οί2ψΒ

ひ 1 + + 1 式 （2 2) α ₂は未知量、

はそれぞれ異なる気泡関数である _c 式（18)に式（22)を代入することにより下記式（23)を得る。

[0064] [数 14] （ 2 3 )

, βί = (ΦΒ, l)n_e + {ΦΒ, l)n_e - 2(φ_Β, ΦΒ)Ω,，

あるいは

βί = ((ΦΒ, l)n_e + (ΦΒ, 1)Ω«一 2(φ_Β, β)η Α ¹ , 5 = ((ΦΒ, + (ΦΒ, 1〉¾― 2{>_B, _c) Α ¹ ,

[0065] 式（20)に式（22)を代入することにより下記式（24)を得る。

[0066] [数 15] ひ 2 = 71ひ 1 + 72 式 （24)

,

7ι

あるいは

一 〈 , Κ¹ - c

H —¹ - c'⁷²一 {ΦΒ, Ω ¹ - C

[0067] 式（23)に式（24)を代入することにより、最終的に _αを未知量とした下記式（25)を

1

得る。

[0068] [数 16] aa + ba_{1 +}c = 0 · . .式 （_{2 5)} α =

+ /¾ + 72

[0069] 式（25)は _α に関する二次方程式なので、解の公式より下記式（26)にて未知量 a を求めることができる。

1

[0070] [数 17]

-δ士 Vも²― 4ac

αι =

2α 式 （26)

[0071] < 3.直交基底を形成する具体的な気泡関数 >

直交基底を形成する具体的な気泡関数について説明する。図 5は、 2次元の気泡 関数要素を示す説明図であり、図 6は、 3次元の気泡関数要素を示す説明図である。

[0072] 直交基底を形成する具体的な気泡関数を導出するため、図 5および図 6に示した 三角形 (四面体)の要素領域を、重心点を用いて 3(4)つの小三角形 (小四面体) w,i = 1···Ν+1に分割する ξ乗気泡関数 (0< ξ <∞)を用いる。 ξ乗気泡関数は小 三角形 (小四面体)ごとにアイソパラメトリック座標系 {r,_S}({r,s,t})を用いて次の式（2 7)および式（28)のように定義される。

[0073] [数 18]

2次元

3^ (1 - r - in W

ΦΒ = ^ in W2

3« in 式 （27)

3次元

『4« (1 t)ミ in wx

4 in W2

4 in wz

4 in W 式 （28)

[0074] これらの気泡関数の要素領域における積分値は、以下の式（29)によって求める: とがでさる。

[0075] [数 19]

式 （29) [0076] 式（22)を、 ξ乗気泡関数を用いて下記式（30)のように定義する。

[0077] [数 20] is l . . .式 （ 30 ) ，ミ，ミは、 各べき乗の値

1

[0078] 式（30)の各べき乗の 65値を下記式（31)のように設定すると二次元（三角形）、三次

一 4 3

元（四面体）の気泡関数に対して α , a は下記式（32)〜式（35)のようになる。

1 2

[0079] [数 21] 式

[0080] [数 22]

2次元 （三角形）

= 3.1871016·· · , ひ 2 = -4.5742685··

― 2α 式 （32) -b― ²― 4ac

ひ 1 = = 3.1457828·· · , ひ ₂ = -3.9426812·

2a 式 （33)

3次元 （四面体）

一 b + ― 4ac

αι = = 3„.1一27一797—4 ,- a₂ = -3.8884541 ·

式 (34)

—b一 δ²― 4ac

ひ i = = 2.6785040· ひ ₂ = -4.0779944·,·

式 （35)

[0081] 式（30)の各べき乗の値を下記式（36)のように設定する（式（31)とは異なった値) と二次元（三角形）、三次元（四面体)の気泡関数に対してひ ，

1 2は下記式（37)〜 式（40)のようになる。

[0082] [数 23] 式 （36

[0083] [数 24] 2次元 （三角形）

-6 + fe² -4ac

。₂ = -0.6433844-

2o 式 （3 7)

—b― o*― 4ac

Qi = = 1.2696544·· · ₂ = -2.2134591··

2a 式 （38)

3次元 （四面体）

— & + _ 4ac

~ ½ = 0.7333511 · ひ ₂ = -1.6387018··

式 （3 9)

—b― ― 4ac

QI = = 1.2578218 - a₂ = -2.2006346··

式 (40)

[0084] 図 7は、式 (32)の a , aを使用した直交基底を形成する三角形気泡関数の形状

1 2

Φ を示す説明図であり、図 8は、式 (32)の a ,αを使用した直交基底を形成する三

Β 1 2

角形気泡関数の形状 Φ ²を示す説明図である。また、図 9は、式 (37)の a ,αを使

Β 1 2 用した直交基底を形成する三角形気泡関数の形状 φ

Βを示す説明図であり、図 10は

、式 (37)の a ,αを使用した直交基底を形成する三角形気泡関数の形状 φ ²を示

1 2 Β す説明図である。

[0085] 式（32)、式（37)および図 7〜図 10に示すように、式（21)を満たす気泡関数 (形状 )は無数に存在するため、直交基底を形成する具体的な気泡関数やその形状には 殆ど意味がなぐ直交基底気泡関数要素の本質は気泡関数要素の基底が直交する 条件式 (21)を発明した部分にある。

[0086] <4.直交基底気泡関数要素の要素レベルの質量行列 >

下記式 (41)〜式 (46)に、 2次元（三角形）、 3次元（四面体）における通常の気泡 関数要素の要素レベルの整合質量行列、集中質量行列、直交基底気泡関数要素 の要素レベルの質量行列 (対角質量行列)を示す。

[0087] [数 25] 二角形気泡関数要素

要素レベルの整合質量行列

一 2 — 2 -2 3 1 1 1 -3

-2 一 2 — 2 3 1 1 1 -3 一 2 — 2 — 2 3 + ΙΚ

1 1 1 一 3

3 3 3 0 一 3 -3 一 3 9

o

式 （4 1 ) 要素レベルの集中質量行列 （集中化した場合)

(Ν+2 \ o

Σ ^Mi†

i=i ノ o o o o

式 （4 2 ) 要素レベルの対角質量行列 （直交基底となる気泡関数)

0 ん.

12 0 0

0 0 ん

12 0

0 0 0

式 （4 3 ) 式（43)の要素レベルの対角質量行列 M^(e)は、生成部 204により、式（41)の要素 レベルの整合質量行列 _M ^(e)のく: L, φ > Qおよび II φ II ² Ωに、気泡関数 φ の積

ij B e B e B 分値 CAを代入することによって生成される。そして、第 1の算出部 205により、各要 e

素の要素レベルの対角質量行列 M^(e)の総和（重ね合せ)をとることによって、解析対 象全域の対角質量行列 Mを算出することができる。また、第 2の算出部 206により、 対角質量行列 Mの逆行列 M—¹を算出することができる。そして、解析部 207により、既 知解析物理量 u、対角質量行列 Mおよび逆行列 M—¹を、上述した式（10)〜式（13) に代入することによって、解析対象の挙動を解析することができる

[数 26] 四面体気泡関数要素

要素レベルの整合質量行列

-2 一 2 一 2 一 2 4 1 1 1 1 -4 一 2 一 2 一 2 一 2 4 ^ o 1 1 1 1 -4

+ 〈1 〉 Ω« 一 2 一 2 一 2 一 2 4 + ^¾ん_ S∞∞ , o 1 1 1 1 一 4

一 2 -2 -2 一 2 4 1 1 1 1 一 4

4 4 4 4 0 一 4 -4 -4 -4 16

o o o o o 式 （4 4 ) 要素レベルの集中質量行列 （集中化した場合）

〈*, 〉 Ω,

« diag ^Mif

Ac

4 0 0 0 0 -1 0 0 0 0

0 0 0 0 0 -1 0 0 0

0 0 0 0 0 0 -1 0 0

0 0 0 0 0 0 0 -1 0

0 0 0 0 0 0 0 0 0 4

式 （4 5 ) 要素レベルの対角質量行列 （直交基底となる気泡関数)

式 ( 4 6 ) 式 (46)の要素レベルの対角質量行列は、生成部 204により、式 (44)の要素レべ ルの整合質量行列の〈1, φ > Qおよび II φ II ² Ω 〖こ、気泡関数 φ の積分値 CA

B e B e B e を代入することによって生成される。そして、第 1の算出部 205により、各要素の要素 レベルの対角質量行列 M^(e)の総和（重ね合せ)をとることによって、解析対象全域の 対角質量行列 Mを算出することができる。また、第 2の算出部 206により、対角質量 行列 Mの逆行列 M—¹を算出することができる。そして、解析部 207により、既知解析物 理量 u、対角質量行列 Mおよび逆行列 M—¹を、上述した式（10)〜式（13)に代入す ること〖こよって、解析対象の挙動を解析することができる。

[0091] < 5.気泡関数について >

気泡関数について説明する。図 11は、 2次元の気泡関数を示す説明図であり、図 12は、 3次元の気泡関数を示す説明図である。気泡関数 φ は、形状関数の定義（

B

重心点が 1、その他の点では 0となる関数)を満たせば任意の関数を選択することが できる。検証の比較計算として、最も多く使用されている以下に示す二つの気泡関数 および直交基底を形成する気泡関数を用いる。図 11および図 12に示したァイソパラ メトリック座標系 {r,s} ( {r,s,t})を用いて、式 (47)〜式（52)のように定義される。

[0092] [数 27]

二次元 （三角形）

多く使用されている気泡関数 1

φ_Β = 27 {l - r - s) r s , (φ_Β, l)si_e =

• · ·式 （4 7 )

• · ·式 （4 8 ) = -A,

• · ·式 （4 9

[0093] [数 28] 次元 （四面体）

多く使

φ_Β =

式 （5 0 ) 多く使用されている気泡関数 2

4 (1— r— s— i) in w\

4 r in u>2

， {φΒ, 1)η. = ^Α_Β , = -A_e

4 s in ΐϋ3

4 t in Wi

式 （5 1 ) 直交基底となる気泡関数

ΟίιψΒ + X2<pB + Φβ

ひ 1 + <¾ + 1

式 （5 2 )

[0094] < 6. Rotating Cone問題 >

発明手法の検証計算として、数値解析手法の性能を比較検討するための有名な ベンチマーク問題である Rotating Cone問題 (物質濃度の移流問題）の解析をおこな う（非特許文献 4参照)。図 13は、 Rotating Cone問題の解析における計算領域を示 す説明図であり、図 14は、 Rotating Cone問題の解析における解析モデルのメッシュ (節点数: 4921、要素数: 9600)を示す説明図であり、図 15は、 Rotating Cone問題の 解析における解析モデルの流速 (流況)を示す説明図である。

[0095] 図 16は、 Rotating Cone問題の解析における初期条件 (初期濃度分布）を示す鳥瞰 図であり、図 17は、 Rotating Cone問題の解析における初期条件 (初期濃度分布）の 等高線を示す説明図である。この Rotating Cone問題は、非定常移流方程式におい て図 13の計算領域に、図 14にて図 15の定常な一定流速 (流況)を仮定し、図 16お よび図 17のようなある物理量 (汚染物質濃度など)の移流状況を解析する問題である

[0096] 実際の物理量 (汚染物質濃度など）は、濃度が拡散しながら流況にそって運ばれる 力 拡散の効果をなくして計算を行った場合 (移流方程式を計算した場合）には、濃 度が流れに乗って計算領域を周回したとき、物理量の分布が、理想的には、物理量 を移流させる作用しかな 、ので元の初期分布と一致しなければなら 、と 、う問題であ る。この問題により数値解析を行えば、対象とした手法が持つ近似誤差を評価するこ とができ、計算精度の検証ができる。

[0097] 図 18は、気泡関数 (式 (47) )での整合質量行列を使用した計算結果を示す鳥瞰 図であり、図 19は、気泡関数 (式 (47) )での整合質量行列を使用した計算結果の等 高線を示す説明図である。また、図 20は、気泡関数 (式 (47) )での集中質量行列を 使用した計算結果を示す鳥瞰図であり、図 21は、気泡関数 (式 (47) )での集中質量 行列を使用した計算結果の等高線を示す説明図である。

[0098] 図 22は、気泡関数 (式 (48) )での整合質量行列を使用した計算結果を示す鳥瞰 図であり、図 23は、気泡関数 (式 48)での整合質量行列を使用した計算結果の等高 線を示す説明図である。また、図 24は、気泡関数 (式 (48) )での集中質量行列を使 用した計算結果を示す鳥瞰図であり、図 25は、気泡関数 (式 (48) )での集中質量行 列を使用した計算結果の等高線を示す説明図である。また、図 26は、気泡関数 (式（ 49) )での対角質量行列を使用した計算結果を示す鳥瞰図であり、図 27は、気泡関 数 (式 (49) )での対角質量行列を使用した計算結果の等高線を示す説明図である。

[0099] 解析方法としては、検証対象としたすベての気泡関数に対して、空間方向の離散 化には数値解析の安定ィ匕を考慮した気泡関数要素安定ィ匕法 (非特許文献 3参照)を 使用し、時間方向の離散化には 4段解法を、微小時間 A tは π Ζ400を採用した。図 18、図 19、図 22および図 23は、気泡関数 (式 (47) )、（式 (48) )での整合質量行列 を使用した計算結果である。

[0100] この計算結果に対して、気泡関数 (式 (47)、式 (48) )での集中質量行列を使用し た計算結果、図 20、図 21、図 24および図 25は、物理量（円錐)の進行方向の後方 に振動が発生しており、円錐の分布状況も崩れている結果となっている。これは、濃 度の物理的に意味のない減衰、広がりが著しぐ濃度が最大になるべき地点も大きく ずれており、数値解析上、信頼性の低い計算結果となっている。

[0101] 一方、気泡関数 (式 (49) )での対角質量行列を使用した計算結果は、円錐の進行 方向の後方に振動が発生しておらず、その分布も崩れていなぐ気泡関数 (式 (47) ) 、式 48) )での整合質量行列を使用した計算結果と同等の計算精度を保っており、数 値解析上、信頼性の高い計算結果が得られている。 [0102] < 7. Rotating Cone問題による 3次元解析の検証 >

発明手法の 3次元解析の検証として、 2次元と同様に Rotating Cone問題 (物質濃度 の移流問題)の解析をおこなう（非特許文献 4参照)。図 28は、 Rotating Cone問題の 解析における解析モデルのメッシュ（節点数: 438413、要素数： 2539200)を示す説明 図である。図 28は、図 13の 2次元解析モデルに鉛直長さ = 2 (— 1〜1)を持たせた モデルである。解析モデルの流速 (流況)および濃度分布は、鉛直方向には変化が なく 2次元解析モデルの場合と同様である。

[0103] 図 29は、 Rotating Cone問題の解析における鉛直位置 =0の断面の初期条件（初 期濃度分布）を示す鳥瞰図であり、図 30は、 Rotating Cone問題の解析における鉛直 位置 =0の部分の初期条件 (初期濃度分布)の等高線を示す説明図である。

[0104] 図 31は、気泡関数 (式（50) )での整合質量行列を使用した鉛直位置 =0の部分の 計算結果を示す鳥瞰図であり、図 32は、気泡関数 (式 (50) )での整合質量行列を使 用した鉛直位置 =0の部分の計算結果の等高線を示す説明図である。また、図 33は 、気泡関数 (式 (50) )での集中質量行列を使用した鉛直位置 = 0の部分の計算結果 を示す鳥瞰図であり、図 34は、気泡関数 (式 (50) )での集中質量行列を使用した鉛 直位置 =0の部分の計算結果の等高線を示す説明図である。

[0105] 図 35は、気泡関数 (式（51) )での整合質量行列を使用した鉛直位置 =0の部分の 計算結果を示す鳥瞰図であり、図 36は、気泡関数 (式 (51) )での整合質量行列を使 用した鉛直位置 =0の部分の計算結果の等高線を示す説明図である。また、図 37は 、気泡関数 (式 (51) )での集中質量行列を使用した鉛直位置 = 0の部分の計算結果 を示す鳥瞰図であり、図 38は、気泡関数 (式 (51) )での集中質量行列を使用した鉛 直位置 =0の部分の計算結果の等高線を示す説明図である。また、図 39は、気泡関 数 (式 (52) )での対角質量行列を使用した鉛直位置 = 0の部分の計算結果を示す 鳥瞰図であり、図 40は、気泡関数 (式 (52) )での対角質量行列を使用した鉛直位置 =0の部分の計算結果の等高線を示す説明図である。

[0106] 解析方法としては、検証対象としたすベての気泡関数に対して、空間方向の離散 化には数値解析の安定ィ匕を考慮した気泡関数要素安定ィ匕法 (非特許文献 3参照)を 使用し、時間方向の離散化には 4段解法を、微小時間 A tは π Ζ600を採用した。図 31、図 32、図 35および図 36は、気泡関数 (式（50)、式（51) )での整合質量行列を 使用した計算結果である。

[0107] この計算結果に対して、気泡関数 (式 (50)、式 (51) )での集中質量行列を使用し た計算結果、図 33、図 34、図 37および図 38は、物理量（円錐)の進行方向の後方 に振動が発生しており、円錐の分布状況も崩れている結果となっている。これは、濃 度の物理的に意味のない減衰、広がりが著しぐ濃度が最大になるべき地点も大きく ずれており、数値解析上、信頼性の低い計算結果となっている。

[0108] 一方、気泡関数 (式 (52) )での対角質量行列を使用した計算結果は、円錐の進行 方向の後方に振動が発生しておらず、その分布も崩れていなぐ気泡関数 (式 (50) 、式 (51) )での整合質量行列を使用した計算結果と同等の計算精度を保っており、 数値解析上、信頼性の高い計算結果が得られている。

[0109] < 8. 1次元の直交基底気泡関数要素と要素レベルの質量行列 >

以上、 2次元、 3次元の直交基底気泡関数要素を発明し、その効果を示した。有限 要素法の数値計算では、実用上、問題となる解析は、殆どの場合、 2次元と 3次元で あるが、直交基底気泡関数要素は、 1次元でも定義することができる。ここでは、 1次 元の気泡関数要素、 ξ乗気泡関数を示し、直交条件式 (14) (または式 (21) )を 1次 元〜 3次元まで空間次元数 Νを用いて統一的に記述する。また、一次元 (線気泡関 数要素）の要素レベルの質量行列を示す。

[0110] 図 41は、アイソパラメトリック座標系 rを示す説明図（その 1)である。図 41に示したァ イソパラメトリック座標系 rを、式（1)、（2)の表現で記述すると、 1次元の Φ^τ, u^T, ¥ a は下記式（53)〜（55)になる。

[0111] [数 29]

• · ·式 （5 3 )

• · ·式 （5 4 )

• · ·式 （5 5 ) 図 42は、アイソパラメトリック座標系 rを示す説明図（その 2)である。図 42に示したァ イソパラメトリック座標系 rにおいて、中間点を用いて、線要素領域を w、 wに分割し、 気泡関数の要素領域での積分値が式 (29)となる 1次元の ξ乗気泡関数は下記式 ( 56)のように定義される。

一 9

[0113] [数 30]

1次元

2ミ (1― 、 in ti»i

一

m υ>2 式 （56)

[0114] 式（56)、 (27)、 (28)の （r?)乗気泡関数 (0< r? <∞)は、下記式（57)に示すよ うに、気泡関数を空間で偏微分した積分値も統一的に記述でき、下記式 (58)、 (59) のような積分値の関係を持っている。

[0115] [数 31] dik ' ^e 1、 dxk dxiメ

式 （5 7)

■Ν+1

N+1 式 （58)

(1,∑ α^)_{Ωί )} a = l-.-N + l , 0< _fci <

N +l k_d=l

式 （59)

[0116] 式（57)の x,x=x · χはそれぞれ空間方向を表し、式（59)の a は実数を、 K

k 1 1 N kd D は正の整数を示している。式（59)は式（58)より誘導している。 C=(N+1)/(N + 2)とし、空間次元数 Nを使用して直交基底気泡関数要素の条件式（14) (または式（ 21))を 1次元〜 3次元まで統一的に記述すると下記式（60)になる。

[0117] [数 32]

N+1

Ν + 2 式 (60)

[0118] 下記式 (61)〜（63)に、 1次元 (線）における通常の気泡関数要素の要素レベルの 整合質量行列、集中質量行列、直交基底気泡関数要素の要素レベルの質量行列（ 対角質量行列)を示す,

[0119] [数 33] 線気泡関数要素

要素 +レベルの整合質量行列

A,,

3 ， 6

6 3

0 0 一 2 一 2 1 1 一 2

一 2 一 2 1 1 一 2

2 2

一 2 一 2 4

式 （6 1 ) 要素レベルの集中質量行列 （集中化した場合)

式 （6 2) 要素レベルの対角質量行列 （直交基底となる気泡関数:

c

6 0

― ―

6 0

0 A,

式 （6

[0120] < 9.拡張された直交基底気泡関数要素 >

気泡関数要素を使用した流体解析、構造解析、固有値解析などでは、一般に一次 要素の形状関数 Ψ αの積分値に加えて、下記式 (64)〜（66)の気泡関数の積分値 が必要になる。

[0121] [数 34] a,< , a= l---N + l . . .式 （6 4)

IliM . . .式 （ 6 5 ) d^xk dxi ^¾ . . .式 （ 6 6 ) [0122] 下記式 (67)と式 (64)より、下記の関係式 (68)を得る。

[0123] [数 35]

• · ·式（6 7 )

• · '式（6 8 )

[0124] 気泡関数として、下記式 (69)の関係を持っている気泡関数の場合には、必要とな る気泡関数の積分値 (64)〜（66)の式 (64)は、式（70)に代用することができる。

[0125] [数 36]

• · '式 （6 9 )

• · ·式（7 0 )

[0126] 式 (69)の関係式は、非特許文献 3で現れるような特殊な気泡関数を導入しない限 り、例えば式 (47)、（48)、（50)、（51)などの一般的によく使用されている気泡関数 で成り立つので、直交条件式（21)の導出仮定で現れる式（17)、（19)についても使 用しており、式（22)は式 (69)を満たす気泡関数であると仮定している。式（22)右辺 の気泡関数として式（30)を採用した場合には、式（59)より式（22)は式 (69)を満足 していることがわ力る。

[0127] 直交基底気泡関数要素を使用した場合には、式（70)、（65)、（66)において、式（ 70)、（65)は直交条件式 (60)よりわかる力 式 (66)は、直交条件により決定するこ とができない。式 (56)、（27)、（28)の ξ乗気泡関数を使用して、直交基底気泡関 数を形成した場合、実際に、式 (66)は、下記式（71)で表される。

[0128] [数 37] ^' 3φ_α、

H , ξ φひ

=ι ^dxk dxi 式 （7 1) 式 (7 2)

3 -4 ,

[0129] 式（71)をみると、 D (

1 ,62,6 )の値が、変数 、 の関数となっている。実

3 1 2 3

際に、式 (31 (36)の値を用いて、式 (72)を計算 ( ξ , ξ ) = ( ξ , ξ ,

1 1 2 3 1+ 1 2

6 )を使用）してみると、下記式（73)〜（75)のように採用する変数によって計算され

3

る値が異なる。

[0130] [数 38]

次元

D = 325.0049920··· , D 5.9124806·

Ϊ0'5 ' 式 （7 3)

2次元

14.3484025·

式 (74)

3次元

£> —,-,- = 4347.7594268· D[3,2,-j = 46.8184447·

,10 0 4/ 式 (7 5)

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

[0132] [数 39] ' dxk，， ⁹ dxi , ··= ^ dxk' , ， , ¾

(N_+ 1)³. d^_a 3 。

N + 2 lム dxk dxi

式 (76)

1次元

[0133]

[0134]

[0135]

[0136]

[0137] ξ の 1変数の問題として、下記式 (82)を定義する。

1

[0138] [数 42]

= («i(fi) + _Q2( ) + 1)² { ）

. · ·式 （8 2)

[0139] 拡張された直交基底気泡関数要素の条件を満たす気泡関数が存在することは、実 数の集合 Rにおいて、 f(6 ) =0となる変数が存在することである。いま、下記式 (83

1

)〜（85)のような実数の集合 Rの部分集合 Sを考える。

[0140] [数 43]

1次元

1.8310 <S< 1.8315 式 （8 3)

2次元

5.2678 < S< 5.2683 式 （84)

3次元

1.8061 <S< 1.8066 式 (8 5)

[0141] ξ では、 1次元〜 3次元において、下記式（86)となる _c

[0142] [数 44]

6²( ) -4α(ξ!) θ0 式 （8 6)

[0143] 従って、 ex a は実数の範囲で値を持つ。また、 さ にて 1次元〜 3次元では、

1 2 1

それぞれ下記式 (87)および式 (88)が成り立つので、式（72)は有界な値をとり、式 ( 82)は実変数の関数となる。

[0144] [数 45] 1 1)+ 2( ) +1≠0 式 (8 7)

ΛΓ-2

Π (^+^+ )≠ 0 , m_d,n_d = 1,2,3

式 (8 8)

[0145] 式 (82)と部分集合 Sを使用すると下記式 (89)〜（91)を得ることが出来る。

[0146] [数 46] 次元

/(1.8310) < 0 ,バ 1.8315) > 0 , ^d f (fr) > 0

d c · · ·式 （ 8 9 )

2次元

/(5.2678) > 0 , /(5.2683) < 0 , く ο

· · ·式 （90)

3次元

(1.8061〉 > 0 , (1.8066) < 0 , ^d f (fr) < 0

· · ·式 （ 9 1 )

[0147] 任意の に対して、 1次元では、式（89)の第 1式が負、第 2式が正で、第 3式

lc

が正（単調増加）、 2次元、 3次元では、式（90)、 （91)の第 1式が正、第 2式が負で、 第 3式が負（単調減少）であるので、 1次元〜 3次元において、部分集合 Sの範囲に は、それぞれ f ( ξ ) =0の解が一意に存在する。

[0148] 具体的に、部分集合 Sの範囲内で、二分法を使用して f( 6 ) =0の解を求めると、

1

下記式（92)〜（94)の値が得られる。

[0149] [数 47]

1次元

^ = ξ = 1.8312889··· , ι = 20.4651928 ··■ , α₂ = -22.5761007· - -

• · ·式 （92)

2次元

ζ_χ =

, ¾ =-2.6201232··· , ο¾ = 3.1609628·--

• · ·式 （93)

3次元

^=^= 1.8063444--- =—17.3739109,·' , ₂ = 17.5493578···

• · ·式 （94)

[0150] 図 43、図 44に、式 (81)、（93)を用いた三角形気泡関数の形状 φ 、三角形気泡

Β

関数の形状 Φ Β ²を示す。

[0151] または、下記式（95)、（96)を使用して、 f(6 ) =0の解の存在を示すことができる

1

[0152] [数 48] 式 （95)

式 （96) d 9(ζι)

/(6)= J (6) = fl( ) =

[0153] 1次元〜 3次元において、 hを正の定数とし、任意の ξ に対して下記式（97)の lc

結果を得る。

[0154] [数 49] < < l 式 (97)

[0155] 平均値の定理から、任意の ξ , ξ に対して、下記式（98)が得られる。

la lb

[0156] [数 50]

, 0<h< l . . .式 （_{9 8)}

[0157] 従って、縮小写像の原理 (定理)より、式 (95)を反復的に逐次計算をさせる下記式

(99)にて得られる数列 ^m}は、部分集合 Sの範囲内で ί(ξ )=0を満足する一意

1 1

の解に収束する。

[0158] [数 51]

. . .式 （ 99 )

[0159] なお、拡張された直交基底気泡関数要素の存在を示すために、式 (80)の三つの 条件を満たす気泡関数として、三つの未知量を α 、 α 、 ξ とし、適当な気泡関数 φ

1 2 1

、 Φ 、 Φ を採用し力、他にも、下記式（100)の気泡関数 φ のように、未知量を _α

B1 Β2 Β3 Β

、 α 、 αとして、適当な気泡関数 φ 、 φ 、 φ 、 φ を採用すれば、式（80)の条

1 2 3 B1 Β2 Β3 Β4

件を満たす気泡関数が得られる。このように、拡張された直交基底気泡関数要素の 条件式を満たす気泡関数は、唯一ではなぐ幾つも存在するため、具体的な気泡関 数やその形状には殆ど意味がなぐ拡張された直交基底気泡関数要素で重要な意 味を成すのは、式 (80)の条件となる気泡関数の積分値である。 [0160] [数 52] ひ !_ +。2 +。3 + 1 式 （1 0 0 )

[0161] < 11.拡張された直交基底気泡関数要素の要素レベルの粘性項 (拡散項)の行列

>

通常の気泡関数、直交基底気泡関数要素の数値解析において式 (66)が必要に なる要素レベルの粘性項 (拡散項)の行列を下記式（101)〜（103)に示す。

[0162] [数 53]

次元

0

9Φ ΘΦ^Τ

ο

0 0 0

…式 （1 0 1 )

2

δ¾.9Φι Q

：

• · ·式 （1 0 2)

3次元

dxk

式 （1 0 3)

[0163] 式（101)〜（103)に式 (76)を代入して整理した、拡張された直交基底気泡関数 要素の要素レベルの粘性項 (拡散項)の行列を下記式（104)〜（106)に示す。

[0164] [数 54] 次元 0

• · ·式 （1 0 4 )

2次元

. 3Φ δΦ^τ. .

• · ·式 （1 0 5 )

3

式 0 以上、拡張された直交基底気泡関数要素の要素レベルの粘性項 (拡散項)の行列

(104)〜（106)を示した力 本発明に関する定式化やプログラムでは、行列（101) 〜（103)に式（76)を代入して整理を行わずに、そのまま行列（101)〜（103)と式（ 76)を使用することも、もちろん可能である。 [0166] <12.直交基底気泡関数要素と正規化気泡関数 >

式（1)は、下記式（107)、（108)のように変形することができる。

[0167] [数 55]

N+1

'式 （1 0 7)

ΛΓ+1

L Φαΐία + <j>S"S

'式 （1 0 8)

1 ( ¹ 、 ' (1,Φ_Β)η. ' N + 1 ,

'式 （1 0 9)

, _A_e ― , N + 2

φ^{3 =} Μ^ ^{φΒ =} ΝΤΪ' 式 （1 1 0)

[0168] 式（108)は、気泡関数自体の積分値が A (1次元では各要素 eの線、 2次元では各

e

要素 eの面積、 3次元では各要素 eの体積)となるような気泡関数を使用して、式変形 を行ったものであり、式（110)は正規化気泡関数と呼ばれている。直交基底気泡関 数要素にお 、て式（1)の表現形式の他に、式（ 107)の表現形式や式（108)の表現 形式を用いて、定式化やプログラミング、計算実行を行うことも可能である。また、式（ 1)、（107)、（108)の表現形式を用いて得られた有限要素方程式 (変分方程式)で は、静的縮約という操作によって要素ごとに重心点の自由度 u、 u '、 u 'を消去する

B B S

ことができるが、重心点の自由度を消去した有限要素方程式 (変分方程式)を用いて 、定式化やプログラミング、計算実行を行うことも可能である。直交基底気泡関数要 素の気泡関数として正規化気泡関数を使用した積分値は下記式 ( 111 )、（ 112)の ようになる。

[0169] [数 56] 式 （1 1 1)

式 （1 1 2)

[0170] さらに、拡張された直交基底気泡関数要素では、下記式（113)が得られる。

[0171] [数 57] 。'

(is•s · ·式)（1 1 3 )

[0172] < 13.直交基底気泡関数要素と安定化作用制御項 >

気泡関数要素の持つ数値不安定性を緩和させるために、下記式（114)の気泡関 数の粘性項 (安定ィ匕作用制御項)を用いて解析精度を向上させる方法がある。

[0173] [数 58] ν'^ί∑ _{n u}'。

'式（1 1 4 )

[0174] 式（114)の安定ィ匕作用制御パラメータ v，（一∞< V，<∞)は、粘性係数の次元（ 粘性係数に相当する変数を無次元化している場合は無次元）を持つパラメータであ り、この値を適切に決定することによって高精度解法を実現する。この方法は、従来、 気泡関数要素に通常の Galerkin法を使用して、得られた有限要素方程式に重心点 のみの粘性項 (安定化作用制御項)を付加して 、たが (松本純一,梅津剛,「気泡関数 を用いた非圧縮性流体の有限要素法解析に関する一考察」，日本応用数理学会平 成 8年度年会, 1996年, 46頁 47頁、松本純一,梅津剛,川原睦人,「線形型気泡関数 を用いた非圧縮性粘性流体解析と適応型有限要素法」，応用力学論文集 (土木学会 ) ,2卷, 1999年 8月，223頁 232頁参照）、現在では、重み関数に関してもう一つの気 泡関数（3レベル気泡関数)を採用して重心点のみの粘性項 (安定化作用制御項)を 導人して ヽる (J.Matsumoto ana M.Kawahara, Shape Iaentincation for Fluid— St ructure Interaction Problem Using Improved Bubble Element", International J ournal of Computational Fluid Dynamics, Vol.15, 2001, pp.33- 45、非特許文 献 3参照)。 3レベル気泡関数は、下記式（115)〜（117)を満たす気泡関数である。

[0175] [数 59] (Φα,^Β)Ω_ε = 0 , ひ = 1·' 式 （1 15)

^(±3レベル気泡関数 式 （1 16) 。 、

¾'お ι· ^Ω* 式 （1 1 7)

[0176] 式 (67)、（115)より下記式（118)を得る。

[0177] [数 60]

Ν+ϊ N+1

(l.^B>n. =《 *a,Vs)n. = L (*«> ΨΒΪα, = 0

ar=l a»l 式 U 18)

[0178] また、下記式（119)の関係を持っている 3レベル気泡関数の場合には、条件式（11 5)〜（： L 17)の式（ 115)は式（ 120)で代用できる。

[0179] [数 61]

N +1 -,ΨΒ)Ω,， a = l---N + l

'式 （1 1 9)

式 （1 20)

[0180] 式（119)の関係を持っている 3レベル気泡関数を仮定し、上記の条件式（120)、 ( 116)、 （117)を満たす具体的な 3レベル気泡関数として、下記式（121)を用いる。

[0181] [数 62] リ δΐψΒΙ + ^ϊψ Ι + ^3ψΒ3 + ψΒί

ΨΒ =

V ίΐ + ί₂ + ¾ + 1 式 （121)

[0182] 式中の νは、粘性係数 (拡散係数)である。式（121)を式（120)、（116)、 (117) に代入して整理すると、下記式（122)を得る。

[0183] [数 63] 一 1

δι Oil 12 13

21 23 h

"31 32 33 • · ·式 （ 1 22)

αιι = (^t PBi)n_e αΐ2 = (^ψΒ2)η_κ 13 == (l)^B3)n_e， α21 = ΦΒ,^Βΐ) _c 022 = {φΒτψΒ2)η_α 23 = <⁽^Β ^^3》Ω_<: d fc οχχ

d<f>Bヽ

32 =

033： >Ω_β，

も 3 =

^dxk ¹ dxi ^{h 1} dxk dxi

あるいは

— 1

11

³¹ (〈 dxk ' dxi 一（ "^

'W+l

₌ ί ,9ψΒ d<fB2 ,9φ_Β― d<i>jB、

\ dxk ' dxi ^lie dxk ' dxi ^∑^ 9xk dxi

_A ―

\ίι ^dxk 9χι

h = -{'ί,ψΒ4 αΛ ¹ , h =—{ΦΒ,ΨΒ (ϊ,Ας'

& ₃

dx_k dxi

[0184] 式（ 121)右辺の四つの気泡関数を具体的に下記式（ 123)のように定義する。

[0185] [数 64]

ΨΒΙ = ΦΒ > ΨΒ2 = Β » ΨΒ3 = , ΨΒ\ = Β . · ·式 （ 1 23 )

[0186] 式（123)、（59)より式（121)は式（119)を満足してレ、ることがわかる。式（123)の 変数を下記式（ 124)のように設定する。

[0187] [数 65] ¾ = 4 , 72 = 3 , τ/3 = 2 , "4 = -式 （1 24)

[0188] α ( ξ , ξ , ξ ) = α ( ξ , ξ , ξ )として、式（31)の変数を用いると、式（122)よ

1 1 2 3 1+ 1 2 3

り下記式（ 125)〜（ 127)を得る。

[0189] [数 66]

1次元

5_{1 =} -7.4533261··· , 5₂ = 13.9037038■· - , 5₃ = -7.4557822 ·· · . . '式 （1 2 5)

2次元

δι = -12.8855927-■ . , = 21.8201295· · · , 5₃ = -9.9378406■- - . . .式 （ 1 2 6 )

3次元

¾ = 73.7388138-·· , ¾ = -102.3508563 ·· · , ¾ = 27.6071956… . . .式 （1 2 7) [0190] α ( ξ , ξ , ξ ) = α ( ξ , ξ , ξ )として、式（36)の変数を用いると、式（122)よ

1 1 2 3 1+ 1 2 3

り下記式（ 128)〜（ 130)を得る。

[0191] [数 67]

1次元

δι = -5.2897450·· · , ¾ = 10.5097066-· · , = -6.2084329· -- . . .式 （ 1 2 8 )

2次元

5ι= -9.2351821-■· , ¾ = 16,2678560… ， =—8.0666407… . · '式 （1 2 9)

3次元

δι = -24.4550047·-· , ¾ = 37.6758731 ·· · ， ¾ =一 14.3507923··. . . '式 （1 3 0)

[0192] α ( ξ , ξ , ξ ) = α ( ξ , ξ , ξ )として、式（92)〜（94)、式 (81)の変数を用い

1 1 2 3 1+ 1 2 3

ると、式（122)より下記式（131)〜（133)を得る。

[0193] [数 68]

1次元

5ι = -6.8577197· · - , ¾ = 12.6283403- -- , δ₃ = -6.8566234· -- . . . （ 1ゥ 1 )

2次元

¾ = -7.0364016-·- , ¾ = 13.3033720··· , ¾ =—7.1674625··. . . .す （1 Q 2)

3次元

δ = -12.0454679··- , = 20.3290654·-· , δ₃ = -9.2229704··· . . ' ( Q ₃ ) [0194] 図 45、図 46に、式（132)を用いた三角形 3レベル気泡関数の形状（ v ' = vとした 場合)、式 (93)、 （81)を用いた三角形気泡関数と式（132)を用いた三角形 3レベル 気泡関数の積の形状（ V ' = Vとした場合)を示す。以上により、直交基底気泡関数 要素を用いた場合にも、式（120) (あるいは式（115) )、 (116)、 (117)を満たす 3レ ベル気泡関数は存在し、安定化作用制御項（114)を導入できる。

[0195] なお、 3レベル気泡関数の存在を示すため、式（120)、 （116)、 （117)の三つの条 件を満足する三つの未知量 δ 、 δ 、 δ と四つの適当な気泡関数を持った 3レベル

1 2 3

気泡関数（121)を採用したが、この他にも、下記式（134)のように、安定化作用制御 ノ メータ V 'と同じ次元を持つパラメータを導入し、未知量を δ、 δ として、式（13

1 2

4)の右辺に適当な三つの気泡関数を用いることによって、式（120)、 （116)の条件 を満たす 3レベル気泡関数を求めることができる。そして、下記式（135)のように安定 化作用制御パラメータ V 'を定義すれば、式（117)も導出でき、式（120) (あるいは 式（115) )、 （116)、 （117)となる 3レベル気泡関数の存在を示すことが出来る。この ように、 3レベル気泡関数の条件式を満たす気泡関数は、唯一ではなぐ幾つも存在 するため、具体的な気泡関数やその形状には殆ど意味がなぐ 3レベル気泡関数で 重要な意味を成すのは、安定化作用制御項を導くために必要な式（120) (あるいは 式（115) )、 (116)、 (117)の条件である。

[0196] [数 69] ρ ν_ψ διψΒΙ + δ₂ψΒ2 + ΨΒ3

ν δ₁ + δ₂ + 1 . . .式 （1 3 4 )

^(-οο < ν_ψ < οο) は と同じ次元を持つパラメータ

〈¾ ' 》 Ω' · · '式 （1 3 5 )

[0197] く 14.直交基底気泡関数要素と正規ィ匕気泡関数による安定ィ匕作用制御項〉

正規化気泡関数を使用した場合にも安定化作用制御項は導入することが出来、下 記式（136)によって表される。

[0198] [数 70] ' ひ 式 （1 3 6 )

[0199] 直交基底気泡関数要素で使用する気泡関数として正規化気泡関数を用いた場合 、 3レベル正規ィ匕気泡関数は下記式（137)と定めることができ、その積分値は下記 式（ 138)〜（ 142)のように記述できる。

[0200] [数 71] ん

Ψ3 ■ψΒ

式 （1 3 7 )

A,

(Φα, <pB) _c = 0 , = l---N+l

〈1 〉 Ω. 式 （1 3 8)

A?

/1 \2 Β ΨΒ

く¹ 〉 )ひ, = 0

式 （1 3 9 ) dx_k ' dx dx_k dxi'"' 式 （1 4 0 )

N +l 〈1 〉 n_eN + l

式 （1 4 1 )

式 ( 1 4 2)

[0201] 上記式（138)〜（142)より、式（115)〜（117)、 （119)、 （120)を満足する 3レべ ル気泡関数が存在すれば、式（138)〜（142)も自動的に満たすので 3レベル正規 化気泡関数も同時に存在することがわかる。

[0202] <15.直交基底気泡関数要素と気泡関数の安定化作用行列 >

式（107)の表現形式を用いた気泡関数要素の重心点の自由度 u 'は、静的縮約と

Β

いう操作によって要素ごとに消去できる。その結果得られる有限要素方程式は、安定 化有限要素法で導かれる有限要素方程式と密接な関係がある (非特許文献 2参照） 。静的縮約によって得られる気泡関数の安定ィ匕作用行列は、一般的に下記式（143 )の形式によって表される _c

[0203] [数 72] r_eS = Κ{ν )_Β

ん 式 （1 4 3)

K

K(y')_B = K(v)_Bmn― m,n= 1 · · · L

K

式 （1 44)

[0204] Lは、求めるべき未知解析物理量の数である。式（144)は、式 (69)を満たす気泡 関数を使用した場合、式 (70)、（65)、（66)の気泡関数の積分値から成る行列であ る。安定ィ匕作用制御パラメータ V 'は、下記式（145)の関係が満たされるように決定 する。

[0205] [数 73] 式 （1 4 5)

TeR = TeRmn =

子 eRLl .. · T_eRL 式 （1 4 6) UB1 式 （1 4 7)

[0206] 式（146)は解析精度が向上することが期待できる安定ィ匕行列である。式（145)の 関係を満たすような安定ィ匕作用制御パラメータ V 'を求めるために、下記式（148)の 評価関数を最小にするような安定ィ匕作用制御パラメータ V 'を決定する問題を考える

[0207] [数 74] 式 （1 4 8) Ι,

式 （1 4 9) [0208] 式（148)を最小にするような安定ィ匕作用制御パラメータ v 'を求めるために、式（14 8)の停留条件である下記関係式（150)を用いて、安定ィ匕作用制御パラメータ を 求める。レ J ,

[0209] [数 75] τ

a v( )

ν{ν')_Β = 0

d v 式 （1 5 0 )

[0210] 3レベル気泡関数は重み関数の中で導入しているので、求めるべき未知解析物理 量の偏微分方程式ごとに異なった 3レベル気泡関数を使用することができる。従って 、安定ィ匕作用制御パラメータ V 'は、下記式（151)の 3レベル気泡関数を用いると、 求めるべき未知解析物理量の数 Lだけ、安定化作用制御パラメータ V 'を定義する m

ことができる。

[0211] [数 76]

式 （1 5 1 )

[0212] 式（151)の係数 δ 、 δ 、 δ は、式（122)によって決定することが出来、式（151)

1 2 3

は式（120)、 （116)を満たす。また、式（117)の代わりに下記式（152)の安定ィ匕作 用制御項を得る。

[0213] [数 77] dxk ' dxi

式 （1 5 2 )

[0214] 式（152)の安定ィ匕作用制御項を用いた場合、式（144)は下記式（153)のようにな る。

[0215] [数 78]

^K{^vm' )B = K{v_m)_Bmn

K{ _L' )BLI · · · K{v_L)_BLL 式 （1 5 3 ) [0216] 式（145)の関係を満たすような安定ィ匕作用制御パラメータ v 'を求めるために、下 m

記式（154)の評価関数を最小にするような安定ィヒ作用制御パラメータ V 'を決定す m る問題を考える。

[0217] [数 79] 式 （1 5 4 )

"(レ _m

式 （1 5 5 )

[0218] 式（154)を最小にするような安定ィ匕作用制御パラメータ v 'を求めるために、式（1 m

54)の停留条件である下記関係式（156)を用いて、安定化作用制御パラメータ V ' m を求める。

[0219] [数 80]

τ

d j_ 9 v{ _m)_B

= o

式 （1 5 6 )

[0220] 式（155)はベクトルであり、ベクトルの成分を下記式（157)のように表現する。

[0221] [数 81]

· · · v{v_L) BL

式 （1 5 7 )

[0222] 式（157)をみるとわ力るように、安定化作用制御パラメータ V 'は各ベクトル成分で m

独立に定義されているので、式（145)の関係を満たすような安定ィ匕作用制御パラメ ータ V 'を決定する問題として、式（157)の成分 mごとに、下記式（158)の評価関 m

数を用いた最小化問題を考えることができる。

[0223] [数 82]

式 ( 1 5 8 )

[0224] 式（158)を最小にするような安定化作用制御パラメータ v 'を求めるために、式（1 m

58)の停留条件である下記関係式（159)を用いて、安定化作用制御パラメータ V ' を求める。

[0225] [数 83]

式 （1 59)

[0226] 式（157)は、安定化作用制御パラメータ V 'が各ベクトル成分で独立に定義され m

ているので、式（154)、 （158)の停留条件力 導かれる安定ィ匕作用制御パラメータ

V 'を求める式（156)、 (159)は、結果的に同じ式になる。

m

[0227] 他の方法として、下記式（160)の関係を満たすように、安定化作用制御パラメータ

V，を決定する問題も考えられる。

[0228] [数 84]

^TeB^UB T_eR^UB 式 （1 60)

[0229] 式（160)の関係を満たすような安定ィ匕作用制御パラメータ V 'を求めるために、下 記式（161)の評価関数を最小にするような安定ィヒ作用制御パラメータ V 'を決定す る問題を考える。

[0230] [数 85] 式 （161)

式 (1 62)

[0231] 式（161)を最小にするような安定ィ匕作用制御パラメータ v 'を求めるために、式（16 1)の停留条件である下記関係式（163)を用いて、安定ィ匕作用制御パラメータ を 求める。

[0232] [数 86] d

W V )B =0

\~dV

式 （163) 安定ィ匕作用制御パラメータ v 'は、式（151)の 3レベル気泡関数を用いると、式（15 2)の安定化作用制御項が導出され、求めるべき未知解析物理量の数 Lだけ、安定 化作用制御パラメータ v 'を定義することができる。式（160)の関係を満たすような

m

安定化作用制御パラメータ V 'を求めるために、下記式（164)の評価関数を最小に

m

するよう -な安定ィ匕作用制御パラメータ V 'を決定する問題を考える。

m

[0234] [数 87] 式 （1 6 4 ) ' 、 ' (Ι, Β)² .., ,

式 （1 6 5 )

[0235] 式（164)を最小にするような安定ィ匕作用制御パラメータ v 'を求めるために、式（1

m

64)の停留条件である下記関係式（166)を用いて、安定化作用制御パラメータ V ' m を求める。

[0236] [数 88]

T

d j d w _m)s

式 （1 6 6 )

[0237] 式（165)はベクトルであり、ベクトルの成分を下記式（167)のように表現する。

[0238] [数 89] ) B = 1

式 （1 6 7 )

[0239] 式（167)をみるとわ力るように、安定化作用制御パラメータ V 'は各ベクトル成分で

m

独立に定義されているので、式（160)の関係を満たすような安定ィ匕作用制御パラメ ータ V 'を決定する問題として、式（167)の成分 mごとに、下記式（168)の評価関 m

数を用いた最小化問題を考えることができる。

[0240] [数 90]

式 （1 6 8 )

[0241] 式（168)を最小にするような安定ィ匕作用制御パラメータ V 'を求めるために、式（1

m

68)の停留条件である下記関係式（169)を用いて、安定化作用制御パラメータ V ' を求める。

[0242] [数 91] d j d w{v_m)_Bm

w{v_m)_Bm = 0

d v_m

式 （1 6 9 )

[0243] 式（167)は、安定ィ匕作用制御パラメータ v 'が各ベクトル成分で独立に定義され

m

ているので、式（164)、 （168)の停留条件力 導かれる安定ィ匕作用制御パラメータ

V 'を求める式（166)、 （169)は、結果的に同じ式になる。

m

[0244] < 16.直交基底気泡関数要素と正規化気泡関数の安定化作用行列 >

式（108)の表現形式を用いた気泡関数要素の重心点の自由度 u 'は、静的縮約

S

という操作によって要素ごとに消去できる。静的縮約によって得られる正規化気泡関 数の安定ィ匕作用行列は、一般的に下記式（170)の形式によって表される。

[0245] [数 92]

. . .式 （ 1 7 0 ) ' 〉 Ω, · · ·式 （1 7 1 )

[0246] 式（170)の第二項は、式（141)の関係を満たす気泡関数を使用した場合、式（11 1)〜（113)各第一項の正規ィ匕気泡関数の積分値力 成る行列であり、式 (70)、 (6 5)、 （66)で表した場合には第三項のように記述できる。安定化作用制御パラメータ

V 'は、下記式（172)の関係が満たされるように決定する。

[0247] [数 93]

TeSU_s f¾ T_&RU₃ · · ·式 ( 1 7 2 )

· . ·式 （1 7 3 )

[0248] 式（172)の関係を満たすような安定ィ匕作用制御パラメータ V 'を求めるために、下 記式（174)の評価関数を最小にするような安定ィヒ作用制御パラメータ V 'を決定す る問題を考える。 [0249] [数 94] 式 （1 74)

_eRu_s - _A ^us

式 （1 75)

[0250] 式（174)を最小にするような安定ィ匕作用制御パラメータ v 'を求めるために、式（17 4)の停留条件である下記関係式（176)を用いて、安定ィ匕作用制御パラメータ を 求める。

[0251] [数 95]

式 （1 76)

[0252] 3レベル正規ィ匕気泡関数は重み関数の中で導入しているので、求めるべき未知解 析物理量の偏微分方程式ごとに異なった 3レベル正規化気泡関数を使用することが できる。従って、安定ィ匕作用制御パラメータ V 'は、下記式（177)の 3レベル正規ィ匕 気泡関数を用いると、求めるべき未知解析物理量の数 Lだけ、安定化作用制御パラ メータ V 'を定義することができる。

m

[0253] [数 96]

Ψβ =

式 （1 77)

[0254] 式（177)の係数 δ ヽ δ 、 δ は、式（122)によって決定することが出来、式（177)

1 2 3

は式（142)、 （139)を満たす。また、式（140)の代わりに下記式（178)の安定ィ匕作 用制御項を得る。

[0255] [数 97]

^=Im(—

ax oxi .-s -)n_e

dxk dxi 式 （ 8) 式（ 178)の安定化作用制御項を用 、た場合、式（ 171)は下記式（ 179)のようにな 5 c

[0257] [数 98] k{v_m)_s = k{u_m)_B

式 （1 7 9 )

[0258] 式（179)の第一項は、式（111)〜（113)各第一項の正規ィ匕気泡関数の積分値か ら成る行列であり、式（70)、 (65) (66)で表した場合には第二項のように記述できる 。式（172)の関係を満たすような安定ィ匕作用制御パラメータ V 'を求めるために、下 m

記式（180)の評価関数を最小にするような安定ィヒ作用制御パラメータ V 'を決定す m る問題を考える。

[0259] [数 99] 式 （1 8 0 ) （ '、 - ' (Ι, Φε)² '■

式 （1 8 1 )

[0260] 式（180)を最小にするような安定ィ匕作用制御パラメータ v 'を求めるために、式（1 m

80)の停留条件である下記関係式（182)を用いて、安定化作用制御パラメータ V ' m を求める。

[0261] [数 100] d J

"iyj' _s

d = o

式 （1 8 2 )

[0262] 式（181)はベクトルであり、ベクトルの成分を下記式（183)のように表現する。

[0263] [数 101]

"( si · · ' V{ _L)SL

式 （1 8 3 ) 式（183)をみるとわ力るように、安定化作用制御パラメータ V 'は各ベクトル成分で m

独立に定義されているので、式（172)の関係を満たすような安定ィ匕作用制御パラメ ータ V 'を決定する問題として、式（183)の成分 mごとに、下記式（184)の評価関 数を用いて最小化問題を考えることができる。

[数 102]

式 （1 8 4 )

[0266] 式（184)を最小にするような安定ィ匕作用制御パラメータ v 'を求めるために、式（1

m

84)の停留条件である下記関係式（185)を用いて、安定化作用制御パラメータ V ' m を求める。

[0267] [数 103]

式 （1 8 5 )

[0268] 式（183)は、安定ィ匕作用制御パラメータ v 'が各ベクトル成分で独立に定義され

m

ているので、式（180)、 （184)の停留条件力も導かれる安定ィ匕作用制御パラメータ

V 'を求める式（182)、（185)は、結果的に同じ式になる。

m

[0269] 他の方法として、下記式（186)の関係を満たすように、安定化作用制御パラメータ

V，を決定する問題も考えられる。

[0270] [数 104] s^us∞干； k^us 式 （1 8 6 )

[0271] 式（186)の関係を満たすような安定ィ匕作用制御パラメータ V 'を求めるために、下 記式（187)の評価関数を最小にするような安定ィヒ作用制御パラメータ V 'を決定す る問題を考える。

[0272] [数 105] 式 （1 8 7 )

式 （1 8 8 )

[0273] 式（187)を最小にするような安定化作用制御パラメータ v 'を求めるために、式（18 7)の停留条件である下記関係式（189)を用いて、安定化作用制御パラメータ V 'を 求める。

[0274] [数 106ノ J]

式 （1 8 9 )

[0275] 安定ィ匕作用制御パラメータ v 'は、式（177)の 3レベル正規ィ匕気泡関数を用いると 、式（178)の安定ィヒ作用制御項が導出され、求めるべき未知解析物理量の数 Lだけ 、安定化作用制御パラメータ V 'を定義することができる。式（186)の関係を満たす m

ような安定ィ匕作用制御パラメータ V 'を求めるために、下記式（190)の評価関数を m

最小にするような安定ィ匕作用制御パラメータ V 'を決定する問題を考える。

m

[0276] [数 107] 式 （1 9 0 )

式 （1 9 1 )

[0277] 式（190)を最小にするような安定ィ匕作用制御パラメータ V 'を求めるために、式（1 m

90)の停留条件である下記関係式（192)を用いて、安定化作用制御パラメータ V ' m を求める。

[0278] [数 108]

τ

d J— 9 w\v_m)_S

式 （1 9 2 )

[0279] 式（191)はベクトルであり、ベクトルの成分を下記式（193)のように表現する。

[0280] [数 109] si L)SL

式 （1 9 3 ) 式（193)をみるとわ力るように、安定化作用制御パラメータ V 'は各ベクトル成分で m

独立に定義されているので、式（186)の関係を満たすような安定ィ匕作用制御パラメ ータ V 'を決定する問題として、式（193)の成分 mごとに、下記式（194)の評価関 m

数を用いた最小化問題を考えることができる。

[0282] [数 110]

L

式 4 )

[0283] 式（194)を最小にするような安定ィ匕作用制御パラメータ v 'を求めるために、式（1 m

94)の停留条件である下記関係式（195)を用いて、安定化作用制御パラメータ V '

m を求める。

[0284] [数 111]

· · ·式 （ 1 9 5 )

[0285] 式（193)は、安定ィ匕作用制御パラメータ v 'が各ベクトル成分で独立に定義され m

ているので、式（190)、 （194)の停留条件力も導かれる安定ィ匕作用制御パラメータ V 'を求める式（192)、 （195)は、結果的に同じ式になる。

m

[0286] < 17.気泡関数と正規ィ匕気泡関数の安定ィ匕作用制御パラメータの関係 >

直交基底気泡関数要素で、式 (1)、 （107)の表現形式を採用した場合には、式 (7 0)、 (65)、 (66)の積分を使用して、式（150)、 （156)、 (159)、 (163)、 （166)、 (1 69)のいずれかを採用し、式（108)の表現形式を採用した場合には、式（111)〜（1 13)各第一項の積分を使用して、式（176)、 (182)、 (185)、 (189)、 (192)、 (195 )のいずれかを採用すれば安定ィ匕作用制御パラメータを計算することが出来る。いま 、式（176)、 （182)、 （185)、 （189)、 （192)、 （195)を、式（70)、（65)、 （66)の積 分を使用して記述すると下記式（196)〜（201)のようになる。

[0287] [数 112]

式 （196) ^_m)s A? d v{v_m' )_B

Ω« L 式 （1 97)

「3 v{v_m' )_Sm d v( _m' )_Bm

〈1， 〉²- d v'

式 （1 98) 式 (1 99)

式 (200)

Al d w{v_m' )B_m

w

d v' 〈1 ( o

Ω. 式 (201)

[0288] 式（70)、（65)、（66)の積分を使用した式（196) （201)の第二項より、安定ィ匕 作用制御パラメータを求めるための関係式は、式（150)、 (156)、 (159)、 (163)、 ( 166) (169)の安定ィ匕作用制御パラメータを求めるための関係式と同じに（等しく） なる。従って、式（150)、（156)、（159)、（163)、（166)、（169)より式（70) (65) （66)の気泡関数の積分を使用して求めた安定ィ匕作用制御パラメータと、式（176) （182)、（185)、（189)、（192)、（195)より式（111)〜（113)各第一項の正規ィ匕 気泡関数の積分を使用して求めた安定ィ匕作用制御パラメータは、同じ値 (等しい値） になることがわかる。

[0289] 上記より、 1)気泡関数、気泡関数自体の積分値が Aとなる正規化気泡関数、 2)気 e

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

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

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

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

産業上の利用可能性

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

Claims

請求の範囲

[1] 解析対象の各要素の整合質量行列を取得する取得工程と、

前記解析対象の各要素の気泡関数に基づいて、前記取得工程によって取得され た各要素の整合質量行列を対角化した、各要素の対角質量行列を生成する生成ェ 程と、

前記解析対象の既知解析物理量と、前記生成工程によって生成された各要素の 対角質量行列と、に基づいて、前記解析対象の挙動を解析する解析工程と、 を含んだことを特徴とする直交基底気泡関数要素数値解析方法。

[2] 前記生成工程は、前記各要素における前記気泡関数の積分値を、前記各要素の 整合質量行列に代入することにより、前記各要素の対角質量行列を生成することを 特徴とする請求項 1に記載の直交基底気泡関数要素数値解析方法。

[3] 前記生成工程によって生成された各要素の対角質量行列に基づいて、前記解析 対象全域の対角質量行列を算出する第 1の算出工程と、

前記第 1の算出工程によって算出された、前記解析対象全域の対角質量行列の逆 行列を算出する第 2の算出工程と、を含み、

前記解析工程は、

前記解析対象の既知解析物理量と、前記解析対象全域の対角質量行列と、前記 第 2の算出工程によって算出された逆行列と、に基づいて、前記解析対象の挙動を 解析することを特徴とする請求項 1または 2に記載の直交基底気泡関数要素数値解 析方法。

[4] 前記解析対象の各要素の整合質量行列を取得する取得工程と、

前記解析対象の各要素の気泡関数に基づいて、前記取得工程によって取得され た各要素の整合質量行列を対角化した、各要素の対角質量行列を生成する生成ェ 程と、

前記解析対象の既知解析物理量と、前記生成工程によって生成された各要素の 対角質量行列と、に基づいて、前記解析対象の挙動を解析する解析工程と、 をコンピュータに実行させることを特徴とする直交基底気泡関数要素数値解析プロ グラム。

[5] 前記生成工程は、前記各要素における前記気泡関数の積分値を、前記各要素の 整合質量行列に代入することにより、前記各要素の対角質量行列を生成することを 特徴とする請求項 4に記載の直交基底気泡関数要素数値解析プログラム。

[6] 前記生成工程によって生成された各要素の対角質量行列に基づいて、前記解析 対象全域の対角質量行列を算出する第 1の算出工程と、

前記第 1の算出工程によって算出された、前記解析対象全域の対角質量行列の逆 行列を算出する第 2の算出工程と、をコンピュータに実行させ、

前記解析工程は、

前記解析対象の既知解析物理量と、前記解析対象全域の対角質量行列と、前記 第 2の算出工程によって算出された逆行列と、に基づいて、前記解析対象の挙動を 解析することを特徴とする請求項 4または 5に記載の直交基底気泡関数要素数値解 析プログラム。

[7] 前記解析対象の各要素の整合質量行列を取得する取得手段と、

前記解析対象の各要素の気泡関数に基づいて、前記取得手段によって取得され た各要素の整合質量行列を対角化した、各要素の対角質量行列を生成する生成手 段と、

前記解析対象の既知解析物理量と、前記生成手段によって生成された各要素の 対角質量行列と、に基づいて、前記解析対象の挙動を解析する解析手段と、 を備えることを特徴とする直交基底気泡関数要素数値解析装置。

[8] 前記生成手段は、前記各要素における前記気泡関数の積分値を、前記各要素の 整合質量行列に代入することにより、前記各要素の対角質量行列を生成することを 特徴とする請求項 7に記載の直交基底気泡関数要素数値解析装置。

[9] 前記生成手段によって生成された各要素の対角質量行列に基づいて、前記解析 対象全域の対角質量行列を算出する第 1の算出手段と、

前記第 1の算出手段によって算出された、前記解析対象全域の対角質量行列の逆 行列を算出する第 2の算出手段と、を備え、

前記解析手段は、

前記解析対象の既知解析物理量と、前記解析対象全域の対角質量行列と、前記 第 2の算出手段によって算出された逆行列と、に基づいて、前記解析対象の挙動を 解析することを特徴とする請求項 7または 8に記載の直交基底気泡関数要素数値解 析装置。