(57)【要約】【課題】  要素分割のデータ入力に際し、データ入力作業から各要素の節点の局所番号付けする作業をなくし、高次要素を用いた計算であっても、データを与えるとほとんど誤差なく短時間に、剛性行列Ｋと質量行列Ｍを直接計算して求め、固有関数を精度良く可視化して出力することができる有限要素法直接計算プログラムおよび解析方法を提供することを目的とする。【解決手段】  有限要素法により解析対象物の数値モデルを作成し、数値解析を行い、解析結果を可視化するプログラムの解析方法において、有限要素法の基底関数から決まる剛性行列Ｋと質量行列Ｍを直接計算することにより、与えられた領域のラプラシアンの固有値問題の固有値と固有関数を求める。【選択図】  図１
2006313400.tif 000002
【特許請求の範囲】
【請求項１】
有限要素法により解析対象物の数値モデルを作成し、数値解析を行い、解析結果を可視化するプログラムにおいて、
前記プログラムは有限要素法の基底関数から決まる剛性行列Ｋと質量行列Ｍを直接計算することにより、与えられた領域のラプラシアンの固有値問題の固有値と固有関数を求める
ことを特徴とするプログラム。
【請求項２】
前記解析対象物が１次元モデルである場合に、解析対象物を１次元モデルによる解析方法より解析し、
前記解析対象物が２次元モデルである場合に、解析対象物を１次元モデルによる解析方法、２次元モデルによる解析方法又はこれらの組み合わせにより解析し、
前記解析対象物が３次元モデルである場合に、前記解析対象物を、１次元モデルによる解析方法、２次元モデルによる解析方法、３次元モデルによる解析方法又はこれらの任意の組み合わせにより解析する
ことを特徴とする請求項１に記載のプログラム。
【請求項３】
有限要素法により解析対象物の数値モデルを作成し、数値解析を行い、解析結果を可視化する解析方法において、
前記解析方法は有限要素法の基底関数から決まる剛性行列Ｋと質量行列Ｍを直接計算することにより、与えられた領域のラプラシアンの固有値問題の固有値と固有関数を求める
ことを特徴とする解析方法。
【請求項４】
前記解析対象物が１次元モデルである場合に、解析対象物を１次元モデルによる解析方法より解析し、
前記解析対象物が２次元モデルである場合に、解析対象物を１次元モデルによる解析方法、２次元モデルによる解析方法又はこれらの組み合わせにより解析し、
前記解析対象物が３次元モデルである場合に、前記解析対象物を、１次元モデルによる解析方法、２次元モデルによる解析方法、３次元モデルによる解析方法又はこれらの任意の組み合わせにより解析する
ことを特徴とする請求項３に記載の解析方法。
【請求項５】
前記解析方法は、剛性行列Ｋと質量行列Ｍを直接計算する場合に、
前記解析対象物が２次元の平面のモデルであるとき、
剛性行列の要素Ｋｉｊについて、
2006313400.tif 000041
の定義により、
剛性行列の要素Ｋｉｊを決定する式が、
ｉ≠ｊのとき、
2006313400.tif 000042
ｉ＝ｊのとき、
2006313400.tif 000043
のように得られる
ことを特徴とする請求項３又は４に記載の解析方法。
【請求項６】
前記解析方法は、剛性行列Ｋと質量行列Ｍを直接計算する場合に、
前記解析対象物が２次元の平面のモデルであるとき、
質量行列の要素Ｍｉｊについて、
2006313400.tif 000044
の定義により、
質量行列の要素Ｍｉｊを決定する式が、
ｉ≠ｊのとき、
2006313400.tif 000045
ｉ＝ｊのとき、
2006313400.tif 000046
のように得られる
ことを特徴とする請求項３又は４に記載の解析方法。
【請求項７】
前記解析方法は、剛性行列Ｋと質量行列Ｍを直接計算する場合に、
前記解析対象物が３次元の中空体のモデルであるとき、
剛性行列の要素Ｋｉｊについて、
2006313400.tif 000047
の定義により、
剛性行列の要素Ｋｉｊを決定する式が、
ｉ≠ｊのとき、
2006313400.tif 000048
ｉ＝ｊのとき、
2006313400.tif 000049
のように得られる
ことを特徴とする請求項３又は４に記載の解析方法。
【請求項８】
前記解析方法は、剛性行列Ｋと質量行列Ｍを直接計算する場合に、
前記解析対象物が３次元の中空体のモデルであるとき、
質量行列の要素Ｍｉｊについて、
2006313400.tif 000050
の定義により、
質量行列の要素Ｍｉｊを決定する式が、
ｉ≠ｊのとき、
2006313400.tif 000051
ｉ＝ｊのとき、
2006313400.tif 000052
のように得られる
ことを特徴とする請求項３又は４に記載の解析方法。
【請求項９】
前記解析方法は、剛性行列Ｋと質量行列Ｍを直接計算する場合に、
前記解析対象物が３次元の中実体のモデルであるとき、
剛性行列の要素Ｋｉｊについて、
2006313400.tif 000053
の定義により、
剛性行列の要素Ｋｉｊを決定する式が、
ｉ≠ｊのとき、
2006313400.tif 000054
ｉ＝ｊのとき、
2006313400.tif 000055
のように得られる
ことを特徴とする請求項３又は４に記載の解析方法。
【請求項１０】
前記解析方法は、剛性行列Ｋと質量行列Ｍを直接計算する場合に、
前記解析対象物が３次元の中実体のモデルであるとき、
質量行列の要素Ｍｉｊについて、
2006313400.tif 000056
の定義により、
質量行列の要素Ｍｉｊを決定する式が、
ｉ≠ｊのとき、
2006313400.tif 000057
ｉ＝ｊのとき、
2006313400.tif 000058
のように得られる
ことを特徴とする請求項３又は４に記載の解析方法。
【請求項１１】
前記解析方法は、剛性行列Ｋと質量行列Ｍを直接計算することにより、与えられた領域のラプラシアンの固有値問題の固有値と固有関数を求める場合に、
固有値問題を、
2006313400.tif 000059
で与えたとき、
第ｋ番目の固有値λｋ
に対応する固有ベクトルｕｋの第ｉ成分を、第ｋ番目の固有関数φｋ の第ｉ番目の節点の値とし、隣接する節点間の値は線形補間によって定めることにより第ｋ番目の固有関数φｋ
の値を求め、
その固有関数φｋ
の値を用いて解析対象物を可視化する
ことを特徴とする請求項３ないし１０のいずれかに記載の解析方法。
【発明の詳細な説明】
【技術分野】
【０００１】
本発明は、有限要素法により解析対象物を解析する場合に、剛性行列Ｋと質量行列Ｍの計算時間を短縮し、精度よく固有関数を可視化することができる有限要素法直接計算プログラムおよび解析方法に関する。
【背景技術】
【０００２】
有限要素法では、連続体である物体を有限の形状と寸法を持つ要素の集合により近似して扱う。従来の有限要素法解析プログラムでは、解析対象物に適した数値解析を行うために、要素の形状、節点の座標、および各要素の局所的な節点番号を定め、これらの要素分割のデータを入力してから、剛性行列Ｋと質量行列Ｍを計算する必要があった。
【０００３】
例えば、解析対象物が２次元の平面のモデルである場合には、剛性行列Ｋと質量行列Ｍは、［数１］と［数２］の定義に従って計算させるプログラムになっている。
【０００４】
２次元の平面領域の三角形分割εに対する剛性行列Ｋ＝（Ｋｉｊ）は、φｉ を点Ｐｉ における基底関数とするとき、次のように定義される。
【数１】
2006313400.tif 000003
ただし、Ｇ（ε）は要素の和集合である。
【０００５】
２次元の平面領域の三角形分割εに対する質量行列Ｍ＝（Ｍｉｊ）は、φｉ を点Ｐｉ における基底関数とするとき、次のように定義される。
【数２】
2006313400.tif 000004
ただし、Ｇ（ε）は要素の和集合である。
【０００６】
また、これらの式は、３次元の中空体においても同様に定義され、解析対象物が３次元の中空体のモデルである場合には、剛性行列Ｋと質量行列Ｍは、［数３］と［数４］の定義に従って計算させるプログラムになっている。
【０００７】
３次元の中空体の三角形分割εに対する剛性行列Ｋ＝（Ｋｉｊ）は、φｉ を点Ｐｉ における基底関数とするとき、次のように定義される。
【数３】
2006313400.tif 000005
ただし、Ｇ（ε）は要素の和集合であり、〈∇φｉ ，∇φｊ 〉ｇ は中空体に誘導される計量ベクトルｇに関する内積であり、ｖｇ
は体積要素を表す。
【０００８】
３次元の中空体の三角形分割εに対する質量行列Ｍ＝（Ｍｉｊ）は、φｉ を点Ｐｉ における基底関数とするとき、次のように定義される。
【数４】
2006313400.tif 000006
ただし、Ｇ（ε）は要素の和集合であり、ｖgは体積要素を表す。
【０００９】
さらに、解析対象物が３次元の中実体のモデルである場合には、剛性行列Ｋと質量行列Ｍは、［数５］と［数６］の定義に従って計算させるプログラムになっている。
【００１０】
３次元の中実体の四面体分割εに対する剛性行列Ｋ＝（Ｋｉｊ）は、φｉ を点Ｐｉ における基底関数とするとき、次のように定義される。
【数５】
2006313400.tif 000007
ただし、Ｇ（ε）は要素の和集合である。
【００１１】
３次元の中実体の四面体分割εに対する質量行列Ｍ＝（Ｍｉｊ）は、φｉ を点Ｐｉ における基底関数とするとき、次のように定義される。
【数６】
2006313400.tif 000008
ただし、Ｇ（ε）は要素の和集合である。
【００１２】
上記従来技術では、計算速度を早くするために要素分割を工夫し、その他有限要素法の使用に当たって留意した場合であっても、計算結果を出力させるまでに、多くのステップの計算を必要とした。このため、多くの時間がかかる上に誤差が累積し、出力された結果は不正確なものになるという問題があった。
また、剛性行列Ｋと質量行列Ｍの計算に関し、直接剛性法により、要素ごとに積分計算することを全要素について行い、全体的な近似方程式を組み立てるという従来技術もあるが、データ入力に際し、各要素の節点番号付けの順序や、局所的な節点番号と全体的な節点番号の対応付けが必要になるなどの制約がある（例えば、非特許文献１参照。）。
さらに、有限要素法解析の精度を高めるために要素を小さくすると、要素数や節点数が増加し、データ入力とステップの計算に多くの作業と時間を要することになる。とくに、解析対象物が３次元中実体の複雑なモデルである場合には、要素数が数千万以上の規模になることがあり、剛性行列Ｋと質量行列Ｍの計算量が増加する。このため多くの時間がかかる上に誤差が累積し、出力された結果は不正確なものになる。さらに、剛性行列Ｋと質量行列Ｍの計算量が膨大になると、大きな記憶容量と高速に計算処理を行うことができる高性能の計算機が必要になり、一般に使用されている計算機では対応できなくなるという問題が生じる。
【００１３】
このため、従来から、剛性行列Ｋと質量行列Ｍの計算に関心が集まり、解析対象物の要素分割を工夫し、計算機における計算処理を分散化させるなどの改良が進められている。例えば、複数のプロセッサからなり、各プロセッサにメモリを持つようなマルチプロセッサ上での並列処理方式を採用し、各プロセッサが可能な限り、格納する剛性行列の成分、荷重ベクトルの成分を計算するのに必要となる要素剛性行列、応力を計算するように各プロセッサへの割り当てを行うことにより、バッファの容量を削減する有限要素法の処理方式（例えば、特許文献１参照。）が開示されている。また、要素分割の対象となる構造物の亀裂先端部分の作成位置と分割数を指定するだけで自動的に要素の頂点とその中間点の節点の座標データを作成し、前記亀裂先端部の基準要素に属する節点番号の集合を有する要素データを求め、前記基準要素の節点番号をもとにして他の要素の節点番号を求めて、その要素の節点番号の集合を有する要素データを次々に求める有限要素法解析用モデルの作成方法（例えば、特許文献２参照。）、有限要素法による解析対象物体の二次元又は三次元の有限要素データの作成方法において、要素の分割数の少ない方向から優先的に節点番号を付ける工程を備えた有限要素データの作成方法（例えば、特許文献３参照。）が開示されている。さらに、有限要素法の定型的手順により構造物の解析モデルから剛性行列および集中質量行列を生成し、これらの行列で表される一般固有値問題を標準的固有値問題Ａｙ＝λｙに変換し、集中質量行列が負の成分を含む場合、アーノルディの反復法によりＡを射影し行列Ｔ＝｛ｈｉｊ｝を得るが、このときの初期ベクトルを、対応する集中質量行列の成分が正であるような成分ｈに対しては実数、負であるような成分ｈに対しては純虚数となるように設定することにより、行列Ａを実ヘッセンベルグ行列Ｔに変換し、このヘッセンベルグ行列Ｔが与える標準的固有値問題の固有値をＱＲ法により求める方法およびシステム（例えば、特許文献４参照。）が開示されている。
【００１４】
【非特許文献１】菊地文雄著「有限要素法概説」１９８０年、Ｐ．５８−６７、株式    会社サイエンス社
【特許文献１】特開平０５−１０８６９６号公報
【特許文献２】特開平１１−２７２６４９号公報
【特許文献３】特開平１１−２８３０５０号公報
【特許文献４】特開２００１−２５６２１６号公報
【発明の開示】
【発明が解決しようとする課題】
【００１５】
しかしながら、特許文献１の開示技術では、解析対象物が３次元の中実体の複雑なモデルになると、要素数が数千万以上の規模になることがあるため、計算処理方法を改良しても、計算機では処理できない膨大な計算量になる可能性がある。また、特許文献２ないし４のいずれの開示技術でも、本発明のように剛性行列の要素Ｋｉｊを決定する式そのものを改良することにより、計算時間を大幅に短縮し、解析精度を向上させるという新たな提案はなされていない。
本発明は上記事情に鑑みてなされたものであり、要素分割のデータ入力に際し、データ入力作業から各要素の節点の局所番号付けする作業をなくし、高次要素を用いた計算であっても、データを与えるとほとんど誤差なく短時間に、剛性行列Ｋと質量行列Ｍを決定する式により計算して求め（以下、「直接計算」と呼ぶ。）、固有関数を精度良く可視化して出力することができる有限要素法直接計算プログラムおよび解析方法を提供することを目的とする。
【課題を解決するための手段】
【００１６】
上記課題を解決するための手段として、本発明は以下の特徴を有している。
本発明のプログラムは、有限要素法により解析対象物の数値モデルを作成し、数値解析を行い、解析結果を可視化するプログラムにおいて、有限要素法の基底関数から決まる剛性行列Ｋと質量行列Ｍを直接計算することにより、与えられた領域のラプラシアンの固有値問題の固有値と固有関数を求めることを特徴とする。
本発明によれば、剛性行列の要素Ｋｉｊを決定する式そのものを改良することにより、要素分割のデータ入力に際し、データ入力作業から各要素の節点の局所番号付けを省くことができ、データを与えるとほとんど誤差なく短時間に剛性行列Ｋと質量行列Ｍを直接計算して求め、固有値と固有関数を出力させることができる。本発明では、各要素の節点のつながりが把握できればよいため、従来必要とされていた各要素の節点番号付けの順序や、局所的な節点番号と全体的な節点番号の対応付けは不要になる。設定できる要素数と節点数の制限をなくしたことにより、１次元の線分、２次元の平面はもちろん、３次元の中空体、又は３次元の中実体の複雑な解析対象物の高次要素を用いた計算にも対応できるようにしている。
【００１７】
本発明のプログラムは、前記解析対象物が１次元モデルである場合に、解析対象物を１次元モデルによる解析方法より解析し、前記解析対象物が２次元モデルである場合に、解析対象物を１次元モデルによる解析方法、２次元モデルによる解析方法又はこれらの組み合わせにより解析し、前記解析対象物が３次元モデルである場合に、前記解析対象物を、１次元モデルによる解析方法、２次元モデルによる解析方法、３次元モデルによる解析方法又はこれらの任意の組み合わせにより解析することを特徴とする。
【００１８】
本発明の解析方法は、有限要素法により解析対象物の数値モデルを作成し、数値解析を行い、解析結果を可視化する解析方法において、有限要素法の基底関数から決まる剛性行列Ｋと質量行列Ｍを直接計算することにより、与えられた領域のラプラシアンの固有値問題の固有値と固有関数を求めることを特徴とする。
【００１９】
本発明の解析方法は、前記解析対象物が１次元モデルである場合に、解析対象物を１次元モデルによる解析方法より解析し、前記解析対象物が２次元モデルである場合に解析対象物を１次元モデルによる解析方法、２次元モデルによる解析方法又はこれらの組み合わせにより解析し、前記解析対象物が３次元モデルである場合に、前記解析対象物を、１次元モデルによる解析方法、２次元モデルによる解析方法、３次元モデルによる解析方法又はこれらの任意の組み合わせにより解析することを特徴とする。
【００２０】
本発明の解析方法は、剛性行列Ｋと質量行列Ｍを直接計算する場合に、前記解析対象物が２次元の平面のモデルであるとき、剛性行列の要素Ｋｉｊについて、
2006313400.tif 000009
の定義により、剛性行列の要素Ｋｉｊを決定する式が、
ｉ≠ｊのとき、
2006313400.tif 000010
ｉ＝ｊのとき、
2006313400.tif 000011
のように得られることを特徴とする。
本発明の剛性行列の要素Ｋｉｊを決定する式は、式そのものを改良したものであり、２次元の平面領域の解析対象物を三角形要素εに分割し、三角形の各線分の長さを分子に、それに対応する三角形の面積を分母にした量で定式化したものである。これらの式は、３次元の中空体においても、同様に適用することができる。
図８は、２次元平面領域および３次元中空体における、ｉ≠ｊのときの剛性行列Ｋｉｊを決定する式を説明するための概念図である。  節点Ｐｉ と節点Ｐｊ  が異なるときは、図８のようになる。上記ｉ≠ｊのときのＫｉｊの式の分母は対応する三角形の面積を表し、ａｒｅａ（ｅμ）は三つの節点Ｐｉ ，Ｐｊ ，Ｐkからなる三角形を、ａｒｅａ（ｅｖ） は三つの節点Ｐｉ ，Ｐｊ  ， Ｐｌ からなる三角形を表す。分子はそれぞれの辺の長さの２乗などを足したり引いたりしたものである。
図９は、２次元平面領域および３次元中空体における、ｉ&shy;＝ｊのときの剛性行列Ｋｉｉを決定する式を説明するための概念図である。  節点Ｐｉ と節点Ｐｊ  が同じときは、図９のようになる。節点Ｐｉ
に接続する節点と三角形が書かれている。上記ｉ&shy;＝ｊのときのＫｉｉの式は、各三角形のＰｉ  を含まない辺の長さを２乗し、三角形の面積で割った量の和で与えられる。
従って、局所的な節点番号と全体的な節点番号を対応づけは不要となる。３次元の中空体の場合も同様である。
【００２１】
本発明の解析方法は、剛性行列Ｋと質量行列Ｍを直接計算する場合に、前記解析対象物が２次元の平面のモデルであるとき、質量行列の要素Ｍｉｊについて、
2006313400.tif 000012
の定義により、質量行列の要素Ｍｉｊを決定する式が、
ｉ≠ｊのとき、
2006313400.tif 000013
ｉ＝ｊのとき、
2006313400.tif 000014
のように得られることを特徴とする。
本発明の質量行列の要素Ｍｉｊを決定する式は、３次元の中空体においても、同様に適用することができる。
【００２２】
本発明の解析方法は、剛性行列Ｋと質量行列Ｍを直接計算する場合に、前記解析対象物が３次元の中空体のモデルであるとき、剛性行列の要素Ｋｉｊについて、
2006313400.tif 000015
の定義により、剛性行列の要素Ｋｉｊを決定する式が、
ｉ≠ｊのとき、
2006313400.tif 000016
ｉ＝ｊのとき、
2006313400.tif 000017
のように得られることを特徴とする。
【００２３】
本発明の解析方法は、剛性行列Ｋと質量行列Ｍを直接計算する場合に、前記解析対象物が３次元の中空体のモデルであるとき、質量行列の要素Ｍｉｊについて、
2006313400.tif 000018
の定義により、質量行列の要素Ｍｉｊを決定する式が、
ｉ≠ｊのとき、
2006313400.tif 000019
ｉ＝ｊのとき、
2006313400.tif 000020
のように得られることを特徴とする。
【００２４】
本発明の解析方法は、剛性行列Ｋと質量行列Ｍを直接計算する場合に、前記解析対象物が３次元の中実体のモデルであるとき、剛性行列の要素Ｋｉｊについて、
2006313400.tif 000021
の定義により、剛性行列の要素Ｋｉｊを決定する式が、
ｉ≠ｊのとき、
2006313400.tif 000022
ｉ＝ｊのとき、
2006313400.tif 000023
のように得られることを特徴とする。
本発明の剛性行列の要素Ｋｉｊを決定する式は、式そのものを改良したものであり、３次元の中実体の解析対象物を四面体要素ｅに分割し、四面体の各線分の長さを分子に、それに対応する四面体の体積を分母にし、四面体要素ｅすべてにわたる和として定式化したものである。
【００２５】
本発明の解析方法は、剛性行列Ｋと質量行列Ｍを直接計算する場合に、前記解析対象物が３次元の中実体のモデルであるとき、質量行列の要素Ｍｉｊについて、
2006313400.tif 000024
の定義により、質量行列の要素Ｍｉｊを決定する式が、
ｉ≠ｊのとき、
2006313400.tif 000025
ｉ＝ｊのとき、
2006313400.tif 000026
のように得られることを特徴とする。
本発明の質量行列の要素Ｍｉｊを決定する式は、四面体要素ｅすべてにわたる和として定式化したものであり、Ｖｏｌ（ｅ）は四面体の体積を表す。
【００２６】
本発明の解析方法は、剛性行列Ｋと質量行列Ｍを直接計算することにより、与えられた領域のラプラシアンの固有値問題の固有値と固有関数を求める場合に、固有値問題を、
2006313400.tif 000027
で与えたとき、第ｋ番目の固有値λｋ
に対応する固有ベクトルｕｋの第ｉ成分を、第ｋ番目の固有関数φｋ の第ｉ番目の節点の値とし、隣接する節点間の値は線形補間によって定めることにより第ｋ番目の固有関数φｋ
の値を求め、その固有関数φｋ の値を用いて解析対象物を可視化することを特徴とする。
本発明の可視化方法は、剛性行列Ｋと質量行列Ｍを直接計算することにより、固有空間に所属する第ｋ番目の固有関数φｋの値を用いて解析対象物を可視化して短時間に精度良く出力することができる。
図１０は、解析対象物が２次元平面および３次元中空体のモデルであるときの、隣接する節点間の値を線形補間によって定める方法の説明をするための概念図である。
【発明の効果】
【００２７】
上記課題を解決するための手段により、本発明によれば、有限要素法により対象物を解析する場合に、要素分割のデータ入力に際し、データ入力作業から各要素の節点の局所番号付けを省くことができた。また、従来必要とされていた各要素の節点番号付けの順序や局所的な節点番号と全体的な節点番号の対応付けを不要にしたことによって、設定できる要素数と節点数の制限をなくすことができた。本発明によれば、１次元の線分、２次元の平面はもちろん、３次元の中空体、又は３次元の中実体の複雑な解析対象物を解析する場合であっても、データを与えると、従来技術に比べて、ほとんど誤差なく短時間に、剛性行列Ｋと質量行列Ｍを直接計算して求めることができた。
また、剛性行列Ｋと質量行列Ｍを直接計算した結果を用いて、ディリクレ境界値問題とノイマン境界値固有問題の固有値と固有関数、および高次の固有値と固有関数を容易かつ精度良く可視化して出力することができることから、定量的で具体的な解析をすることができ、取り扱いが容易で実用的な有限要素法直接計算プログラムおよび解析方法を提供することができた。
【発明を実施するための最良の形態】
【００２８】
本発明を実施するための最良の形態を、図１から図３に基づいて説明する。図１から図３では、連続体である３次元の中実体の解析対象物を、有限の形状と寸法を持つ四面体要素ｅに分割し、各要素の集合により近似して扱う。
ただし、これらは一実施形態にすぎず、本発明の特許請求の範囲を限定するものではない。
【００２９】
図１は、有限要素法直接計算プログラムの解析方法の手順を示したフローチャートである。ステップＳ０１〜Ｓ１４により剛性行列Ｋと質量行列Ｍを決定して出力し、ステップＳ１５〜Ｓ１６により固有値と固有関数を求めて出力する手順により、解析する。
【００３０】
図２は、図４における領域Ω１ の四面体要素の３５個の節点の座標を与える、節点の座標ファイルの例である。節点とは各四面体要素の頂点をいう。図２の節点の座標ファイルは、前もって準備しておく。節点数を１行目に入力し、２行目以降は節点の座標を１行ずつ入力する。この順番が節点番号となるので、節点番号を与える必要はない。座標は、左からｘ座標、ｙ座標およびｚ座標とする。例えば、１行目は節点数が３５であり、２行目は節点番号１のｘ座標が０．５、ｙ座標が０．５、ｚ座標が０．５であることを意味する。
【００３１】
図３は、図４における領域Ω１
の９６個の四面体要素の節点の連結関係を与える、連結ファイルの例である。連結ファイルとは、各要素をなす節点の番号のファイルをいう。節点番号は、３５個の節点Ｐ１
〜Ｐ35 の添字を用いて１〜３５のように表す。図３の連結ファイルは、前もって準備しておく。図３の連結ファイルは、連結関係にある４点の節点番号の組の総数を１行目に入力し、２行目以降は連結関係にある４点の節点番号の組を１行ずつ入力する。この場合、四面体をなす節点番号が分かればよく、時計、反時計まわりのどちらの順番で入力してもよい。図３の２行目は節点番号１６，３５，２３，２からなる四面体の要素を表している。
【００３２】
本発明の剛性行列の要素Ｋｉｊと質量行列の要素Ｍｉｊを決定する式の計算方法について述べる。
【００３３】
ｉ≠ｊのとき、剛性行列Ｋと質量行列Ｍの（ｉ，ｊ）成分は、［数７］、［数８］および［数９］により計算される。ここで、ｉ、ｊ、ｋ、ｌは自然数である。
【００３４】
節点Ｐの座標を（ｘ（Ｐ），ｙ（Ｐ），ｚ（Ｐ））で表したとき、始点Ｐ、終点ＱのベクトルＰＱと始点Ｒ、終点ＳのベクトルＲＳの内積は次のように表される。
【数７】
2006313400.tif 000028
【００３５】
【数８】
2006313400.tif 000029
ここで、右辺の和は辺Ｐｉ
Ｐｊ を含む四面体ｅすべてにわたる和であり、四面体ｅの４つの節点Ｐｉ ，Ｐｊ ，Ｐk （ｅ），Ｐl （ｅ）で表す。また、ＰｉＰl （ｅ）等は、始点をＰｉ 、終点をＰl （ｅ）とするベクトルを表し、（ＰｉＰl （ｅ），Ｐｊ Ｐl （ｅ））は２つのベクトルＰｉＰl （ｅ）とＰｊ Ｐl （ｅ）との、［数７］における内積を表す。Ｖｏｌ（ｅ）は四面体ｅの体積を表している。
【００３６】
【数９】
2006313400.tif 000030
ここで、右辺の和は辺Ｐｉ
Ｐｊ を含む四面体ｅすべてにわたる和であり、Ｖｏｌ（ｅ）は四面体ｅの体積を表している。
【００３７】
ｉ＝ｊのとき、剛性行列Ｋと質量行列Ｍの（ｉ，ｉ）成分は、［数７］、［数１０］および［数１１］により計算される。ここで、ｉ、ｊ、ｋ、ｌは自然数である。
【００３８】
【数１０】
2006313400.tif 000031
ここで、右辺の和は節点Ｐｉ
を含む四面体ｅすべてにわたる和である。Ｖｏｌ（ｅ）は四面体ｅの体積を表す。
【００３９】
【数１１】
2006313400.tif 000032
ここで、右辺の和は節点Ｐｉ
を含む四面体ｅすべてにわたる和である。Ｖｏｌ（ｅ）は四面体ｅの体積を表す。
【００４０】
以下、図１に示すステップＳ０１〜Ｓ１４により剛性行列Ｋと質量行列Ｍを決定して出力する解析方法について説明する。
ステップＳ０１では、節点の座標ファイルを読み込む。節点の座標ファイルは、図２のように前もって準備しておく。
ステップＳ０２では、連結ファイルを読み込む。連結ファイルは、図３のように前もって準備しておく。
ステップＳ０３では、Ｌ×Ｌ行列のすべての成分を零とする初期化を行っている。Ｌ×Ｌ行列は、解決すべき問題に対し、必要な行列とする。ただし、Ｌは１≦Ｌ≦Ｎの定数とし、Ｌはプログラムの中で与える。ここで、Ｎは自然数であり、節点数を表す。
【００４１】
ステップＳ０４〜Ｓ０７では、ｉ≠ｊのときの、剛性行列の要素Ｋｉｊと質量行列の要素Ｍｉｊを決定する式の計算をしている。
ステップＳ０４では、ステップＳ０２で読み込んだ節点Ｐｉ
およびＰｊ（ｉ＜ｊ）をもつ四面体ｅの残りの節点がＰｋ ，Ｐｌ  ならば、ステップＳ０１で読み込んだ各節点の座標（ｘ（ｉ），ｙ（ｉ），ｚ（ｉ））、（ｘ（ｊ），ｙ（ｊ），ｚ（ｊ））、（ｘ（ｋ），ｙ（ｋ），ｚ（ｋ））、（ｘ（ｌ）,ｙ（ｌ）,ｚ（ｌ））を用いて、内積（ＰｉＰｌ，ＰｊＰｌ）、（ＰｉＰｋ，ＰｊＰｋ）、（ＰｉＰｌ，ＰｊＰｋ）、（ＰｉＰｋ，ＰｊＰｌ）を計算する。
ステップＳ０５では、四面体ｅの体積Ｖｏｌ（ｅ）を求める。
ステップＳ０６では、ＡＫ（ｉ，ｊ）、ＡＭ（ｉ，ｊ）を求める。
［数８］より
2006313400.tif 000033
の値をＡＫ（ｉ，ｊ）、
［数７］より
2006313400.tif 000034
の値をＡＭ（ｉ，ｊ）とする。
ステップＳ０７では、ステップＳ０６と同様に、四面体ｅの辺Ｐｉ
Ｐｊ 以外についても計算する。具体的には、
ＡＫ（ｉ，ｋ）、ＡＫ（ｉ，ｌ）、ＡＫ（ｊ，ｋ）、ＡＫ（ｊ，ｌ）、ＡＫ（ｋ，ｌ）、
ＡＭ（ｉ，ｋ）、ＡＭ（ｉ，ｌ）、ＡＭ（ｊ，ｋ）、ＡＭ（ｊ，ｌ）、ＡＭ（ｋ，ｌ）
である。
【００４２】
ステップＳ０８〜Ｓ１１では、ｉ＝ｊのときの、剛性行列の要素Ｋｉｊと質量行列の要素Ｍｉｊを決定する式の計算をしている。
ステップＳ０８では、節点Ｐｉ
を含む四面体ｅの残りの節点がＰｊ ，Ｐk  ，Ｐl ならばステップＳ０１で読み込んだ各節点の座標（ｘ（ｉ），ｙ（ｉ），ｚ（ｉ））、（ｘ（ｊ），ｙ（ｊ），ｚ（ｊ））、（ｘ（ｋ），ｙ（ｋ），ｚ（ｋ））（ｘ（ｌ），ｙ（ｌ），ｚ（ｌ））を用いて、内積（Ｐｊ  Ｐｌ，Ｐｊ
Ｐｋ  ）と‖ＰｊＰｌ‖２、‖ＰｊＰｋ‖２
を計算する。
ステップＳ０９では、四面体ｅの体積Ｖｏｌ（ｅ）を求める。
ステップＳ１０では、ＡＫ（ｉ，ｉ）、ＡＭ（ｉ，ｉ）を求める。
［数８］より
2006313400.tif 000035
の値をＡＫ（ｉ，ｉ）、
［数９］より
2006313400.tif 000036
の値をＡＭ（ｉ，ｉ）とする。
ステップＳ１１では、ステップＳ１０と同様に、四面体ｅの点Ｐｉ
以外についても計算する。具体的には、
ＡＫ（ｊ，ｊ）、ＡＫ（ｋ，ｋ）、ＡＫ（ｌ，ｌ）、ＡＭ（ｊ，ｊ）、ＡＭ（ｋ，ｋ）、
ＡＭ（ｌ，ｌ）
である。
【００４３】
ステップＳ１２では、すべての四面体要素に対してＳ０４〜Ｓ１１で行ったことを繰り返す。そして、各要素を加えた和集合が、行列ＫおよびＭの要素となる。
ステップＳ１３では、ステップＳ１２で求めた行列ＫおよびＭの成分を用いて、（ｊ，ｉ）成分（ｉ＞ｊ）を求める。連結関係のない成分はすべて零であり、ステップＳ０３で初期化した結果をそのまま用いる。以上で、行列ＫおよびＭの成分が決定できる。
ステップＳ１４では、行列ＫおよびＭのすべての成分を出力する。このステップは省略してもよい。
【００４４】
次に、ステップＳ１５〜Ｓ１６により固有値と固有関数を求めて出力することにより、解析対象物を可視化する解析方法について説明する。
ステップＳ１５では、［数１２］で表される行列の一般固有値問題を解く。ここでは既に使用されている行列の一般固有値問題を解くプログラムを組み込んで解いているが、ステップＳ１４で出力してある行列ＫおよびＭの成分を用いて別のプログラムを使用しても構わない。
【数１２】
2006313400.tif 000037
【００４５】
【数１３】
2006313400.tif 000038
ただし、Δはラプラシアンであって、φのラプラシアンは次のように表される。
2006313400.tif 000039
【００４６】
［数１２］における第ｋ番目の固有値を、［数１３］における第ｋ番目の固有値λｋと定める。また、［数１２］の第ｋ番目の固有値λｋ
に対応する固有ベクトルｕｋから、［数１３］における第ｋ番目の固有関数φｋを得るには、固有ベクトルｕｋの第ｉ成分を固有関数φｋ
の第ｉ番目の節点の値とし、隣接する節点間の値を線形補間によって固有関数φｋ の値を定める。
すなわち、（ｘ，ｙ，ｚ）∈四面体Ｐｉ Ｐｊ Ｐk Ｐｌ において、
φｋ
（ｘ，ｙ，ｚ）＝ａ×（ｕの第ｉ成分）＋ｂ×（ｕの第ｊ成分）＋ｃ×（ｕの第ｋ成分）＋ｄ×（ｕの第ｌ成分）
によって与える。
ただし、０≦ａ，ｂ，ｃ，ｄ≦１、ａ＋ｂ＋ｃ＋ｄ≦１である。
【００４７】
ステップＳ１６では、［数１３］における第ｋ番目の固有値λｋ
、固有関数φｋ の出力を行う。ただし、ｋは１≦ｋ≦Ｌの定数とする。
以上で、図１に示すステップＳ０１〜Ｓ１９の一連の処理を終了する。
【００４８】
次に、本発明の剛性行列Ｋの計算時間について、従来の計算方法と比較するために、２次元の場合について説明する。
従来の計算方法については、非特許文献１に記載の剛性行列の計算法に従って演算回数を調べる。三角形要素の１つをＰ１ （ｘ１ ，ｙ１ ），Ｐ２
（ｘ２ ，ｙ２ ），Ｐ３ （ｘ３ ，ｙ３ ）とし、反時計回りに並んでいるものとする。従来の方法では、３点がＰ１ ，Ｐ２ ，Ｐ３ の順に反時計回りに並んでいる必要があるが、本発明の方法では向きを考慮する必要はない。
【００４９】
従来の計算方法によると、
（１）  面積Ｓの２倍に等しい行列式Ｄの計算
Ｄ＝２×Ｓ＝ｘ１×（ｙ２
−ｙ３ ）＋ｘ２×（ｙ３ −ｙ１ ）＋ｘ３×（ｙ１
−ｙ２ ）
と計算できるので、加算２回、減算３回、乗算３回の合計８回の演算が必要となる。
（２）  面積座標での近似関数の係数の計算
2006313400.tif 000040
は減算１回、除算１回でのそれぞれで行うので、合計１２回の演算が必要となる。
（３）Ｓ＝Ｄ／２より、除算１回が必要となる。
（４）要素係数行列の計算
（２）で求めた係数を用いて、Ａｉｊ＝Ｓ×（ｂｊ×ｂｉ＋ｃｊ×ｃｉ）を計算する。
すなわち、Ａ11 ，Ａ22，Ａ33 ，Ａ12（＝Ａ21），Ａ23（＝Ａ32），Ａ31（＝Ａ13）のそれぞれで乗算３回、加算１回を６回行うので、合計２４回の演算が必要となる。
（５）９つの要素を加える
合計９回の演算が必要となる。
以上より、要素がＮならば、（８＋１２＋１＋２４＋９）Ｎ＝５４Ｎの演算回数が必要となる。
【００５０】
本発明の計算方法によると、
（１）面積Ｓの２倍に等しい行列式Ｄの計算
従来の方法と同様に８回の演算が必要となる。
（２）線分の２乗の計算
ＰｉＰj２ ＝（ｘｉ
−ｘｊ ）２＋（ｙｉ −ｙｊ）２
と計算できるので、加算１回、減算２回、乗算２回となり、
ＰｊＰｋ２，ＰｋＰｉ２のそれぞれで行うので、合計１５回の演算が必要となる。
（３）Ｓ１ ＝−０．２５／Ｄ、Ｓ２ ＝０．５／Ｄより、除算２回が必要となる。
（４）要素係数行列の計算
Ａｉｊ＝Ｓ１×（ＰｉＰk２＋ＰｊＰk２−ＰｉＰj２）より、加算１回、減算１回、乗算１回
Ａｉｉ＝Ｓ２×ＰｊＰk２より、乗算１回を用いて、
Ａ11 ，Ａ22，Ａ33 ，Ａ12（＝Ａ21），Ａ23（＝Ａ32），Ａ31（＝Ａ13）のそれぞれを計算するので、合計１２回の演算が必要となる。
（５）９つの要素を加える
合計９回の演算が必要となる。
以上より、要素がＮならば、（８＋１５＋２＋１２＋９）Ｎ＝４６Ｎの演算回数が必要となる。
【００５１】
従って、従来の計算方法と本発明の計算方法とを比較すると、
（１−４６Ｎ／５４Ｎ）×１００＝１４．８％
となり、およそ１５％の計算時間の短縮が見込まれる。
１５％の計算時間の短縮により、従来の計算方法では３ギガバイトの記憶容量が必要であったものが、本発明の計算方法では２．５ギガバイトの記憶容量で足りることになる。経費にすると、およそ５、６万円の削減が可能になる。
また、要素数５０００（節点２６０１）の場合、従来の方法では０．００１３８秒という結果が得られている。このことから、要素の数が数千万以上の規模となった場合、本発明の計算方法を適用すれば、計算時間の短縮に寄与することが可能になる。
【００５２】
また、３次元の場合、中空体、中実体のいずれも、剛性行列Ｋの計算時間について厳密に計算ステップを記載してある従来の計算方法は見られないようであるが、３次元の中実体の複雑な解析対象物になると、従来の計算方法では計算に１日ほど要することが見込まれる。
しかしながら、本発明の計算方法を適用すれば、３次元の中空体、又は３次元の中実体の複雑な解析対象物を解析する場合であっても、データを与えると、従来の計算方法に比べて、ほとんど誤差なく短時間に、剛性行列Ｋと質量行列Ｍを直接計算することにより解析することができる。
例えば、節点２３７９、要素数８８７１の直方体の場合には、
（１） 
行列ＫおよびＭの決定に、８４秒を要し、
（２） 
固有値、固有ベクトルの決定に、４４０秒を要した。
以上より、８４秒＋４４０秒＝５２４秒＝８分４４秒の計算時間で処理できる。
ただし、（１）（２）いずれもデータの入出力等の時間が含まれており、計算のみに比べれば時間は余分にかかっている。
【実施例１】
【００５３】
［実施例１］は、３次元の中実体の例である。
図４は、領域Ω１
における立方体領域を四面体要素に分割したときの節点および要素の状態を示した図である。
領域Ω１
＝｛０＜ｘ＜１，０＜ｙ＜１，０＜ｚ＜１｝における
固有値問題：Δｕ＝λｕ    ｕ＝ｕ（ｘ，ｙ，ｚ）；（ｘ，ｙ，ｚ）∈Ω１  、
ディリクレ境界条件：ｕ（ｘ，ｙ，ｚ）＝０    （ｘ，ｙ，ｚ）∈∂Ω１   
を図４の四面体要素に分割することにより求める。この四面体分割の節点の座標と連結情報は図２および図３で与えるものとする。
図５は、本発明の解析方法を図４に適用した場合の行列ＫおよびＭの各成分の一覧である。
図６は、本発明の解析方法を図４に適用した場合の固有値および固有ベクトルの数値解である。
図７は、本発明の解析方法を図４に適用した場合の第１〜第５固有関数を可視化したものである。右図は左図の内部を取り出したものである。
【実施例２】
【００５４】
［実施例２］は、２次元平面領域の例である。
図１１は、領域Ω１の節点および要素の状態を示した図である。
領域Ω１
＝｛０＜ｘ＜１，０＜ｙ＜１｝における
固有値問題：Δｕ＝λｕ    ｕ＝ｕ（ｘ，ｙ）；（ｘ，ｙ）∈Ω１  、
ディリクレ境界条件：ｕ（ｘ，ｙ）＝０    （ｘ，ｙ）∈∂Ω１   
を図１１の三角形要素に分割することにより求める。
【実施例３】
【００５５】
［実施例３］は、２次元平面領域の例である。
領域Ω２ ついて三角形分割を行い、Ω２ の固有値問題およびディリクレ境界条件及びノイマン境界条件に対する固有値問題に対する第１９の固有関数を計算し、可視化した結果を図１５と１６に示す。
【実施例４】
【００５６】
［実施例４］は、３次元の中空体の例である。
３次元中空体の曲面における固有関数の結果を図１７〜１９に示す。
固有値については、Ωを上面Ω＋ と下面Ω− の２重に三角形分割して得られる、潰れた曲面とみなし、固有値を出力する。固有関数については、出力された固有関数の上面Ω＋と下面Ω−
の節点の値の差を出力する。
【実施例５】
【００５７】
熱方程式に対して、温度分布を解析した実施例を示す。１ｍ四方の領域の頂点を、それぞれ（０，０）、（１，０）、（１，１）、（０，１）とおく。初期条件として、ｕ（ｔ，ｘ，ｙ）＝ｕ（０，０．５，０．２５）＝１０４を与える。この領域の時刻ｔ＝０．０２秒のディリクレ境界条件とノイマン境界条件の温度分布を図２０〜２３にそれぞれ示す。
【実施例６】
【００５８】
波動方程式に対して、波の動きを解析した実施例を示す。１ｍ四方の領域の頂点をそれぞれ（０，０）、（１，０）、（１，１）、（０，１）とおく。初期条件として、ｕ（ｔ，ｘ，ｙ）＝ｕ（０，０．０１，０．０１）＝ｕ（０，０．９９，０．９９）＝１０４  を与える。
図２４は、この領域の時刻ｔ＝０．０４秒のディリクレ境界条件の時の波の動きを表したものである。
図２５は、この領域の時刻ｔ＝０．０３秒のノイマン境界条件の時の波の動きを表したものである。
【図面の簡単な説明】
【００５９】
【図１】有限要素法直接計算プログラムの解析方法の手順を示したフローチャートで    ある。
【図２】図４における領域Ω１の四面体要素の３５個の節点の座標を与える、節点の    座標ファイルの例である。
【図３】図４における領域Ω１の９６個の四面体要素の節点の連結関係を与える、連    結ファイルの例である。
【図４】領域Ω１における立方体領域を四面体要素に分割したときの節点および要素    の状態を示した図である。
【図５】本発明の解析方法を図４に適用した場合の行列ＫおよびＭの各成分の一覧である。
【図６】本発明の解析方法を図４に適用した場合の固有値および固有ベクトルの数値解である。
【図７】本発明の解析方法を図４に適用した場合の第１〜第５固有関数を可視化したものである。
【図８】２次元平面領域および３次元中空体における、ｉ&shy;＝ｊのときの剛性行列Ｋｉｊを決定する式を説明するための概念図である。
【図９】２次元平面領域および３次元中空体における、ｉ≠ｊのときの剛性行列Ｋｉｉを決定する式を説明するための概念図である。
【図１０】解析対象物が２次元平面および３次元中空体のモデルであるときの、隣接する節点間の値を線形補間によって定める方法の説明をするための概念図である。
【図１１】領域Ω１ の節点および要素の状態を示した図である。
【図１２】本発明の解析方法を図１１に適用した場合の行列ＫおよびＭの各成分の一覧である。
【図１３】本発明の解析方法を図１１に適用した場合の固有値および固有ベクトルの数値解である。
【図１４】領域Ω２ の節点および要素の状態を示した図である。
【図１５】本発明の解析方法を図１４に適用した場合のディリクレ境界条件における第１９番目の固有関数を可視化した図である。
【図１６】本発明の解析方法を図１４に適用した場合のノイマン境界条件における第１９番目の固有関数を可視化した図である。
【図１７】本発明の解析方法を球面に適用した場合の第２４番目の固有関数を可視化した図である。
【図１８】本発明の解析方法をトーラスに適用した場合の第２７番目の固有関数を可視化した図である。
【図１９】本発明の解析方法を楕円面に適用した場合の第２６番目の固有関数を可視化した図である。
【図２０】本発明の解析方法を正方形領域に適用した場合のディリクレ境界条件における時刻ｔ＝０．０２での温度分布の正面図である。
【図２１】本発明の解析方法を正方形領域に適用した場合のディリクレ境界条件における時刻ｔ＝０．０２での温度分布の立体図である。
【図２２】本発明の解析方法を正方形領域に適用した場合のノイマン境界条件における時刻ｔ＝０．０２での温度分布の正面図である。
【図２３】本発明の解析方法を正方形領域に適用した場合のノイマン境界条件における時刻ｔ＝０．０２での温度分布の立体図である。
【図２４】本発明の解析方法を正方形領域に適用した場合のディリクレ境界条件における時刻ｔ＝０．０４での波の動きを表す正面図である。
【図２５】本発明の解析方法を正方形領域に適用した場合のノイマン境界条件における時刻ｔ＝０．０３での波の動きを表す立体図である。
【図１】
2006313400.tif 000060
【図２】
2006313400.tif 000061
【図３】
2006313400.tif 000062
【図４】
2006313400.tif 000063
【図５】
2006313400.tif 000064
【図６】
2006313400.tif 000065
【図７】
2006313400.tif 000066
【図８】
2006313400.tif 000067
【図９】
2006313400.tif 000068
【図１０】
2006313400.tif 000069
【図１１】
2006313400.tif 000070
【図１２】
2006313400.tif 000071
【図１３】
2006313400.tif 000072
【図１４】
2006313400.tif 000073
【図１５】
2006313400.tif 000074
【図１６】
2006313400.tif 000075
【図１７】
2006313400.tif 000076
【図１８】
2006313400.tif 000077
【図１９】
2006313400.tif 000078
【図２０】
2006313400.tif 000079
【図２１】
2006313400.tif 000080
【図２２】
2006313400.tif 000081
【図２３】
2006313400.tif 000082
【図２４】
2006313400.tif 000083
【図２５】
2006313400.tif 000084