TOP > 国内特許検索 > 直交基底気泡関数要素数値解析方法、直交基底気泡関数要素数値解析プログラムおよび直交基底気泡関数要素数値解析装置 > 明細書

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

発行国 日本国特許庁(JP)
公報種別 特許公報(B2)
特許番号 特許第4729767号 (P4729767)
登録日 平成23年4月28日(2011.4.28)
発行日 平成23年7月20日(2011.7.20)
発明の名称または考案の名称 直交基底気泡関数要素数値解析方法、直交基底気泡関数要素数値解析プログラムおよび直交基底気泡関数要素数値解析装置
国際特許分類 G06F  17/50        (2006.01)
FI G06F 17/50 612H
請求項の数または発明の数 9
全頁数 60
出願番号 特願2006-547876 (P2006-547876)
出願日 平成17年11月25日(2005.11.25)
国際出願番号 PCT/JP2005/021727
国際公開番号 WO2006/057359
国際公開日 平成18年6月1日(2006.6.1)
優先権出願番号 2004343213
2005071239
優先日 平成16年11月26日(2004.11.26)
平成17年3月14日(2005.3.14)
優先権主張国 日本国(JP)
日本国(JP)
審査請求日 平成19年9月21日(2007.9.21)
特許権者または実用新案権者 【識別番号】503360115
【氏名又は名称】独立行政法人科学技術振興機構
【識別番号】301021533
【氏名又は名称】独立行政法人産業技術総合研究所
発明者または考案者 【氏名】松本 純一
個別代理人の代理人 【識別番号】100104190、【弁理士】、【氏名又は名称】酒井 昭徳
審査官 【審査官】鹿野 博嗣
参考文献・文献 特開2003-220808(JP,A)
特開2002-328955(JP,A)
国際公開第2004/088563(WO,A1)
特開2005-115934(JP,A)
特開2006-099549(JP,A)
調査した分野 G06F 17/50
特許請求の範囲 【請求項1】
生成手段と解析手段とを備えたコンピュータが、
前記生成手段が、解析対象範囲に分割形成されたメッシュごとに、当該メッシュの解析対象となる有限要素方程式内の整合質量行列に対して、気泡関数要素の基底が直交する条件を満たす気泡関数を求め、前記整合質量行列に、前記気泡関数の積分値を代入することによって、前記整合質量行列を対角化した対角質量行列を生成する生成工程と、
前記解析手段が、前記解析対象となる物理量の物性を表す既知物性値および前記解析対象となる物理量の解析条件を表す境界解析物理量からなる既知解析物理量を、前記生成工程によって生成されたメッシュごとの対角質量行列に代入することによって、前記解析対象となる物理量の変化を解析する解析工程と、
を実行することを特徴とする直交基底気泡関数要素数値解析方法。
【請求項2】
前記解析手段が、前記既知物性値、前記境界解析物理量および前記解析対象となる物理量の解析開始時の値を表す初期解析物理量からなる既知解析物理量を、前記生成工程によって生成されたメッシュごとに対角質量行列に代入することによって、前記解析対象となる物理量の微少時間における変化を解析することを特徴とする請求項1に記載の直交基底気泡関数要素数値解析方法。
【請求項3】
前記生成工程は、
前記整合質量行列に対して、気泡関数要素の基底が直交する下記式(1)を満たす気泡関数を求め、前記整合質量行列に、前記気泡関数の積分値を代入することによって、前記整合質量行列を対角化した対角質量行列を生成することを特徴とする請求項1または2に記載の直交基底気泡関数要素数値解析方法。
【数1】
JP0004729767B2_000114t.gif

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

【請求項9】
前記生成手段によって生成された前記解析対象範囲のメッシュごとの対角質量行列の総和を求めることによって、前記解析対象範囲全域の対角質量行列を算出する第1の算出手段と、
前記第1の算出手段によって算出された前記解析対象範囲全域の対角質量行列の逆行列を算出する第2の算出手段と、を備え、
前記解析手段は、
前記解析対象範囲全域についての前記既知解析物理量と前記第2の算出手段によって算出された前記対角質量行列の逆行列とを乗算することによって、前記解析対象となる物理量の変化を解析することを特徴とする請求項6または8に記載の直交基底気泡関数要素数値解析装置。
発明の詳細な説明 【技術分野】
【0001】
この発明は、気泡関数要素を用いた有限要素法による解析(有限要素解析)について、計算効率の良い対角項のみとなる質量行列を用いて、信頼性の高い数値シミュレーションをおこなうための直交基底気泡関数要素数値解析方法、直交基底気泡関数要素数値解析プログラム、および直交基底気泡関数要素数値解析装置に関する。
【背景技術】
【0002】
従来の気泡関数要素について説明する。図47は、従来の2次元の気泡関数要素を示す説明図であり、図48は、従来の3次元の気泡関数要素を示す説明図である。図47および図48のように、三角形(四面体)要素を用いた気泡関数要素は、各要素において三角形(四面体)を形成する3(4)点と重心点の4(5)つの節点を用いて、アイソパラメトリック座標系[r,s]({r,s,t})で下記式(1)のように表される(たとえば、下記非特許文献1、非特許文献2、非特許文献3参照)。
【0003】
【数1】
JP0004729767B2_000002t.gif

【0004】
式(1)のΦα,φBは、気泡関数要素の形状関数、uα,uBは三角形(四面体)の各頂点の値(解析物理量)、重心点の値(解析物理量)、Nは空間次元数を示している。ベクトル形式で記述すれば、形状関数は下記式(3)~式(6)になる。
【0005】
【数2】
JP0004729767B2_000003t.gif

【0006】
式(2)のΨαは、2次元および3次元の一次要素を用いた形状関数であり、下記式(7)、(8)であらわされる。
【0007】
【数3】
JP0004729767B2_000004t.gif

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

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

【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.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 (ジェー・シー・シモ、エフ・アルメロ、シー・エー・テイラー著 「ステイブル アンド タイム-ディシペイティブ ファイナイト エレメント メソーズ フォー ザ インコンプレッシブル ナビエ-ストークス イクエイションズ イン アドベクション ドミネーティド フロウズ」 インターナショナル ジャーナル フォー ニューメリカル メソーズ イン エンジニアリング 38巻 1995年 1475頁-1506頁)
【非特許文献3】松本純一,「気泡関数を用いた非圧縮性粘性流れ解析のための2レベル-3レベル有限要素法」,応用力学論文集(土木学会),7巻,2004年8月,339頁-346頁
【非特許文献4】畑中勝守,「多段階有限要素法による非圧縮粘性流体の順・逆解析に関する計算力学的研究」,中央大学博士論文,1993年3月
【発明の開示】
【発明が解決しようとする課題】
【0014】
ここで、上述した式(10)~式(13)に示したように、4段解法を使用して、この既知解析物理量unから未知解析物理量un+1を求めるためには、質量行列の逆行列が必要になる。図49は、従来のRotating Cone問題の解析に用いる解析モデルを示す説明図であり、図50は、従来のRotating Cone問題の解析に用いる初期条件の等高線を示す説明図であり、図51は、従来の整合質量行列によって計算されたRotating Cone問題の解析結果の等高線を示す説明図である(非特許文献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)を満たす気泡関数を用いる。
【数6】
JP0004729767B2_000007t.gif

【0023】
上記式(15-1)~式(15-3)は、気泡関数要素の定式化をおこなうために必要となる積分であり、<・,・>Ωeは要素領域Ωeの積分、Aeは要素領域Ωeの面積(体積)を示している。これにより、各要素の整合質量行列内の成分を足し合わせることによって近似された集中質量行列を用いずに、気泡関数要素の基底(形状関数)が直交する条件を導入して、各要素の高精度な対角質量行列を用いて解析することができる。また、解析対象全域の対角質量行列やその逆行列も簡単に算出することができ、解析対象の挙動を、未知の解析対象物理量を用いて高精度に解析することができる。
【発明の効果】
【0024】
本発明にかかる直交基底気泡関数要素数値解析方法、直交基底気泡関数要素数値解析プログラム、および直交基底気泡関数要素数値解析装置によれば、簡単かつ高精度に解析処理をおこなうことにより、解析対象の挙動解析の信頼性の向上を図ることができるという効果を奏する。また、記憶容量の低減化および解析時間の短縮化など、効率的な解析手法(直交基底気泡関数要素数値解析方法)を実現することができるという効果を奏する。
【図面の簡単な説明】
【0025】
【図1】図1は、この発明の実施の形態にかかる数値解析装置のハードウェア構成を示すブロック図である。
【図2】図2は、この発明の実施の形態にかかる数値解析装置の機能構成を示すブロック図である。
【図3】図3は、この発明の実施の形態にかかる数値解析装置のデータ格納部が格納している解析物理量情報テーブルを示す説明図である。
【図4】図4は、この発明の実施の形態にかかる数値解析装置における数値解析処理手順を示すフローチャートである。
【図5】図5は、2次元の気泡関数要素を示す説明図である。
【図6】図6は、3次元の気泡関数要素を示す説明図である。
【図7】図7は、直交基底を形成する三角形気泡関数の形状φBを示す説明図である。
【図8】図8は、直交基底を形成する三角形気泡関数の形状φB2を示す説明図である。
【図9】図9は、直交基底を形成する三角形気泡関数の形状φBを示す説明図である。
【図10】図10は、直交基底を形成する三角形気泡関数の形状φB2を示す説明図である。
【図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問題の解析における解析モデルのメッシュ(節点数:438413、要素数: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)を用いた三角形気泡関数の形状φB2を示す説明図である。
【図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と、I/F(インターフェース)109と、キーボード110と、マウス111と、スキャナ112と、プリンタ113と、を備えている。また、各構成部は、バス100によってそれぞれ接続されている。
【0029】
ここで、CPU101は、数値解析装置の全体の制御を司る。ROM102は、ブートプログラムなどのプログラムを記憶している。RAM103は、CPU101のワークエリアとして使用される。HDD104は、CPU101の制御にしたがってHD105に対するデータのリード/ライトを制御する。HD105は、HDD104の制御で書き込まれたデータを記憶する。
【0030】
FDD106は、CPU101の制御にしたがってFD107に対するデータのリード/ライトを制御する。FD107は、FDD106の制御で書き込まれたデータを記憶したり、FD107に記憶されたデータを数値解析装置に読み取らせたりする。着脱可能な記録媒体として、FD107のほか、CD-ROM(CD-R、CD-RW)、MO、DVD(Digital Versatile Disk)、メモリーカードなどであってもよい。ディスプレイ108は、カーソル、アイコンあるいはツールボックスをはじめ、文書、画像、機能情報などのデータを表示する。このディスプレイ108は、たとえば、CRT、TFT液晶ディスプレイ、プラズマディスプレイ等を採用することができる。
【0031】
I/F109は、通信回線を通じてインターネットなどのネットワークに接続され、このネットワークを介して他の装置に接続される。そして、I/F109は、ネットワークと内部のインターフェースを司り、外部装置からのデータの入出力を制御する。I/F109には、たとえばモデムやLANアダプタなどを採用することができる。
【0032】
キーボード110は、文字、数字、各種指示などの入力のためのキーを備え、データの入力をおこなう。また、タッチパネル式の入力パッドやテンキーなどであってもよい。マウス111は、カーソルの移動や範囲選択、あるいはウィンドウの移動やサイズの変更などをおこなう。ポインティングデバイスとして同様に機能を備えるものであれば、トラックボールやジョイスティックなどであってもよい。
【0033】
スキャナ112は、画像を光学的に読み取り、数値解析装置内に画像データを取り込む。また、プリンタ113は、画像データや文書データを印刷する。プリンタ113には、たとえば、レーザプリンタやインクジェットプリンタを採用することができる。
【0034】
(数値解析装置の機能的構成)
この発明の実施の形態にかかる数値解析装置の機能的構成について説明する。図2はこの発明の実施の形態にかかる数値解析装置の機能構成を示すブロック図である。図2に示すように、数値解析装置200は、データ格納部201と、第1の取得部202と、第2の取得部203と、生成部204と、第1の算出部205と、第2の算出部206と、解析部207と、から構成されている。
【0035】
データ格納部201は、解析対象の数値解析に用いる解析物理量情報テーブルを有する。ここで、解析物理量情報テーブルについて説明する。図3は、この発明の実施の形態にかかる解析物理量情報テーブルを示す説明図である。
【0036】
図3において、解析物理量情報テーブルには、解析対象についての、解析物理量IDと、解析物理量名(汚染物質、熱、流体、構造など)と、既知解析物理量とが格納されている。また、既知解析物理量としては、既知物性値と、境界解析物理量と、初期解析物理量とが格納されている。
【0037】
解析物理量IDおよび解析物理量名は、ユーザの入力操作によって指定される情報であり、解析物理量IDや解析物理量名が指定されると、関連付けられている既知物性値、境界解析物理量および初期解析物理量を抽出することができる。
【0038】
既知物性値は、ユーザの入力操作によって解析物理量IDまたは解析物理量名が指定された解析対象が、どのような物質であるかを決定する情報である。既知物性値としては、たとえば、拡散係数、密度、比熱、熱伝導係数、粘性係数、水底面摩擦係数、ヤング率、ポアソン比などが挙げられる。
【0039】
境界解析物理量とは、解析対象の各要素からなる解析領域(三角形や四面体に分割したメッシュ)がどのような条件下で解析がおこなわれるかという解析条件をあらわす情報である。境界解析物理量としては、たとえば、物質濃度、温度、流速、圧力、流量、水深、変位などが挙げられる。
【0040】
境界解析物理量については、解析対象がたとえば汚染物質や熱の場合、汚染物質や熱の発生(供給)源の部分で、物質濃度=発生(供給)源値、温度=発生(供給)源値となる。また、解析対象が流体で境界が壁ならば、流速=0、流量=0、解析対象が構造で境界が物体を固定ならば、変位=0に設定することができる。定常問題(解析物理量が時間に依存せず時々刻々と変化しない問題)であれば、上述した既知物性値および境界解析物理量によって、解析することができる。
【0041】
非定常問題(解析物理量が時間に依存し時々刻々と変化する問題)の場合は、さらに初期条件となる初期解析物理量が必要になる。初期解析物理量は、解析開始時の解析物理量の値である。また、データ格納部201において、図示はしないが、初期設定情報として、時間方向に対して、逐次解析計算をおこなうために時間増分量、計算ステップ数などの情報も、解析物理量情報テーブルに格納されている。
【0042】
さらに、解析対象範囲に分割形成されるメッシュの分割幅、節点数、要素数、節点番号に対応した節点座標値および要素番号に対応した要素結合情報などのメッシュ情報も格納されている。このメッシュは、解析対象範囲が2次元空間である場合は、三角形の形状、解析対象範囲が3次元空間である場合は、四面体の有限な形状の要素となる。メッシュは、図1に示したスキャナ112により読み取ることによっても取り込むことができる。このデータ格納部201は、具体的には、例えば、図1に示したRAM103、HD105またはFD107によりその機能を実現する。
【0043】
第1の取得部202は、解析対象の既知解析物理量を取得する。具体的には、入力装置(たとえば図1に示したキーボード110またはマウス111)を操作することにより、解析対象となる物理量の解析物理量IDを選択し、データ格納部201から既知解析物理量を抽出する。この第1の取得部202は、具体的には、たとえば図1に示したROM102、RAM103、HD105、FD107などに格納されたプログラムをCPU101が実行することによって、また図1に示したI/F109によってその機能を実現する。
【0044】
第2の取得部203は、解析対象の各要素(メッシュ)の整合質量行列を取得する。第2の取得部203は、具体的には、たとえば、図1に示したROM102、RAM103、HD105、FD107などに格納されたプログラムをCPU101が実行することによって、また図1に示したI/F109によってその機能を実現する。
【0045】
生成部204は、解析対象の各要素の気泡関数に基づいて、第2の取得部203によって取得された各要素の整合質量行列を対角化した、各要素の対角質量行列を生成する。生成部204は、具体的には、たとえば、図1に示したROM102、RAM103、HD105、FD107などに格納されたプログラムをCPU101が実行することによって、また図1に示したI/F109によってその機能を実現する。
【0046】
第1の算出部205は、生成部204によって生成された各要素の対角質量行列に基づいて、解析対象全域の対角質量行列を算出する。第1の算出部205は、具体的には、たとえば、図1に示したROM102、RAM103、HD105、FD107などに格納されたプログラムをCPU101が実行することによって、また図1に示したI/F109によってその機能を実現する。
【0047】
第2の算出部206は、第1の算出部205によって算出された、解析対象全域の対角質量行列の逆行列を算出する。第2の算出部206は、具体的には、たとえば、図1に示したROM102、RAM103、HD105、FD107などに格納されたプログラムをCPU101が実行することによって、また図1に示したI/F109によってその機能を実現する。
【0048】
解析部207は、第1の取得部202によって取得された既知解析物理量と、生成部204によって生成された各要素の対角質量行列と、に基づいて、解析対象の挙動を解析する。具体的には、第1の取得部202によって取得された既知解析物理量と、第1の算出部205によって算出された解析対象全域の対角質量行列と、第2の算出部206によって算出された逆行列と、に基づいて、解析対象の挙動を解析する。なお、解析対象全域の対角質量行列を必要とせず、各要素の対角質量行列のみを用いて解析対象の挙動を解析する解法があるが(たとえば、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参照)、この場合には、生成部204と第1の算出部205を同一視して、第2の算出部206へ進むことも可能である。解析部207は、具体的には、たとえば、図1に示したROM102、RAM103、HD105、FD107などに格納されたプログラムをCPU101が実行することによって、また図1に示したI/F109によってその機能を実現する。
【0049】
本実施の形態にかかる数値解析装置200の数値解析処理手順について説明する。図4はこの発明の実施の形態にかかる数値解析装置200の数値解析処理手順を示すフローチャートである。まず、第1の取得部202により、解析対象の既知解析物理量を取得する(ステップS401)。つぎに、第2の取得部203により、各要素(要素レベル)の整合質量行列を取得する(ステップS402)。
【0050】
続いて、要素ごとに、気泡関数を積分し(ステップS403)、各要素(要素レベル)の整合質量行列に、ステップS403で積分された値を代入することによって、各要素(要素レベル)の対角質量行列を算出する(ステップS404)。つぎに、各要素(要素レベル)の対角質量行列の総和(重ね合せ)により、解析対象全域の対角質量行列を算出する(ステップS405)。
【0051】
そして、この対角質量行列の逆行列を算出し(ステップS406)、解析対象の既知解析物理量と、解析対象全域の対角質量行列と、その逆行列とに基づいて、解析対象の挙動を解析する。具体的には、たとえば、上記式(10)~式(13)に代入することにより、解析することができる。
【0052】
数値解析装置200の具体的な解析処理内容について説明する。
<1.直交基底気泡関数要素>
質量行列Mは、有限要素法による解析(有限要素解析)について、ある解析物理量を補間した関数に重み関数を掛け、対象領域で積分を行った場合に現れる。気泡関数要素を使用した場合には、解析対象となる任意領域を区分的に三角形(四面体)形状にて要素数Neで分割した要素群を全体系に重ね合せる作業によって形成され、要素領域での積分を用いて下記式(16)のようにあらわすことができる。
【0053】
【数7】
JP0004729767B2_000008t.gif

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

【0056】
【数9】
JP0004729767B2_000010t.gif

【0057】
【数10】
JP0004729767B2_000011t.gif

【0058】
【数11】
JP0004729767B2_000012t.gif

【0059】
式(18)、(20)より以下の下記関係式(21)を得る。
【0060】
【数12】
JP0004729767B2_000013t.gif

【0061】
<2.直交条件を満たす気泡関数>
式(18)、式(20)を満たすことができる気泡関数として下記式(22)で表される気泡関数を提案する。
【0062】
【数13】
JP0004729767B2_000014t.gif

【0063】
式(18)に式(22)を代入することにより下記式(23)を得る。
【0064】
【数14】
JP0004729767B2_000015t.gif

【0065】
式(20)に式(22)を代入することにより下記式(24)を得る。
【0066】
【数15】
JP0004729767B2_000016t.gif

【0067】
式(23)に式(24)を代入することにより、最終的にα1を未知量とした下記式(25)を得る。
【0068】
【数16】
JP0004729767B2_000017t.gif

【0069】
式(25)はα1に関する二次方程式なので、解の公式より下記式(26)にて未知量α1を求めることができる。
【0070】
【数17】
JP0004729767B2_000018t.gif

【0071】
<3.直交基底を形成する具体的な気泡関数>
直交基底を形成する具体的な気泡関数について説明する。図5は、2次元の気泡関数要素を示す説明図であり、図6は、3次元の気泡関数要素を示す説明図である。
【0072】
直交基底を形成する具体的な気泡関数を導出するため、図5および図6に示した三角形(四面体)の要素領域を、重心点を用いて3(4)つの小三角形(小四面体)wi,i=1・・・N+1に分割するξ乗気泡関数(0<ξ<∞)を用いる。ξ乗気泡関数は小三角形(小四面体)ごとにアイソパラメトリック座標系{r,s}({r,s,t})を用いて次の式(27)および式(28)のように定義される。
【0073】
【数18】
JP0004729767B2_000019t.gif

【0074】
これらの気泡関数の要素領域における積分値は、以下の式(29)によって求めることができる。
【0075】
【数19】
JP0004729767B2_000020t.gif

【0076】
式(22)を、ξ乗気泡関数を用いて下記式(30)のように定義する。
【0077】
【数20】
JP0004729767B2_000021t.gif

【0078】
式(30)の各べき乗の値を下記式(31)のように設定すると二次元(三角形)、三次元(四面体)の気泡関数に対してα1,α2は下記式(32)~式(35)のようになる。
【0079】
【数21】
JP0004729767B2_000022t.gif

【0080】
【数22】
JP0004729767B2_000023t.gif

【0081】
式(30)の各べき乗の値を下記式(36)のように設定する(式(31)とは異なった値)と二次元(三角形)、三次元(四面体)の気泡関数に対してα1, α2は下記式(37)~式(40)のようになる。
【0082】
【数23】
JP0004729767B2_000024t.gif

【0083】
【数24】
JP0004729767B2_000025t.gif

【0084】
図7は、式(32)のα12を使用した直交基底を形成する三角形気泡関数の形状φBを示す説明図であり、図8は、式(32)のα12を使用した直交基底を形成する三角形気泡関数の形状φB2を示す説明図である。また、図9は、式(37)のα12を使用した直交基底を形成する三角形気泡関数の形状φBを示す説明図であり、図10は、式(37)のα12を使用した直交基底を形成する三角形気泡関数の形状φB2を示す説明図である。
【0085】
式(32)、式(37)および図7~図10に示すように、式(21)を満たす気泡関数(形状)は無数に存在するため、直交基底を形成する具体的な気泡関数やその形状には殆ど意味がなく、直交基底気泡関数要素の本質は気泡関数要素の基底が直交する条件式(21)を発明した部分にある。
【0086】
<4.直交基底気泡関数要素の要素レベルの質量行列>
下記式(41)~式(46)に、2次元(三角形)、3次元(四面体)における通常の気泡関数要素の要素レベルの整合質量行列、集中質量行列、直交基底気泡関数要素の要素レベルの質量行列(対角質量行列)を示す。
【0087】
【数25】
JP0004729767B2_000026t.gif

【0088】
式(43)の要素レベルの対角質量行列M(e)ijは、生成部204により、式(41)の要素レベルの整合質量行列M(e)ijの<1,φeおよび∥φB2Ωeに、気泡関数φBの積分値CAeを代入することによって生成される。そして、第1の算出部205により、各要素の要素レベルの対角質量行列M(e)ijの総和(重ね合せ)をとることによって、解析対象全域の対角質量行列Mを算出することができる。また、第2の算出部206により、対角質量行列Mの逆行列M-1を算出することができる。そして、解析部207により、既知解析物理量u、対角質量行列Mおよび逆行列M-1を、上述した式(10)~式(13)に代入することによって、解析対象の挙動を解析することができる。
【0089】
【数26】
JP0004729767B2_000027t.gif

【0090】
式(46)の要素レベルの対角質量行列は、生成部204により、式(44)の要素レベルの整合質量行列の<1,φeおよび∥φB2Ωeに、気泡関数φBの積分値CAeを代入することによって生成される。そして、第1の算出部205により、各要素の要素レベルの対角質量行列M(e)ijの総和(重ね合せ)をとることによって、解析対象全域の対角質量行列Mを算出することができる。また、第2の算出部206により、対角質量行列Mの逆行列M-1を算出することができる。そして、解析部207により、既知解析物理量u、対角質量行列Mおよび逆行列M-1を、上述した式(10)~式(13)に代入することによって、解析対象の挙動を解析することができる。
【0091】
<5.気泡関数について>
気泡関数について説明する。図11は、2次元の気泡関数を示す説明図であり、図12は、3次元の気泡関数を示す説明図である。気泡関数φBは、形状関数の定義(重心点が1、その他の点では0となる関数)を満たせば任意の関数を選択することができる。検証の比較計算として、最も多く使用されている以下に示す二つの気泡関数および直交基底を形成する気泡関数を用いる。図11および図12に示したアイソパラメトリック座標系{r,s}({r,s,t})を用いて、式(47)~式(52)のように定義される。
【0092】
【数27】
JP0004729767B2_000028t.gif

【0093】
【数28】
JP0004729767B2_000029t.gif

【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段解法を、微小時間Δ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段解法を、微小時間Δ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次元まで空間次元数Nを用いて統一的に記述する。また、一次元(線気泡関数要素)の要素レベルの質量行列を示す。
【0110】
図41は、アイソパラメトリック座標系rを示す説明図(その1)である。図41に示したアイソパラメトリック座標系rを、式(1)、(2)の表現で記述すると、1次元のΦT, uT,Ψαは下記式(53)~(55)になる。
【0111】
【数29】
JP0004729767B2_000030t.gif

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

【0114】
式(56)、(27)、(28)のξ(η)乗気泡関数(0<η<∞)は、下記式(57)に示すように、気泡関数を空間で偏微分した積分値も統一的に記述でき、下記式(58)、(59)のような積分値の関係を持っている。
【0115】
【数31】
JP0004729767B2_000032t.gif

【0116】
式(57)のxk,xl=x1・・・xNはそれぞれ空間方向を表し、式(59)のαkdは実数を、KDは正の整数を示している。式(59)は式(58)より誘導している。C=(N+1)/(N+2)とし、空間次元数Nを使用して直交基底気泡関数要素の条件式(14)(または式(21))を1次元~3次元まで統一的に記述すると下記式(60)になる。
【0117】
【数32】
JP0004729767B2_000033t.gif

【0118】
下記式(61)~(63)に、1次元(線)における通常の気泡関数要素の要素レベルの整合質量行列、集中質量行列、直交基底気泡関数要素の要素レベルの質量行列(対角質量行列)を示す。
【0119】
【数33】
JP0004729767B2_000034t.gif

【0120】
<9.拡張された直交基底気泡関数要素>
気泡関数要素を使用した流体解析、構造解析、固有値解析などでは、一般に一次要素の形状関数Ψαの積分値に加えて、下記式(64)~(66)の気泡関数の積分値が必要になる。
【0121】
【数34】
JP0004729767B2_000035t.gif

【0122】
下記式(67)と式(64)より、下記の関係式(68)を得る。
【0123】
【数35】
JP0004729767B2_000036t.gif

【0124】
気泡関数として、下記式(69)の関係を持っている気泡関数の場合には、必要となる気泡関数の積分値(64)~(66)の式(64)は、式(70)に代用することができる。
【0125】
【数36】
JP0004729767B2_000037t.gif

【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】
JP0004729767B2_000038t.gif

【0129】
式(71)をみると、D(ξ123)の値が、変数ξ1、ξ2、ξ3の関数となっている。実際に、式(31)、(36)の値を用いて、式(72)を計算(α1(ξ123)=α1+(ξ123)を使用)してみると、下記式(73)~(75)のように採用する変数によって計算される値が異なる。
【0130】
【数38】
JP0004729767B2_000039t.gif

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

【0133】
基底が直交する条件式(60)に、式(76)の条件を追加した下記式(80)を拡張された直交基底気泡関数要素の条件式とする。
【0134】
【数40】
JP0004729767B2_000041t.gif

【0135】
<10.拡張された直交基底気泡関数要素を満たす気泡関数>
式(71)からわかるように、式(56)、(27)、(28)のξ乗気泡関数を使用して直交基底気泡関数を形成した式(66)の積分値は、変数ξ1、ξ2、ξ3の値によって決まる。本特許では、拡張された直交基底気泡関数要素の条件式(80)を満足するような変数ξ1、ξ2、ξ3を考え、そのような変数が実際に存在することを示す。変数はξ1、ξ2、ξ3の三つあるが、その内の2変数を下記式(81)のように固定し、ξ1の1変数の問題として、拡張された直交基底気泡関数要素の条件を満足する気泡関数を求めることを考える。
【0136】
【数41】
JP0004729767B2_000042t.gif

【0137】
ξ1の1変数の問題として、下記式(82)を定義する。
【0138】
【数42】
JP0004729767B2_000043t.gif

【0139】
拡張された直交基底気泡関数要素の条件を満たす気泡関数が存在することは、実数の集合Rにおいて、f(ξ1)=0となる変数が存在することである。いま、下記式(83)~(85)のような実数の集合Rの部分集合Sを考える。
【0140】
【数43】
JP0004729767B2_000044t.gif

【0141】
ξ1∈Sでは、1次元~3次元において、下記式(86)となる。
【0142】
【数44】
JP0004729767B2_000045t.gif

【0143】
従って、α1、α2は実数の範囲で値を持つ。また、ξ1∈Sにて1次元~3次元では、それぞれ下記式(87)および式(88)が成り立つので、式(72)は有界な値をとり、式(82)は実変数の関数となる。
【0144】
【数45】
JP0004729767B2_000046t.gif

【0145】
式(82)と部分集合Sを使用すると下記式(89)~(91)を得ることが出来る。
【0146】
【数46】
JP0004729767B2_000047t.gif

【0147】
任意のξ1c∈Sに対して、1次元では、式(89)の第1式が負、第2式が正で、第3式が正(単調増加)、2次元、3次元では、式(90)、(91)の第1式が正、第2式が負で、第3式が負(単調減少)であるので、1次元~3次元において、部分集合Sの範囲には、それぞれf(ξ1)=0の解が一意に存在する。
【0148】
具体的に、部分集合Sの範囲内で、二分法を使用してf(ξ1)=0の解を求めると、下記式(92)~(94)の値が得られる。
【0149】
【数47】
JP0004729767B2_000048t.gif

【0150】
図43、図44に、式(81)、(93)を用いた三角形気泡関数の形状φB、三角形気泡関数の形状φB2を示す。
【0151】
または、下記式(95)、(96)を使用して、f(ξ1)=0の解の存在を示すことができる。
【0152】
【数48】
JP0004729767B2_000049t.gif

【0153】
1次元~3次元において、hを正の定数とし、任意のξ1c∈Sに対して下記式(97)の結果を得る。
【0154】
【数49】
JP0004729767B2_000050t.gif

【0155】
平均値の定理から、任意のξ1a, ξ1b∈Sに対して、下記式(98)が得られる。
【0156】
【数50】
JP0004729767B2_000051t.gif

【0157】
従って、縮小写像の原理(定理)より、式(95)を反復的に逐次計算をさせる下記式(99)にて得られる数列{ξ1m}は、部分集合Sの範囲内でf(ξ1)=0を満足する一意の解に収束する。
【0158】
【数51】
JP0004729767B2_000052t.gif

【0159】
なお、拡張された直交基底気泡関数要素の存在を示すために、式(80)の三つの条件を満たす気泡関数として、三つの未知量をα1、α2、ξ1とし、適当な気泡関数φB1、φB2、φB3を採用しが、他にも、下記式(100)の気泡関数φBのように、未知量をα1、α2、α3として、適当な気泡関数φB1、φB2、φB3、φB4を採用すれば、式(80)の条件を満たす気泡関数が得られる。このように、拡張された直交基底気泡関数要素の条件式を満たす気泡関数は、唯一ではなく、幾つも存在するため、具体的な気泡関数やその形状には殆ど意味がなく、拡張された直交基底気泡関数要素で重要な意味を成すのは、式(80)の条件となる気泡関数の積分値である。
【0160】
【数52】
JP0004729767B2_000053t.gif

【0161】
<11.拡張された直交基底気泡関数要素の要素レベルの粘性項(拡散項)の行列>
通常の気泡関数、直交基底気泡関数要素の数値解析において式(66)が必要になる要素レベルの粘性項(拡散項)の行列を下記式(101)~(103)に示す。
【0162】
【数53】
JP0004729767B2_000054t.gif

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

【0165】
以上、拡張された直交基底気泡関数要素の要素レベルの粘性項(拡散項)の行列(104)~(106)を示したが、本発明に関する定式化やプログラムでは、行列(101)~(103)に式(76)を代入して整理を行わずに、そのまま行列(101)~(103)と式(76)を使用することも、もちろん可能である。
【0166】
<12.直交基底気泡関数要素と正規化気泡関数>
式(1)は、下記式(107)、(108)のように変形することができる。
【0167】
【数55】
JP0004729767B2_000056t.gif

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

【0170】
さらに、拡張された直交基底気泡関数要素では、下記式(113)が得られる。
【0171】
【数57】
JP0004729767B2_000058t.gif

【0172】
<13.直交基底気泡関数要素と安定化作用制御項>
気泡関数要素の持つ数値不安定性を緩和させるために、下記式(114)の気泡関数の粘性項(安定化作用制御項)を用いて解析精度を向上させる方法がある。
【0173】
【数58】
JP0004729767B2_000059t.gif

【0174】
式(114)の安定化作用制御パラメータν’(-∞<ν’<∞)は、粘性係数の次元(粘性係数に相当する変数を無次元化している場合は無次元)を持つパラメータであり、この値を適切に決定することによって高精度解法を実現する。この方法は、従来、気泡関数要素に通常のGalerkin法を使用して、得られた有限要素方程式に重心点のみの粘性項(安定化作用制御項)を付加していたが(松本純一,梅津剛,「気泡関数を用いた非圧縮性流体の有限要素法解析に関する一考察」,日本応用数理学会平成8年度年会,1996年,46頁-47頁、松本純一,梅津剛,川原睦人,「線形型気泡関数を用いた非圧縮性粘性流体解析と適応型有限要素法」,応用力学論文集(土木学会),2巻,1999年8月,223頁-232頁参照)、現在では、重み関数に関してもう一つの気泡関数(3レベル気泡関数)を採用して重心点のみの粘性項(安定化作用制御項)を導入している(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、非特許文献3参照)。3レベル気泡関数は、下記式(115)~(117)を満たす気泡関数である。
【0175】
【数59】
JP0004729767B2_000060t.gif

【0176】
式(67)、(115)より下記式(118)を得る。
【0177】
【数60】
JP0004729767B2_000061t.gif

【0178】
また、下記式(119)の関係を持っている3レベル気泡関数の場合には、条件式(115)~(117)の式(115)は式(120)で代用できる。
【0179】
【数61】
JP0004729767B2_000062t.gif

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

【0182】
式中のνは、粘性係数(拡散係数)である。式(121)を式(120)、(116)、(117)に代入して整理すると、下記式(122)を得る。
【0183】
【数63】
JP0004729767B2_000064t.gif

【0184】
式(121)右辺の四つの気泡関数を具体的に下記式(123)のように定義する。
【0185】
【数64】
JP0004729767B2_000065t.gif

【0186】
式(123)、(59)より式(121)は式(119)を満足していることがわかる。式(123)の変数を下記式(124)のように設定する。
【0187】
【数65】
JP0004729767B2_000066t.gif

【0188】
α1(ξ123)=α1+(ξ123)として、式(31)の変数を用いると、式(122)より下記式(125)~(127)を得る。
【0189】
【数66】
JP0004729767B2_000067t.gif

【0190】
α1(ξ123)=α1+(ξ123)として、式(36)の変数を用いると、式(122)より下記式(128)~(130)を得る。
【0191】
【数67】
JP0004729767B2_000068t.gif

【0192】
α1(ξ123)=α1+(ξ123)として、式(92)~(94)、式(81)の変数を用いると、式(122)より下記式(131)~(133)を得る。
【0193】
【数68】
JP0004729767B2_000069t.gif

【0194】
図45、図46に、式(132)を用いた三角形3レベル気泡関数の形状(ν’=νとした場合)、式(93)、(81)を用いた三角形気泡関数と式(132)を用いた三角形3レベル気泡関数の積の形状(ν’=νとした場合)を示す。以上により、直交基底気泡関数要素を用いた場合にも、式(120)(あるいは式(115))、(116)、(117)を満たす3レベル気泡関数は存在し、安定化作用制御項(114)を導入できる。
【0195】
なお、3レベル気泡関数の存在を示すため、式(120)、(116)、(117)の三つの条件を満足する三つの未知量δ1、δ2、δ3と四つの適当な気泡関数を持った3レベル気泡関数(121)を採用したが、この他にも、下記式(134)のように、安定化作用制御パラメータν’と同じ次元を持つパラメータを導入し、未知量をδ1、δ2として、式(134)の右辺に適当な三つの気泡関数を用いることによって、式(120)、(116)の条件を満たす3レベル気泡関数を求めることができる。そして、下記式(135)のように安定化作用制御パラメータν’を定義すれば、式(117)も導出でき、式(120)(あるいは式(115))、(116)、(117)となる3レベル気泡関数の存在を示すことが出来る。このように、3レベル気泡関数の条件式を満たす気泡関数は、唯一ではなく、幾つも存在するため、具体的な気泡関数やその形状には殆ど意味がなく、3レベル気泡関数で重要な意味を成すのは、安定化作用制御項を導くために必要な式(120)(あるいは式(115))、(116)、(117)の条件である。
【0196】
【数69】
JP0004729767B2_000070t.gif

【0197】
<14.直交基底気泡関数要素と正規化気泡関数による安定化作用制御項>
正規化気泡関数を使用した場合にも安定化作用制御項は導入することが出来、下記式(136)によって表される。
【0198】
【数70】
JP0004729767B2_000071t.gif

【0199】
直交基底気泡関数要素で使用する気泡関数として正規化気泡関数を用いた場合、3レベル正規化気泡関数は下記式(137)と定めることができ、その積分値は下記式(138)~(142)のように記述できる。
【0200】
【数71】
JP0004729767B2_000072t.gif

【0201】
上記式(138)~(142)より、式(115)~(117)、(119)、(120)を満足する3レベル気泡関数が存在すれば、式(138)~(142)も自動的に満たすので3レベル正規化気泡関数も同時に存在することがわかる。
【0202】
<15.直交基底気泡関数要素と気泡関数の安定化作用行列>
式(107)の表現形式を用いた気泡関数要素の重心点の自由度uB’は、静的縮約という操作によって要素ごとに消去できる。その結果得られる有限要素方程式は、安定化有限要素法で導かれる有限要素方程式と密接な関係がある(非特許文献2参照)。静的縮約によって得られる気泡関数の安定化作用行列は、一般的に下記式(143)の形式によって表される。
【0203】
【数72】
JP0004729767B2_000073t.gif

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

【0206】
式(146)は解析精度が向上することが期待できる安定化行列である。式(145)の関係を満たすような安定化作用制御パラメータν’を求めるために、下記式(148)の評価関数を最小にするような安定化作用制御パラメータν’を決定する問題を考える。
【0207】
【数74】
JP0004729767B2_000075t.gif

【0208】
式(148)を最小にするような安定化作用制御パラメータν’を求めるために、式(148)の停留条件である下記関係式(150)を用いて、安定化作用制御パラメータν’を求める。
【0209】
【数75】
JP0004729767B2_000076t.gif

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

【0212】
式(151)の係数δ1、δ2、δ3は、式(122)によって決定することが出来、式(151)は式(120)、(116)を満たす。また、式(117)の代わりに下記式(152)の安定化作用制御項を得る。
【0213】
【数77】
JP0004729767B2_000078t.gif

【0214】
式(152)の安定化作用制御項を用いた場合、式(144)は下記式(153)のようになる。
【0215】
【数78】
JP0004729767B2_000079t.gif

【0216】
式(145)の関係を満たすような安定化作用制御パラメータνm’を求めるために、下記式(154)の評価関数を最小にするような安定化作用制御パラメータνm’を決定する問題を考える。
【0217】
【数79】
JP0004729767B2_000080t.gif

【0218】
式(154)を最小にするような安定化作用制御パラメータνm’を求めるために、式(154)の停留条件である下記関係式(156)を用いて、安定化作用制御パラメータνm’を求める。
【0219】
【数80】
JP0004729767B2_000081t.gif

【0220】
式(155)はベクトルであり、ベクトルの成分を下記式(157)のように表現する。
【0221】
【数81】
JP0004729767B2_000082t.gif

【0222】
式(157)をみるとわかるように、安定化作用制御パラメータνm’は各ベクトル成分で独立に定義されているので、式(145)の関係を満たすような安定化作用制御パラメータνm’を決定する問題として、式(157)の成分mごとに、下記式(158)の評価関数を用いた最小化問題を考えることができる。
【0223】
【数82】
JP0004729767B2_000083t.gif

【0224】
式(158)を最小にするような安定化作用制御パラメータνm’を求めるために、式(158)の停留条件である下記関係式(159)を用いて、安定化作用制御パラメータνm’を求める。
【0225】
【数83】
JP0004729767B2_000084t.gif

【0226】
式(157)は、安定化作用制御パラメータνm’が各ベクトル成分で独立に定義されているので、式(154)、(158)の停留条件から導かれる安定化作用制御パラメータνm’を求める式(156)、(159)は、結果的に同じ式になる。
【0227】
他の方法として、下記式(160)の関係を満たすように、安定化作用制御パラメータν’を決定する問題も考えられる。
【0228】
【数84】
JP0004729767B2_000085t.gif

【0229】
式(160)の関係を満たすような安定化作用制御パラメータν’を求めるために、下記式(161)の評価関数を最小にするような安定化作用制御パラメータν’を決定する問題を考える。
【0230】
【数85】
JP0004729767B2_000086t.gif

【0231】
式(161)を最小にするような安定化作用制御パラメータν’を求めるために、式(161)の停留条件である下記関係式(163)を用いて、安定化作用制御パラメータν’を求める。
【0232】
【数86】
JP0004729767B2_000087t.gif

【0233】
安定化作用制御パラメータν’は、式(151)の3レベル気泡関数を用いると、式(152)の安定化作用制御項が導出され、求めるべき未知解析物理量の数Lだけ、安定化作用制御パラメータνm’を定義することができる。式(160)の関係を満たすような安定化作用制御パラメータνm’を求めるために、下記式(164)の評価関数を最小にするような安定化作用制御パラメータνm’を決定する問題を考える。
【0234】
【数87】
JP0004729767B2_000088t.gif

【0235】
式(164)を最小にするような安定化作用制御パラメータνm’を求めるために、式(164)の停留条件である下記関係式(166)を用いて、安定化作用制御パラメータνm’を求める。
【0236】
【数88】
JP0004729767B2_000089t.gif

【0237】
式(165)はベクトルであり、ベクトルの成分を下記式(167)のように表現する。
【0238】
【数89】
JP0004729767B2_000090t.gif

【0239】
式(167)をみるとわかるように、安定化作用制御パラメータνm’は各ベクトル成分で独立に定義されているので、式(160)の関係を満たすような安定化作用制御パラメータνm’を決定する問題として、式(167)の成分mごとに、下記式(168)の評価関数を用いた最小化問題を考えることができる。
【0240】
【数90】
JP0004729767B2_000091t.gif

【0241】
式(168)を最小にするような安定化作用制御パラメータνm’を求めるために、式(168)の停留条件である下記関係式(169)を用いて、安定化作用制御パラメータνm’を求める。
【0242】
【数91】
JP0004729767B2_000092t.gif

【0243】
式(167)は、安定化作用制御パラメータνm’が各ベクトル成分で独立に定義されているので、式(164)、(168)の停留条件から導かれる安定化作用制御パラメータνm’を求める式(166)、(169)は、結果的に同じ式になる。
【0244】
<16.直交基底気泡関数要素と正規化気泡関数の安定化作用行列>
式(108)の表現形式を用いた気泡関数要素の重心点の自由度uS’は、静的縮約という操作によって要素ごとに消去できる。静的縮約によって得られる正規化気泡関数の安定化作用行列は、一般的に下記式(170)の形式によって表される。
【0245】
【数92】
JP0004729767B2_000093t.gif

【0246】
式(170)の第二項は、式(141)の関係を満たす気泡関数を使用した場合、式(111)~(113)各第一項の正規化気泡関数の積分値から成る行列であり、式(70)、(65)、(66)で表した場合には第三項のように記述できる。安定化作用制御パラメータν’は、下記式(172)の関係が満たされるように決定する。
【0247】
【数93】
JP0004729767B2_000094t.gif

【0248】
式(172)の関係を満たすような安定化作用制御パラメータν’を求めるために、下記式(174)の評価関数を最小にするような安定化作用制御パラメータν’を決定する問題を考える。
【0249】
【数94】
JP0004729767B2_000095t.gif

【0250】
式(174)を最小にするような安定化作用制御パラメータν’を求めるために、式(174)の停留条件である下記関係式(176)を用いて、安定化作用制御パラメータν’を求める。
【0251】
【数95】
JP0004729767B2_000096t.gif

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

【0254】
式(177)の係数δ1、δ2、δ3は、式(122)によって決定することが出来、式(177)は式(142)、(139)を満たす。また、式(140)の代わりに下記式(178)の安定化作用制御項を得る。
【0255】
【数97】
JP0004729767B2_000098t.gif

【0256】
式(178)の安定化作用制御項を用いた場合、式(171)は下記式(179)のようになる。
【0257】
【数98】
JP0004729767B2_000099t.gif

【0258】
式(179)の第一項は、式(111)~(113)各第一項の正規化気泡関数の積分値から成る行列であり、式(70)、(65)(66)で表した場合には第二項のように記述できる。式(172)の関係を満たすような安定化作用制御パラメータνm’を求めるために、下記式(180)の評価関数を最小にするような安定化作用制御パラメータνm’を決定する問題を考える。
【0259】
【数99】
JP0004729767B2_000100t.gif

【0260】
式(180)を最小にするような安定化作用制御パラメータνm’を求めるために、式(180)の停留条件である下記関係式(182)を用いて、安定化作用制御パラメータνm’を求める。
【0261】
【数100】
JP0004729767B2_000101t.gif

【0262】
式(181)はベクトルであり、ベクトルの成分を下記式(183)のように表現する。
【0263】
【数101】
JP0004729767B2_000102t.gif

【0264】
式(183)をみるとわかるように、安定化作用制御パラメータνm’は各ベクトル成分で独立に定義されているので、式(172)の関係を満たすような安定化作用制御パラメータνm’を決定する問題として、式(183)の成分mごとに、下記式(184)の評価関数を用いて最小化問題を考えることができる。
【0265】
【数102】
JP0004729767B2_000103t.gif

【0266】
式(184)を最小にするような安定化作用制御パラメータνm’を求めるために、式(184)の停留条件である下記関係式(185)を用いて、安定化作用制御パラメータνm’を求める。
【0267】
【数103】
JP0004729767B2_000104t.gif

【0268】
式(183)は、安定化作用制御パラメータνm’が各ベクトル成分で独立に定義されているので、式(180)、(184)の停留条件から導かれる安定化作用制御パラメータνm’を求める式(182)、(185)は、結果的に同じ式になる。
【0269】
他の方法として、下記式(186)の関係を満たすように、安定化作用制御パラメータν’を決定する問題も考えられる。
【0270】
【数104】
JP0004729767B2_000105t.gif

【0271】
式(186)の関係を満たすような安定化作用制御パラメータν’を求めるために、下記式(187)の評価関数を最小にするような安定化作用制御パラメータν’を決定する問題を考える。
【0272】
【数105】
JP0004729767B2_000106t.gif

【0273】
式(187)を最小にするような安定化作用制御パラメータν’を求めるために、式(187)の停留条件である下記関係式(189)を用いて、安定化作用制御パラメータν’を求める。
【0274】
【数106】
JP0004729767B2_000107t.gif

【0275】
安定化作用制御パラメータν’は、式(177)の3レベル正規化気泡関数を用いると、式(178)の安定化作用制御項が導出され、求めるべき未知解析物理量の数Lだけ、安定化作用制御パラメータνm’を定義することができる。式(186)の関係を満たすような安定化作用制御パラメータνm’を求めるために、下記式(190)の評価関数を最小にするような安定化作用制御パラメータνm’を決定する問題を考える。
【0276】
【数107】
JP0004729767B2_000108t.gif

【0277】
式(190)を最小にするような安定化作用制御パラメータνm’を求めるために、式(190)の停留条件である下記関係式(192)を用いて、安定化作用制御パラメータνm’を求める。
【0278】
【数108】
JP0004729767B2_000109t.gif

【0279】
式(191)はベクトルであり、ベクトルの成分を下記式(193)のように表現する。
【0280】
【数109】
JP0004729767B2_000110t.gif

【0281】
式(193)をみるとわかるように、安定化作用制御パラメータνm’は各ベクトル成分で独立に定義されているので、式(186)の関係を満たすような安定化作用制御パラメータνm’を決定する問題として、式(193)の成分mごとに、下記式(194)の評価関数を用いた最小化問題を考えることができる。
【0282】
【数110】
JP0004729767B2_000111t.gif

【0283】
式(194)を最小にするような安定化作用制御パラメータνm’を求めるために、式(194)の停留条件である下記関係式(195)を用いて、安定化作用制御パラメータνm’を求める。
【0284】
【数111】
JP0004729767B2_000112t.gif

【0285】
式(193)は、安定化作用制御パラメータνm’が各ベクトル成分で独立に定義されているので、式(190)、(194)の停留条件から導かれる安定化作用制御パラメータνm’を求める式(192)、(195)は、結果的に同じ式になる。
【0286】
<17.気泡関数と正規化気泡関数の安定化作用制御パラメータの関係>
直交基底気泡関数要素で、式(1)、(107)の表現形式を採用した場合には、式(70)、(65)、(66)の積分を使用して、式(150)、(156)、(159)、(163)、(166)、(169)のいずれかを採用し、式(108)の表現形式を採用した場合には、式(111)~(113)各第一項の積分を使用して、式(176)、(182)、(185)、(189)、(192)、(195)のいずれかを採用すれば安定化作用制御パラメータを計算することが出来る。いま、式(176)、(182)、(185)、(189)、(192)、(195)を、式(70)、(65)、(66)の積分を使用して記述すると下記式(196)~(201)のようになる。
【0287】
【数112】
JP0004729767B2_000113t.gif

【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)気泡関数、気泡関数自体の積分値がAeとなる正規化気泡関数、2)気泡関数の持つ数値不安定を緩和させるために導入する安定化作用制御項、および3レベル気泡関数、3レベル正規化気泡関数、3)安定化作用制御パラメータの計算で用いる安定化作用行列、の定式化やプログラミングにおいて、式(70)、(65)、(66)の気泡関数の積分値、あるいは式(111)~(113)各第一項の正規化気泡関数の積分値が必要になる。拡張された直交基底気泡関数要素は、必要となる積分値が全て事前に(複雑な積分を行うことなく)式(80)によって与えられているので、整合質量行列と同等の解析精度を持つ対角質量行列の作成や気泡関数の持つ数値不安定性を改良する安定化作用制御項の導入で現れる一連の作業において、非常に合理的に処理を行うことが出来る。
【0290】
直交基底気泡関数要素は、2次元、3次元の解析結果からわかるように、記憶容量、計算時間の効率の良い解法(例えば、4段解法)を採用しながら、信頼性のある結果を得ることができる。さらに、拡張された直交基底気泡関数要素の条件式(80)により、インクジェットの液滴挙動解析、都市や河川の氾濫解析、地震等で発生する津波の解析のみならず、構造物の応力解析、自動車のエンジン固有振動解析などで必要となる気泡関数要素の積分値や気泡関数の持つ数値不安定性を改良するために必要となる処理を合理的、かつ統一的に扱うことができるので、上記のような多岐にわたる数値解析への適用が容易に行えるという効果をもたらす。
【0291】
以上説明したように、本実施の形態によれば、気泡関数要素を使用した数値解析において、記憶容量、計算時間などの計算効率の良い解析手法(直交基底気泡関数要素数値解析方法)を実現することができる。
【0292】
なお、本実施の形態で説明した直交基底気泡関数要素数値解析方法は、予め用意されたプログラムをパーソナル・コンピュータ、ワークステーション、PCクラスタ、スーパーコンピュータ等のコンピュータで実行することにより実現できる。このプログラムは、ハードディスク、フレキシブルディスク、CD-ROM、MO、DVD等のコンピュータで読み取り可能な記録媒体に記録され、コンピュータにて記録媒体から読み出されることによって実行される。またこのプログラムは、インターネット等のネットワークを介して配布することが可能な伝送媒体であってもよい。
【産業上の利用可能性】
【0293】
以上のように、本発明にかかる直交基底気泡関数要素数値解析方法、直交基底気泡関数要素数値解析プログラムおよび直交基底気泡関数要素数値解析装置は、線、三角形および四面体形状のメッシュを使用した有限要素法に基づく数値解析に有用である。

図面
【図1】
0
【図2】
1
【図3】
2
【図4】
3
【図5】
4
【図6】
5
【図7】
6
【図8】
7
【図9】
8
【図10】
9
【図11】
10
【図12】
11
【図13】
12
【図14】
13
【図15】
14
【図16】
15
【図17】
16
【図18】
17
【図19】
18
【図20】
19
【図21】
20
【図22】
21
【図23】
22
【図24】
23
【図25】
24
【図26】
25
【図27】
26
【図28】
27
【図29】
28
【図30】
29
【図31】
30
【図32】
31
【図33】
32
【図34】
33
【図35】
34
【図36】
35
【図37】
36
【図38】
37
【図39】
38
【図40】
39
【図41】
40
【図42】
41
【図43】
42
【図44】
43
【図45】
44
【図46】
45
【図47】
46
【図48】
47
【図49】
48
【図50】
49
【図51】
50
【図52】
51