JP2006313400A - 有限要素法直接計算プログラムおよび解析方法 - Google Patents
有限要素法直接計算プログラムおよび解析方法 Download PDFInfo
- Publication number
- JP2006313400A JP2006313400A JP2005134797A JP2005134797A JP2006313400A JP 2006313400 A JP2006313400 A JP 2006313400A JP 2005134797 A JP2005134797 A JP 2005134797A JP 2005134797 A JP2005134797 A JP 2005134797A JP 2006313400 A JP2006313400 A JP 2006313400A
- Authority
- JP
- Japan
- Prior art keywords
- analysis
- analysis method
- matrix
- dimensional
- dimensional model
- Prior art date
- Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
- Granted
Links
Landscapes
- Complex Calculations (AREA)
Abstract
【課題】 要素分割のデータ入力に際し、データ入力作業から各要素の節点の局所番号付けする作業をなくし、高次要素を用いた計算であっても、データを与えるとほとんど誤差なく短時間に、剛性行列Ｋと質量行列Ｍを直接計算して求め、固有関数を精度良く可視化して出力することができる有限要素法直接計算プログラムおよび解析方法を提供することを目的とする。
【解決手段】 有限要素法により解析対象物の数値モデルを作成し、数値解析を行い、解析結果を可視化するプログラムの解析方法において、有限要素法の基底関数から決まる剛性行列Ｋと質量行列Ｍを直接計算することにより、与えられた領域のラプラシアンの固有値問題の固有値と固有関数を求める。
【選択図】 図１
【解決手段】 有限要素法により解析対象物の数値モデルを作成し、数値解析を行い、解析結果を可視化するプログラムの解析方法において、有限要素法の基底関数から決まる剛性行列Ｋと質量行列Ｍを直接計算することにより、与えられた領域のラプラシアンの固有値問題の固有値と固有関数を求める。
【選択図】 図１
Description
本発明は、有限要素法により解析対象物を解析する場合に、剛性行列Ｋと質量行列Ｍの計算時間を短縮し、精度よく固有関数を可視化することができる有限要素法直接計算プログラムおよび解析方法に関する。
有限要素法では、連続体である物体を有限の形状と寸法を持つ要素の集合により近似して扱う。従来の有限要素法解析プログラムでは、解析対象物に適した数値解析を行うために、要素の形状、節点の座標、および各要素の局所的な節点番号を定め、これらの要素分割のデータを入力してから、剛性行列Ｋと質量行列Ｍを計算する必要があった。
例えば、解析対象物が２次元の平面のモデルである場合には、剛性行列Ｋと質量行列Ｍは、［数１］と［数２］の定義に従って計算させるプログラムになっている。
２次元の平面領域の三角形分割εに対する剛性行列Ｋ＝（Ｋｉｊ）は、φｉ を点Ｐｉ における基底関数とするとき、次のように定義される。
ただし、Ｇ（ε）は要素の和集合である。
ただし、Ｇ（ε）は要素の和集合である。
２次元の平面領域の三角形分割εに対する質量行列Ｍ＝（Ｍｉｊ）は、φｉ を点Ｐｉ における基底関数とするとき、次のように定義される。
ただし、Ｇ（ε）は要素の和集合である。
ただし、Ｇ（ε）は要素の和集合である。
また、これらの式は、３次元の中空体においても同様に定義され、解析対象物が３次元の中空体のモデルである場合には、剛性行列Ｋと質量行列Ｍは、［数３］と［数４］の定義に従って計算させるプログラムになっている。
３次元の中空体の三角形分割εに対する剛性行列Ｋ＝（Ｋｉｊ）は、φｉ を点Ｐｉ における基底関数とするとき、次のように定義される。
ただし、Ｇ（ε）は要素の和集合であり、〈∇φｉ ，∇φｊ 〉ｇ は中空体に誘導される計量ベクトルｇに関する内積であり、ｖｇ
は体積要素を表す。
ただし、Ｇ（ε）は要素の和集合であり、〈∇φｉ ，∇φｊ 〉ｇ は中空体に誘導される計量ベクトルｇに関する内積であり、ｖｇ
は体積要素を表す。
３次元の中空体の三角形分割εに対する質量行列Ｍ＝（Ｍｉｊ）は、φｉ を点Ｐｉ における基底関数とするとき、次のように定義される。
ただし、Ｇ（ε）は要素の和集合であり、ｖgは体積要素を表す。
ただし、Ｇ（ε）は要素の和集合であり、ｖgは体積要素を表す。
さらに、解析対象物が３次元の中実体のモデルである場合には、剛性行列Ｋと質量行列Ｍは、［数５］と［数６］の定義に従って計算させるプログラムになっている。
３次元の中実体の四面体分割εに対する剛性行列Ｋ＝（Ｋｉｊ）は、φｉ を点Ｐｉ における基底関数とするとき、次のように定義される。
ただし、Ｇ（ε）は要素の和集合である。
ただし、Ｇ（ε）は要素の和集合である。
３次元の中実体の四面体分割εに対する質量行列Ｍ＝（Ｍｉｊ）は、φｉ を点Ｐｉ における基底関数とするとき、次のように定義される。
ただし、Ｇ（ε）は要素の和集合である。
ただし、Ｇ（ε）は要素の和集合である。
上記従来技術では、計算速度を早くするために要素分割を工夫し、その他有限要素法の使用に当たって留意した場合であっても、計算結果を出力させるまでに、多くのステップの計算を必要とした。このため、多くの時間がかかる上に誤差が累積し、出力された結果は不正確なものになるという問題があった。
また、剛性行列Ｋと質量行列Ｍの計算に関し、直接剛性法により、要素ごとに積分計算することを全要素について行い、全体的な近似方程式を組み立てるという従来技術もあるが、データ入力に際し、各要素の節点番号付けの順序や、局所的な節点番号と全体的な節点番号の対応付けが必要になるなどの制約がある（例えば、非特許文献１参照。）。
さらに、有限要素法解析の精度を高めるために要素を小さくすると、要素数や節点数が増加し、データ入力とステップの計算に多くの作業と時間を要することになる。とくに、解析対象物が３次元中実体の複雑なモデルである場合には、要素数が数千万以上の規模になることがあり、剛性行列Ｋと質量行列Ｍの計算量が増加する。このため多くの時間がかかる上に誤差が累積し、出力された結果は不正確なものになる。さらに、剛性行列Ｋと質量行列Ｍの計算量が膨大になると、大きな記憶容量と高速に計算処理を行うことができる高性能の計算機が必要になり、一般に使用されている計算機では対応できなくなるという問題が生じる。
また、剛性行列Ｋと質量行列Ｍの計算に関し、直接剛性法により、要素ごとに積分計算することを全要素について行い、全体的な近似方程式を組み立てるという従来技術もあるが、データ入力に際し、各要素の節点番号付けの順序や、局所的な節点番号と全体的な節点番号の対応付けが必要になるなどの制約がある（例えば、非特許文献１参照。）。
さらに、有限要素法解析の精度を高めるために要素を小さくすると、要素数や節点数が増加し、データ入力とステップの計算に多くの作業と時間を要することになる。とくに、解析対象物が３次元中実体の複雑なモデルである場合には、要素数が数千万以上の規模になることがあり、剛性行列Ｋと質量行列Ｍの計算量が増加する。このため多くの時間がかかる上に誤差が累積し、出力された結果は不正確なものになる。さらに、剛性行列Ｋと質量行列Ｍの計算量が膨大になると、大きな記憶容量と高速に計算処理を行うことができる高性能の計算機が必要になり、一般に使用されている計算機では対応できなくなるという問題が生じる。
このため、従来から、剛性行列Ｋと質量行列Ｍの計算に関心が集まり、解析対象物の要素分割を工夫し、計算機における計算処理を分散化させるなどの改良が進められている。例えば、複数のプロセッサからなり、各プロセッサにメモリを持つようなマルチプロセッサ上での並列処理方式を採用し、各プロセッサが可能な限り、格納する剛性行列の成分、荷重ベクトルの成分を計算するのに必要となる要素剛性行列、応力を計算するように各プロセッサへの割り当てを行うことにより、バッファの容量を削減する有限要素法の処理方式（例えば、特許文献１参照。）が開示されている。また、要素分割の対象となる構造物の亀裂先端部分の作成位置と分割数を指定するだけで自動的に要素の頂点とその中間点の節点の座標データを作成し、前記亀裂先端部の基準要素に属する節点番号の集合を有する要素データを求め、前記基準要素の節点番号をもとにして他の要素の節点番号を求めて、その要素の節点番号の集合を有する要素データを次々に求める有限要素法解析用モデルの作成方法（例えば、特許文献２参照。）、有限要素法による解析対象物体の二次元又は三次元の有限要素データの作成方法において、要素の分割数の少ない方向から優先的に節点番号を付ける工程を備えた有限要素データの作成方法（例えば、特許文献３参照。）が開示されている。さらに、有限要素法の定型的手順により構造物の解析モデルから剛性行列および集中質量行列を生成し、これらの行列で表される一般固有値問題を標準的固有値問題Ａｙ＝λｙに変換し、集中質量行列が負の成分を含む場合、アーノルディの反復法によりＡを射影し行列Ｔ＝｛ｈｉｊ｝を得るが、このときの初期ベクトルを、対応する集中質量行列の成分が正であるような成分ｈに対しては実数、負であるような成分ｈに対しては純虚数となるように設定することにより、行列Ａを実ヘッセンベルグ行列Ｔに変換し、このヘッセンベルグ行列Ｔが与える標準的固有値問題の固有値をＱＲ法により求める方法およびシステム（例えば、特許文献４参照。）が開示されている。
菊地文雄著「有限要素法概説」１９８０年、Ｐ．５８−６７、株式 会社サイエンス社
特開平０５−１０８６９６号公報
特開平１１−２７２６４９号公報
特開平１１−２８３０５０号公報
特開２００１−２５６２１６号公報
しかしながら、特許文献１の開示技術では、解析対象物が３次元の中実体の複雑なモデルになると、要素数が数千万以上の規模になることがあるため、計算処理方法を改良しても、計算機では処理できない膨大な計算量になる可能性がある。また、特許文献２ないし４のいずれの開示技術でも、本発明のように剛性行列の要素Ｋｉｊを決定する式そのものを改良することにより、計算時間を大幅に短縮し、解析精度を向上させるという新たな提案はなされていない。
本発明は上記事情に鑑みてなされたものであり、要素分割のデータ入力に際し、データ入力作業から各要素の節点の局所番号付けする作業をなくし、高次要素を用いた計算であっても、データを与えるとほとんど誤差なく短時間に、剛性行列Ｋと質量行列Ｍを決定する式により計算して求め（以下、「直接計算」と呼ぶ。）、固有関数を精度良く可視化して出力することができる有限要素法直接計算プログラムおよび解析方法を提供することを目的とする。
本発明は上記事情に鑑みてなされたものであり、要素分割のデータ入力に際し、データ入力作業から各要素の節点の局所番号付けする作業をなくし、高次要素を用いた計算であっても、データを与えるとほとんど誤差なく短時間に、剛性行列Ｋと質量行列Ｍを決定する式により計算して求め（以下、「直接計算」と呼ぶ。）、固有関数を精度良く可視化して出力することができる有限要素法直接計算プログラムおよび解析方法を提供することを目的とする。
上記課題を解決するための手段として、本発明は以下の特徴を有している。
本発明のプログラムは、有限要素法により解析対象物の数値モデルを作成し、数値解析を行い、解析結果を可視化するプログラムにおいて、有限要素法の基底関数から決まる剛性行列Ｋと質量行列Ｍを直接計算することにより、与えられた領域のラプラシアンの固有値問題の固有値と固有関数を求めることを特徴とする。
本発明によれば、剛性行列の要素Ｋｉｊを決定する式そのものを改良することにより、要素分割のデータ入力に際し、データ入力作業から各要素の節点の局所番号付けを省くことができ、データを与えるとほとんど誤差なく短時間に剛性行列Ｋと質量行列Ｍを直接計算して求め、固有値と固有関数を出力させることができる。本発明では、各要素の節点のつながりが把握できればよいため、従来必要とされていた各要素の節点番号付けの順序や、局所的な節点番号と全体的な節点番号の対応付けは不要になる。設定できる要素数と節点数の制限をなくしたことにより、１次元の線分、２次元の平面はもちろん、３次元の中空体、又は３次元の中実体の複雑な解析対象物の高次要素を用いた計算にも対応できるようにしている。
本発明のプログラムは、有限要素法により解析対象物の数値モデルを作成し、数値解析を行い、解析結果を可視化するプログラムにおいて、有限要素法の基底関数から決まる剛性行列Ｋと質量行列Ｍを直接計算することにより、与えられた領域のラプラシアンの固有値問題の固有値と固有関数を求めることを特徴とする。
本発明によれば、剛性行列の要素Ｋｉｊを決定する式そのものを改良することにより、要素分割のデータ入力に際し、データ入力作業から各要素の節点の局所番号付けを省くことができ、データを与えるとほとんど誤差なく短時間に剛性行列Ｋと質量行列Ｍを直接計算して求め、固有値と固有関数を出力させることができる。本発明では、各要素の節点のつながりが把握できればよいため、従来必要とされていた各要素の節点番号付けの順序や、局所的な節点番号と全体的な節点番号の対応付けは不要になる。設定できる要素数と節点数の制限をなくしたことにより、１次元の線分、２次元の平面はもちろん、３次元の中空体、又は３次元の中実体の複雑な解析対象物の高次要素を用いた計算にも対応できるようにしている。
本発明のプログラムは、前記解析対象物が１次元モデルである場合に、解析対象物を１次元モデルによる解析方法より解析し、前記解析対象物が２次元モデルである場合に、解析対象物を１次元モデルによる解析方法、２次元モデルによる解析方法又はこれらの組み合わせにより解析し、前記解析対象物が３次元モデルである場合に、前記解析対象物を、１次元モデルによる解析方法、２次元モデルによる解析方法、３次元モデルによる解析方法又はこれらの任意の組み合わせにより解析することを特徴とする。
本発明の解析方法は、有限要素法により解析対象物の数値モデルを作成し、数値解析を行い、解析結果を可視化する解析方法において、有限要素法の基底関数から決まる剛性行列Ｋと質量行列Ｍを直接計算することにより、与えられた領域のラプラシアンの固有値問題の固有値と固有関数を求めることを特徴とする。
本発明の解析方法は、前記解析対象物が１次元モデルである場合に、解析対象物を１次元モデルによる解析方法より解析し、前記解析対象物が２次元モデルである場合に解析対象物を１次元モデルによる解析方法、２次元モデルによる解析方法又はこれらの組み合わせにより解析し、前記解析対象物が３次元モデルである場合に、前記解析対象物を、１次元モデルによる解析方法、２次元モデルによる解析方法、３次元モデルによる解析方法又はこれらの任意の組み合わせにより解析することを特徴とする。
本発明の解析方法は、剛性行列Ｋと質量行列Ｍを直接計算する場合に、前記解析対象物が２次元の平面のモデルであるとき、剛性行列の要素Ｋｉｊについて、
の定義により、剛性行列の要素Ｋｉｊを決定する式が、
ｉ≠ｊのとき、
ｉ＝ｊのとき、
のように得られることを特徴とする。
本発明の剛性行列の要素Ｋｉｊを決定する式は、式そのものを改良したものであり、２次元の平面領域の解析対象物を三角形要素εに分割し、三角形の各線分の長さを分子に、それに対応する三角形の面積を分母にした量で定式化したものである。これらの式は、３次元の中空体においても、同様に適用することができる。
図８は、２次元平面領域および３次元中空体における、ｉ≠ｊのときの剛性行列Ｋｉｊを決定する式を説明するための概念図である。 節点Ｐｉ と節点Ｐｊ が異なるときは、図８のようになる。上記ｉ≠ｊのときのＫｉｊの式の分母は対応する三角形の面積を表し、ａｒｅａ（ｅμ）は三つの節点Ｐｉ ，Ｐｊ ，Ｐkからなる三角形を、ａｒｅａ（ｅｖ） は三つの節点Ｐｉ ，Ｐｊ ， Ｐｌ からなる三角形を表す。分子はそれぞれの辺の長さの２乗などを足したり引いたりしたものである。
図９は、２次元平面領域および３次元中空体における、ｉ­＝ｊのときの剛性行列Ｋｉｉを決定する式を説明するための概念図である。 節点Ｐｉ と節点Ｐｊ が同じときは、図９のようになる。節点Ｐｉ
に接続する節点と三角形が書かれている。上記ｉ­＝ｊのときのＫｉｉの式は、各三角形のＰｉ を含まない辺の長さを２乗し、三角形の面積で割った量の和で与えられる。
従って、局所的な節点番号と全体的な節点番号を対応づけは不要となる。３次元の中空体の場合も同様である。
の定義により、剛性行列の要素Ｋｉｊを決定する式が、
ｉ≠ｊのとき、
ｉ＝ｊのとき、
のように得られることを特徴とする。
本発明の剛性行列の要素Ｋｉｊを決定する式は、式そのものを改良したものであり、２次元の平面領域の解析対象物を三角形要素εに分割し、三角形の各線分の長さを分子に、それに対応する三角形の面積を分母にした量で定式化したものである。これらの式は、３次元の中空体においても、同様に適用することができる。
図８は、２次元平面領域および３次元中空体における、ｉ≠ｊのときの剛性行列Ｋｉｊを決定する式を説明するための概念図である。 節点Ｐｉ と節点Ｐｊ が異なるときは、図８のようになる。上記ｉ≠ｊのときのＫｉｊの式の分母は対応する三角形の面積を表し、ａｒｅａ（ｅμ）は三つの節点Ｐｉ ，Ｐｊ ，Ｐkからなる三角形を、ａｒｅａ（ｅｖ） は三つの節点Ｐｉ ，Ｐｊ ， Ｐｌ からなる三角形を表す。分子はそれぞれの辺の長さの２乗などを足したり引いたりしたものである。
図９は、２次元平面領域および３次元中空体における、ｉ­＝ｊのときの剛性行列Ｋｉｉを決定する式を説明するための概念図である。 節点Ｐｉ と節点Ｐｊ が同じときは、図９のようになる。節点Ｐｉ
に接続する節点と三角形が書かれている。上記ｉ­＝ｊのときのＫｉｉの式は、各三角形のＰｉ を含まない辺の長さを２乗し、三角形の面積で割った量の和で与えられる。
従って、局所的な節点番号と全体的な節点番号を対応づけは不要となる。３次元の中空体の場合も同様である。
本発明の解析方法は、剛性行列Ｋと質量行列Ｍを直接計算する場合に、前記解析対象物が２次元の平面のモデルであるとき、質量行列の要素Ｍｉｊについて、
の定義により、質量行列の要素Ｍｉｊを決定する式が、
ｉ≠ｊのとき、
ｉ＝ｊのとき、
のように得られることを特徴とする。
本発明の質量行列の要素Ｍｉｊを決定する式は、３次元の中空体においても、同様に適用することができる。
の定義により、質量行列の要素Ｍｉｊを決定する式が、
ｉ≠ｊのとき、
ｉ＝ｊのとき、
のように得られることを特徴とする。
本発明の質量行列の要素Ｍｉｊを決定する式は、３次元の中空体においても、同様に適用することができる。
本発明の解析方法は、剛性行列Ｋと質量行列Ｍを直接計算する場合に、前記解析対象物が３次元の中空体のモデルであるとき、剛性行列の要素Ｋｉｊについて、
の定義により、剛性行列の要素Ｋｉｊを決定する式が、
ｉ≠ｊのとき、
ｉ＝ｊのとき、
のように得られることを特徴とする。
の定義により、剛性行列の要素Ｋｉｊを決定する式が、
ｉ≠ｊのとき、
ｉ＝ｊのとき、
のように得られることを特徴とする。
本発明の解析方法は、剛性行列Ｋと質量行列Ｍを直接計算する場合に、前記解析対象物が３次元の中空体のモデルであるとき、質量行列の要素Ｍｉｊについて、
の定義により、質量行列の要素Ｍｉｊを決定する式が、
ｉ≠ｊのとき、
ｉ＝ｊのとき、
のように得られることを特徴とする。
の定義により、質量行列の要素Ｍｉｊを決定する式が、
ｉ≠ｊのとき、
ｉ＝ｊのとき、
のように得られることを特徴とする。
本発明の解析方法は、剛性行列Ｋと質量行列Ｍを直接計算する場合に、前記解析対象物が３次元の中実体のモデルであるとき、剛性行列の要素Ｋｉｊについて、
の定義により、剛性行列の要素Ｋｉｊを決定する式が、
ｉ≠ｊのとき、
ｉ＝ｊのとき、
のように得られることを特徴とする。
本発明の剛性行列の要素Ｋｉｊを決定する式は、式そのものを改良したものであり、３次元の中実体の解析対象物を四面体要素ｅに分割し、四面体の各線分の長さを分子に、それに対応する四面体の体積を分母にし、四面体要素ｅすべてにわたる和として定式化したものである。
の定義により、剛性行列の要素Ｋｉｊを決定する式が、
ｉ≠ｊのとき、
ｉ＝ｊのとき、
のように得られることを特徴とする。
本発明の剛性行列の要素Ｋｉｊを決定する式は、式そのものを改良したものであり、３次元の中実体の解析対象物を四面体要素ｅに分割し、四面体の各線分の長さを分子に、それに対応する四面体の体積を分母にし、四面体要素ｅすべてにわたる和として定式化したものである。
本発明の解析方法は、剛性行列Ｋと質量行列Ｍを直接計算する場合に、前記解析対象物が３次元の中実体のモデルであるとき、質量行列の要素Ｍｉｊについて、
の定義により、質量行列の要素Ｍｉｊを決定する式が、
ｉ≠ｊのとき、
ｉ＝ｊのとき、
のように得られることを特徴とする。
本発明の質量行列の要素Ｍｉｊを決定する式は、四面体要素ｅすべてにわたる和として定式化したものであり、Ｖｏｌ（ｅ）は四面体の体積を表す。
の定義により、質量行列の要素Ｍｉｊを決定する式が、
ｉ≠ｊのとき、
ｉ＝ｊのとき、
のように得られることを特徴とする。
本発明の質量行列の要素Ｍｉｊを決定する式は、四面体要素ｅすべてにわたる和として定式化したものであり、Ｖｏｌ（ｅ）は四面体の体積を表す。
本発明の解析方法は、剛性行列Ｋと質量行列Ｍを直接計算することにより、与えられた領域のラプラシアンの固有値問題の固有値と固有関数を求める場合に、固有値問題を、
で与えたとき、第ｋ番目の固有値λｋ
に対応する固有ベクトルｕｋの第ｉ成分を、第ｋ番目の固有関数φｋ の第ｉ番目の節点の値とし、隣接する節点間の値は線形補間によって定めることにより第ｋ番目の固有関数φｋ
の値を求め、その固有関数φｋ の値を用いて解析対象物を可視化することを特徴とする。
本発明の可視化方法は、剛性行列Ｋと質量行列Ｍを直接計算することにより、固有空間に所属する第ｋ番目の固有関数φｋの値を用いて解析対象物を可視化して短時間に精度良く出力することができる。
図１０は、解析対象物が２次元平面および３次元中空体のモデルであるときの、隣接する節点間の値を線形補間によって定める方法の説明をするための概念図である。
で与えたとき、第ｋ番目の固有値λｋ
に対応する固有ベクトルｕｋの第ｉ成分を、第ｋ番目の固有関数φｋ の第ｉ番目の節点の値とし、隣接する節点間の値は線形補間によって定めることにより第ｋ番目の固有関数φｋ
の値を求め、その固有関数φｋ の値を用いて解析対象物を可視化することを特徴とする。
本発明の可視化方法は、剛性行列Ｋと質量行列Ｍを直接計算することにより、固有空間に所属する第ｋ番目の固有関数φｋの値を用いて解析対象物を可視化して短時間に精度良く出力することができる。
図１０は、解析対象物が２次元平面および３次元中空体のモデルであるときの、隣接する節点間の値を線形補間によって定める方法の説明をするための概念図である。
上記課題を解決するための手段により、本発明によれば、有限要素法により対象物を解析する場合に、要素分割のデータ入力に際し、データ入力作業から各要素の節点の局所番号付けを省くことができた。また、従来必要とされていた各要素の節点番号付けの順序や局所的な節点番号と全体的な節点番号の対応付けを不要にしたことによって、設定できる要素数と節点数の制限をなくすことができた。本発明によれば、１次元の線分、２次元の平面はもちろん、３次元の中空体、又は３次元の中実体の複雑な解析対象物を解析する場合であっても、データを与えると、従来技術に比べて、ほとんど誤差なく短時間に、剛性行列Ｋと質量行列Ｍを直接計算して求めることができた。
また、剛性行列Ｋと質量行列Ｍを直接計算した結果を用いて、ディリクレ境界値問題とノイマン境界値固有問題の固有値と固有関数、および高次の固有値と固有関数を容易かつ精度良く可視化して出力することができることから、定量的で具体的な解析をすることができ、取り扱いが容易で実用的な有限要素法直接計算プログラムおよび解析方法を提供することができた。
また、剛性行列Ｋと質量行列Ｍを直接計算した結果を用いて、ディリクレ境界値問題とノイマン境界値固有問題の固有値と固有関数、および高次の固有値と固有関数を容易かつ精度良く可視化して出力することができることから、定量的で具体的な解析をすることができ、取り扱いが容易で実用的な有限要素法直接計算プログラムおよび解析方法を提供することができた。
本発明を実施するための最良の形態を、図１から図３に基づいて説明する。図１から図３では、連続体である３次元の中実体の解析対象物を、有限の形状と寸法を持つ四面体要素ｅに分割し、各要素の集合により近似して扱う。
ただし、これらは一実施形態にすぎず、本発明の特許請求の範囲を限定するものではない。
ただし、これらは一実施形態にすぎず、本発明の特許請求の範囲を限定するものではない。
図１は、有限要素法直接計算プログラムの解析方法の手順を示したフローチャートである。ステップＳ０１〜Ｓ１４により剛性行列Ｋと質量行列Ｍを決定して出力し、ステップＳ１５〜Ｓ１６により固有値と固有関数を求めて出力する手順により、解析する。
図２は、図４における領域Ω１ の四面体要素の３５個の節点の座標を与える、節点の座標ファイルの例である。節点とは各四面体要素の頂点をいう。図２の節点の座標ファイルは、前もって準備しておく。節点数を１行目に入力し、２行目以降は節点の座標を１行ずつ入力する。この順番が節点番号となるので、節点番号を与える必要はない。座標は、左からｘ座標、ｙ座標およびｚ座標とする。例えば、１行目は節点数が３５であり、２行目は節点番号１のｘ座標が０．５、ｙ座標が０．５、ｚ座標が０．５であることを意味する。
図３は、図４における領域Ω１
の９６個の四面体要素の節点の連結関係を与える、連結ファイルの例である。連結ファイルとは、各要素をなす節点の番号のファイルをいう。節点番号は、３５個の節点Ｐ１
〜Ｐ35 の添字を用いて１〜３５のように表す。図３の連結ファイルは、前もって準備しておく。図３の連結ファイルは、連結関係にある４点の節点番号の組の総数を１行目に入力し、２行目以降は連結関係にある４点の節点番号の組を１行ずつ入力する。この場合、四面体をなす節点番号が分かればよく、時計、反時計まわりのどちらの順番で入力してもよい。図３の２行目は節点番号１６，３５，２３，２からなる四面体の要素を表している。
の９６個の四面体要素の節点の連結関係を与える、連結ファイルの例である。連結ファイルとは、各要素をなす節点の番号のファイルをいう。節点番号は、３５個の節点Ｐ１
〜Ｐ35 の添字を用いて１〜３５のように表す。図３の連結ファイルは、前もって準備しておく。図３の連結ファイルは、連結関係にある４点の節点番号の組の総数を１行目に入力し、２行目以降は連結関係にある４点の節点番号の組を１行ずつ入力する。この場合、四面体をなす節点番号が分かればよく、時計、反時計まわりのどちらの順番で入力してもよい。図３の２行目は節点番号１６，３５，２３，２からなる四面体の要素を表している。
本発明の剛性行列の要素Ｋｉｊと質量行列の要素Ｍｉｊを決定する式の計算方法について述べる。
ｉ≠ｊのとき、剛性行列Ｋと質量行列Ｍの（ｉ，ｊ）成分は、［数７］、［数８］および［数９］により計算される。ここで、ｉ、ｊ、ｋ、ｌは自然数である。
節点Ｐの座標を（ｘ（Ｐ），ｙ（Ｐ），ｚ（Ｐ））で表したとき、始点Ｐ、終点ＱのベクトルＰＱと始点Ｒ、終点ＳのベクトルＲＳの内積は次のように表される。
ここで、右辺の和は辺Ｐｉ
Ｐｊ を含む四面体ｅすべてにわたる和であり、四面体ｅの４つの節点Ｐｉ ，Ｐｊ ，Ｐk （ｅ），Ｐl （ｅ）で表す。また、ＰｉＰl （ｅ）等は、始点をＰｉ 、終点をＰl （ｅ）とするベクトルを表し、（ＰｉＰl （ｅ），Ｐｊ Ｐl （ｅ））は２つのベクトルＰｉＰl （ｅ）とＰｊ Ｐl （ｅ）との、［数７］における内積を表す。Ｖｏｌ（ｅ）は四面体ｅの体積を表している。
ここで、右辺の和は辺Ｐｉ
Ｐｊ を含む四面体ｅすべてにわたる和であり、Ｖｏｌ（ｅ）は四面体ｅの体積を表している。
ｉ＝ｊのとき、剛性行列Ｋと質量行列Ｍの（ｉ，ｉ）成分は、［数７］、［数１０］および［数１１］により計算される。ここで、ｉ、ｊ、ｋ、ｌは自然数である。
ここで、右辺の和は節点Ｐｉ
を含む四面体ｅすべてにわたる和である。Ｖｏｌ（ｅ）は四面体ｅの体積を表す。
ここで、右辺の和は節点Ｐｉ
を含む四面体ｅすべてにわたる和である。Ｖｏｌ（ｅ）は四面体ｅの体積を表す。
以下、図１に示すステップＳ０１〜Ｓ１４により剛性行列Ｋと質量行列Ｍを決定して出力する解析方法について説明する。
ステップＳ０１では、節点の座標ファイルを読み込む。節点の座標ファイルは、図２のように前もって準備しておく。
ステップＳ０２では、連結ファイルを読み込む。連結ファイルは、図３のように前もって準備しておく。
ステップＳ０３では、Ｌ×Ｌ行列のすべての成分を零とする初期化を行っている。Ｌ×Ｌ行列は、解決すべき問題に対し、必要な行列とする。ただし、Ｌは１≦Ｌ≦Ｎの定数とし、Ｌはプログラムの中で与える。ここで、Ｎは自然数であり、節点数を表す。
ステップＳ０１では、節点の座標ファイルを読み込む。節点の座標ファイルは、図２のように前もって準備しておく。
ステップＳ０２では、連結ファイルを読み込む。連結ファイルは、図３のように前もって準備しておく。
ステップＳ０３では、Ｌ×Ｌ行列のすべての成分を零とする初期化を行っている。Ｌ×Ｌ行列は、解決すべき問題に対し、必要な行列とする。ただし、Ｌは１≦Ｌ≦Ｎの定数とし、Ｌはプログラムの中で与える。ここで、Ｎは自然数であり、節点数を表す。
ステップＳ０４〜Ｓ０７では、ｉ≠ｊのときの、剛性行列の要素Ｋｉｊと質量行列の要素Ｍｉｊを決定する式の計算をしている。
ステップＳ０４では、ステップＳ０２で読み込んだ節点Ｐｉ
およびＰｊ（ｉ＜ｊ）をもつ四面体ｅの残りの節点がＰｋ ，Ｐｌ ならば、ステップＳ０１で読み込んだ各節点の座標（ｘ（ｉ），ｙ（ｉ），ｚ（ｉ））、（ｘ（ｊ），ｙ（ｊ），ｚ（ｊ））、（ｘ（ｋ），ｙ（ｋ），ｚ（ｋ））、（ｘ（ｌ）,ｙ（ｌ）,ｚ（ｌ））を用いて、内積（ＰｉＰｌ，ＰｊＰｌ）、（ＰｉＰｋ，ＰｊＰｋ）、（ＰｉＰｌ，ＰｊＰｋ）、（ＰｉＰｋ，ＰｊＰｌ）を計算する。
ステップＳ０５では、四面体ｅの体積Ｖｏｌ（ｅ）を求める。
ステップＳ０６では、ＡＫ（ｉ，ｊ）、ＡＭ（ｉ，ｊ）を求める。
［数８］より
の値をＡＫ（ｉ，ｊ）、
［数７］より
の値をＡＭ（ｉ，ｊ）とする。
ステップＳ０７では、ステップＳ０６と同様に、四面体ｅの辺Ｐｉ
Ｐｊ 以外についても計算する。具体的には、
ＡＫ（ｉ，ｋ）、ＡＫ（ｉ，ｌ）、ＡＫ（ｊ，ｋ）、ＡＫ（ｊ，ｌ）、ＡＫ（ｋ，ｌ）、
ＡＭ（ｉ，ｋ）、ＡＭ（ｉ，ｌ）、ＡＭ（ｊ，ｋ）、ＡＭ（ｊ，ｌ）、ＡＭ（ｋ，ｌ）
である。
ステップＳ０４では、ステップＳ０２で読み込んだ節点Ｐｉ
およびＰｊ（ｉ＜ｊ）をもつ四面体ｅの残りの節点がＰｋ ，Ｐｌ ならば、ステップＳ０１で読み込んだ各節点の座標（ｘ（ｉ），ｙ（ｉ），ｚ（ｉ））、（ｘ（ｊ），ｙ（ｊ），ｚ（ｊ））、（ｘ（ｋ），ｙ（ｋ），ｚ（ｋ））、（ｘ（ｌ）,ｙ（ｌ）,ｚ（ｌ））を用いて、内積（ＰｉＰｌ，ＰｊＰｌ）、（ＰｉＰｋ，ＰｊＰｋ）、（ＰｉＰｌ，ＰｊＰｋ）、（ＰｉＰｋ，ＰｊＰｌ）を計算する。
ステップＳ０５では、四面体ｅの体積Ｖｏｌ（ｅ）を求める。
ステップＳ０６では、ＡＫ（ｉ，ｊ）、ＡＭ（ｉ，ｊ）を求める。
［数８］より
の値をＡＫ（ｉ，ｊ）、
［数７］より
の値をＡＭ（ｉ，ｊ）とする。
ステップＳ０７では、ステップＳ０６と同様に、四面体ｅの辺Ｐｉ
Ｐｊ 以外についても計算する。具体的には、
ＡＫ（ｉ，ｋ）、ＡＫ（ｉ，ｌ）、ＡＫ（ｊ，ｋ）、ＡＫ（ｊ，ｌ）、ＡＫ（ｋ，ｌ）、
ＡＭ（ｉ，ｋ）、ＡＭ（ｉ，ｌ）、ＡＭ（ｊ，ｋ）、ＡＭ（ｊ，ｌ）、ＡＭ（ｋ，ｌ）
である。
ステップＳ０８〜Ｓ１１では、ｉ＝ｊのときの、剛性行列の要素Ｋｉｊと質量行列の要素Ｍｉｊを決定する式の計算をしている。
ステップＳ０８では、節点Ｐｉ
を含む四面体ｅの残りの節点がＰｊ ，Ｐk ，Ｐl ならばステップＳ０１で読み込んだ各節点の座標（ｘ（ｉ），ｙ（ｉ），ｚ（ｉ））、（ｘ（ｊ），ｙ（ｊ），ｚ（ｊ））、（ｘ（ｋ），ｙ（ｋ），ｚ（ｋ））（ｘ（ｌ），ｙ（ｌ），ｚ（ｌ））を用いて、内積（Ｐｊ Ｐｌ，Ｐｊ
Ｐｋ ）と‖ＰｊＰｌ‖２、‖ＰｊＰｋ‖２
を計算する。
ステップＳ０９では、四面体ｅの体積Ｖｏｌ（ｅ）を求める。
ステップＳ１０では、ＡＫ（ｉ，ｉ）、ＡＭ（ｉ，ｉ）を求める。
［数８］より
の値をＡＫ（ｉ，ｉ）、
［数９］より
の値をＡＭ（ｉ，ｉ）とする。
ステップＳ１１では、ステップＳ１０と同様に、四面体ｅの点Ｐｉ
以外についても計算する。具体的には、
ＡＫ（ｊ，ｊ）、ＡＫ（ｋ，ｋ）、ＡＫ（ｌ，ｌ）、ＡＭ（ｊ，ｊ）、ＡＭ（ｋ，ｋ）、
ＡＭ（ｌ，ｌ）
である。
ステップＳ０８では、節点Ｐｉ
を含む四面体ｅの残りの節点がＰｊ ，Ｐk ，Ｐl ならばステップＳ０１で読み込んだ各節点の座標（ｘ（ｉ），ｙ（ｉ），ｚ（ｉ））、（ｘ（ｊ），ｙ（ｊ），ｚ（ｊ））、（ｘ（ｋ），ｙ（ｋ），ｚ（ｋ））（ｘ（ｌ），ｙ（ｌ），ｚ（ｌ））を用いて、内積（Ｐｊ Ｐｌ，Ｐｊ
Ｐｋ ）と‖ＰｊＰｌ‖２、‖ＰｊＰｋ‖２
を計算する。
ステップＳ０９では、四面体ｅの体積Ｖｏｌ（ｅ）を求める。
ステップＳ１０では、ＡＫ（ｉ，ｉ）、ＡＭ（ｉ，ｉ）を求める。
［数８］より
の値をＡＫ（ｉ，ｉ）、
［数９］より
の値をＡＭ（ｉ，ｉ）とする。
ステップＳ１１では、ステップＳ１０と同様に、四面体ｅの点Ｐｉ
以外についても計算する。具体的には、
ＡＫ（ｊ，ｊ）、ＡＫ（ｋ，ｋ）、ＡＫ（ｌ，ｌ）、ＡＭ（ｊ，ｊ）、ＡＭ（ｋ，ｋ）、
ＡＭ（ｌ，ｌ）
である。
ステップＳ１２では、すべての四面体要素に対してＳ０４〜Ｓ１１で行ったことを繰り返す。そして、各要素を加えた和集合が、行列ＫおよびＭの要素となる。
ステップＳ１３では、ステップＳ１２で求めた行列ＫおよびＭの成分を用いて、（ｊ，ｉ）成分（ｉ＞ｊ）を求める。連結関係のない成分はすべて零であり、ステップＳ０３で初期化した結果をそのまま用いる。以上で、行列ＫおよびＭの成分が決定できる。
ステップＳ１４では、行列ＫおよびＭのすべての成分を出力する。このステップは省略してもよい。
ステップＳ１３では、ステップＳ１２で求めた行列ＫおよびＭの成分を用いて、（ｊ，ｉ）成分（ｉ＞ｊ）を求める。連結関係のない成分はすべて零であり、ステップＳ０３で初期化した結果をそのまま用いる。以上で、行列ＫおよびＭの成分が決定できる。
ステップＳ１４では、行列ＫおよびＭのすべての成分を出力する。このステップは省略してもよい。
次に、ステップＳ１５〜Ｓ１６により固有値と固有関数を求めて出力することにより、解析対象物を可視化する解析方法について説明する。
ステップＳ１５では、［数１２］で表される行列の一般固有値問題を解く。ここでは既に使用されている行列の一般固有値問題を解くプログラムを組み込んで解いているが、ステップＳ１４で出力してある行列ＫおよびＭの成分を用いて別のプログラムを使用しても構わない。
ステップＳ１５では、［数１２］で表される行列の一般固有値問題を解く。ここでは既に使用されている行列の一般固有値問題を解くプログラムを組み込んで解いているが、ステップＳ１４で出力してある行列ＫおよびＭの成分を用いて別のプログラムを使用しても構わない。
ただし、Δはラプラシアンであって、φのラプラシアンは次のように表される。
［数１２］における第ｋ番目の固有値を、［数１３］における第ｋ番目の固有値λｋと定める。また、［数１２］の第ｋ番目の固有値λｋ
に対応する固有ベクトルｕｋから、［数１３］における第ｋ番目の固有関数φｋを得るには、固有ベクトルｕｋの第ｉ成分を固有関数φｋ
の第ｉ番目の節点の値とし、隣接する節点間の値を線形補間によって固有関数φｋ の値を定める。
すなわち、（ｘ，ｙ，ｚ）∈四面体Ｐｉ Ｐｊ Ｐk Ｐｌ において、
φｋ
（ｘ，ｙ，ｚ）＝ａ×（ｕの第ｉ成分）＋ｂ×（ｕの第ｊ成分）＋ｃ×（ｕの第ｋ成分）＋ｄ×（ｕの第ｌ成分）
によって与える。
ただし、０≦ａ，ｂ，ｃ，ｄ≦１、ａ＋ｂ＋ｃ＋ｄ≦１である。
に対応する固有ベクトルｕｋから、［数１３］における第ｋ番目の固有関数φｋを得るには、固有ベクトルｕｋの第ｉ成分を固有関数φｋ
の第ｉ番目の節点の値とし、隣接する節点間の値を線形補間によって固有関数φｋ の値を定める。
すなわち、（ｘ，ｙ，ｚ）∈四面体Ｐｉ Ｐｊ Ｐk Ｐｌ において、
φｋ
（ｘ，ｙ，ｚ）＝ａ×（ｕの第ｉ成分）＋ｂ×（ｕの第ｊ成分）＋ｃ×（ｕの第ｋ成分）＋ｄ×（ｕの第ｌ成分）
によって与える。
ただし、０≦ａ，ｂ，ｃ，ｄ≦１、ａ＋ｂ＋ｃ＋ｄ≦１である。
ステップＳ１６では、［数１３］における第ｋ番目の固有値λｋ
、固有関数φｋ の出力を行う。ただし、ｋは１≦ｋ≦Ｌの定数とする。
以上で、図１に示すステップＳ０１〜Ｓ１９の一連の処理を終了する。
、固有関数φｋ の出力を行う。ただし、ｋは１≦ｋ≦Ｌの定数とする。
以上で、図１に示すステップＳ０１〜Ｓ１９の一連の処理を終了する。
次に、本発明の剛性行列Ｋの計算時間について、従来の計算方法と比較するために、２次元の場合について説明する。
従来の計算方法については、非特許文献１に記載の剛性行列の計算法に従って演算回数を調べる。三角形要素の１つをＰ１ （ｘ１ ，ｙ１ ），Ｐ２
（ｘ２ ，ｙ２ ），Ｐ３ （ｘ３ ，ｙ３ ）とし、反時計回りに並んでいるものとする。従来の方法では、３点がＰ１ ，Ｐ２ ，Ｐ３ の順に反時計回りに並んでいる必要があるが、本発明の方法では向きを考慮する必要はない。
従来の計算方法については、非特許文献１に記載の剛性行列の計算法に従って演算回数を調べる。三角形要素の１つをＰ１ （ｘ１ ，ｙ１ ），Ｐ２
（ｘ２ ，ｙ２ ），Ｐ３ （ｘ３ ，ｙ３ ）とし、反時計回りに並んでいるものとする。従来の方法では、３点がＰ１ ，Ｐ２ ，Ｐ３ の順に反時計回りに並んでいる必要があるが、本発明の方法では向きを考慮する必要はない。
従来の計算方法によると、
（１） 面積Ｓの２倍に等しい行列式Ｄの計算
Ｄ＝２×Ｓ＝ｘ１×（ｙ２
−ｙ３ ）＋ｘ２×（ｙ３ −ｙ１ ）＋ｘ３×（ｙ１
−ｙ２ ）
と計算できるので、加算２回、減算３回、乗算３回の合計８回の演算が必要となる。
（２） 面積座標での近似関数の係数の計算
は減算１回、除算１回でのそれぞれで行うので、合計１２回の演算が必要となる。
（３）Ｓ＝Ｄ／２より、除算１回が必要となる。
（４）要素係数行列の計算
（２）で求めた係数を用いて、Ａｉｊ＝Ｓ×（ｂｊ×ｂｉ＋ｃｊ×ｃｉ）を計算する。
すなわち、Ａ11 ，Ａ22，Ａ33 ，Ａ12（＝Ａ21），Ａ23（＝Ａ32），Ａ31（＝Ａ13）のそれぞれで乗算３回、加算１回を６回行うので、合計２４回の演算が必要となる。
（５）９つの要素を加える
合計９回の演算が必要となる。
以上より、要素がＮならば、（８＋１２＋１＋２４＋９）Ｎ＝５４Ｎの演算回数が必要となる。
（１） 面積Ｓの２倍に等しい行列式Ｄの計算
Ｄ＝２×Ｓ＝ｘ１×（ｙ２
−ｙ３ ）＋ｘ２×（ｙ３ −ｙ１ ）＋ｘ３×（ｙ１
−ｙ２ ）
と計算できるので、加算２回、減算３回、乗算３回の合計８回の演算が必要となる。
（２） 面積座標での近似関数の係数の計算
は減算１回、除算１回でのそれぞれで行うので、合計１２回の演算が必要となる。
（３）Ｓ＝Ｄ／２より、除算１回が必要となる。
（４）要素係数行列の計算
（２）で求めた係数を用いて、Ａｉｊ＝Ｓ×（ｂｊ×ｂｉ＋ｃｊ×ｃｉ）を計算する。
すなわち、Ａ11 ，Ａ22，Ａ33 ，Ａ12（＝Ａ21），Ａ23（＝Ａ32），Ａ31（＝Ａ13）のそれぞれで乗算３回、加算１回を６回行うので、合計２４回の演算が必要となる。
（５）９つの要素を加える
合計９回の演算が必要となる。
以上より、要素がＮならば、（８＋１２＋１＋２４＋９）Ｎ＝５４Ｎの演算回数が必要となる。
本発明の計算方法によると、
（１）面積Ｓの２倍に等しい行列式Ｄの計算
従来の方法と同様に８回の演算が必要となる。
（２）線分の２乗の計算
ＰｉＰj ２ ＝（ｘｉ
−ｘｊ ）２＋（ｙｉ −ｙｊ）２
と計算できるので、加算１回、減算２回、乗算２回となり、
ＰｊＰｋ ２，ＰｋＰｉ ２のそれぞれで行うので、合計１５回の演算が必要となる。
（３）Ｓ１ ＝−０．２５／Ｄ、Ｓ２ ＝０．５／Ｄより、除算２回が必要となる。
（４）要素係数行列の計算
Ａｉｊ＝Ｓ１×（ＰｉＰk ２＋ＰｊＰk ２−ＰｉＰj ２）より、加算１回、減算１回、乗算１回
Ａｉｉ＝Ｓ２×ＰｊＰk ２より、乗算１回を用いて、
Ａ11 ，Ａ22，Ａ33 ，Ａ12（＝Ａ21），Ａ23（＝Ａ32），Ａ31（＝Ａ13）のそれぞれを計算するので、合計１２回の演算が必要となる。
（５）９つの要素を加える
合計９回の演算が必要となる。
以上より、要素がＮならば、（８＋１５＋２＋１２＋９）Ｎ＝４６Ｎの演算回数が必要となる。
（１）面積Ｓの２倍に等しい行列式Ｄの計算
従来の方法と同様に８回の演算が必要となる。
（２）線分の２乗の計算
ＰｉＰj ２ ＝（ｘｉ
−ｘｊ ）２＋（ｙｉ −ｙｊ）２
と計算できるので、加算１回、減算２回、乗算２回となり、
ＰｊＰｋ ２，ＰｋＰｉ ２のそれぞれで行うので、合計１５回の演算が必要となる。
（３）Ｓ１ ＝−０．２５／Ｄ、Ｓ２ ＝０．５／Ｄより、除算２回が必要となる。
（４）要素係数行列の計算
Ａｉｊ＝Ｓ１×（ＰｉＰk ２＋ＰｊＰk ２−ＰｉＰj ２）より、加算１回、減算１回、乗算１回
Ａｉｉ＝Ｓ２×ＰｊＰk ２より、乗算１回を用いて、
Ａ11 ，Ａ22，Ａ33 ，Ａ12（＝Ａ21），Ａ23（＝Ａ32），Ａ31（＝Ａ13）のそれぞれを計算するので、合計１２回の演算が必要となる。
（５）９つの要素を加える
合計９回の演算が必要となる。
以上より、要素がＮならば、（８＋１５＋２＋１２＋９）Ｎ＝４６Ｎの演算回数が必要となる。
従って、従来の計算方法と本発明の計算方法とを比較すると、
（１−４６Ｎ／５４Ｎ）×１００＝１４．８％
となり、およそ１５％の計算時間の短縮が見込まれる。
１５％の計算時間の短縮により、従来の計算方法では３ギガバイトの記憶容量が必要であったものが、本発明の計算方法では２．５ギガバイトの記憶容量で足りることになる。経費にすると、およそ５、６万円の削減が可能になる。
また、要素数５０００（節点２６０１）の場合、従来の方法では０．００１３８秒という結果が得られている。このことから、要素の数が数千万以上の規模となった場合、本発明の計算方法を適用すれば、計算時間の短縮に寄与することが可能になる。
（１−４６Ｎ／５４Ｎ）×１００＝１４．８％
となり、およそ１５％の計算時間の短縮が見込まれる。
１５％の計算時間の短縮により、従来の計算方法では３ギガバイトの記憶容量が必要であったものが、本発明の計算方法では２．５ギガバイトの記憶容量で足りることになる。経費にすると、およそ５、６万円の削減が可能になる。
また、要素数５０００（節点２６０１）の場合、従来の方法では０．００１３８秒という結果が得られている。このことから、要素の数が数千万以上の規模となった場合、本発明の計算方法を適用すれば、計算時間の短縮に寄与することが可能になる。
また、３次元の場合、中空体、中実体のいずれも、剛性行列Ｋの計算時間について厳密に計算ステップを記載してある従来の計算方法は見られないようであるが、３次元の中実体の複雑な解析対象物になると、従来の計算方法では計算に１日ほど要することが見込まれる。
しかしながら、本発明の計算方法を適用すれば、３次元の中空体、又は３次元の中実体の複雑な解析対象物を解析する場合であっても、データを与えると、従来の計算方法に比べて、ほとんど誤差なく短時間に、剛性行列Ｋと質量行列Ｍを直接計算することにより解析することができる。
例えば、節点２３７９、要素数８８７１の直方体の場合には、
（１）
行列ＫおよびＭの決定に、８４秒を要し、
（２）
固有値、固有ベクトルの決定に、４４０秒を要した。
以上より、８４秒＋４４０秒＝５２４秒＝８分４４秒の計算時間で処理できる。
ただし、（１）（２）いずれもデータの入出力等の時間が含まれており、計算のみに比べれば時間は余分にかかっている。
しかしながら、本発明の計算方法を適用すれば、３次元の中空体、又は３次元の中実体の複雑な解析対象物を解析する場合であっても、データを与えると、従来の計算方法に比べて、ほとんど誤差なく短時間に、剛性行列Ｋと質量行列Ｍを直接計算することにより解析することができる。
例えば、節点２３７９、要素数８８７１の直方体の場合には、
（１）
行列ＫおよびＭの決定に、８４秒を要し、
（２）
固有値、固有ベクトルの決定に、４４０秒を要した。
以上より、８４秒＋４４０秒＝５２４秒＝８分４４秒の計算時間で処理できる。
ただし、（１）（２）いずれもデータの入出力等の時間が含まれており、計算のみに比べれば時間は余分にかかっている。
［実施例１］は、３次元の中実体の例である。
図４は、領域Ω１
における立方体領域を四面体要素に分割したときの節点および要素の状態を示した図である。
領域Ω１
＝｛０＜ｘ＜１，０＜ｙ＜１，０＜ｚ＜１｝における
固有値問題：Δｕ＝λｕ ｕ＝ｕ（ｘ，ｙ，ｚ）；（ｘ，ｙ，ｚ）∈Ω１ 、
ディリクレ境界条件：ｕ（ｘ，ｙ，ｚ）＝０ （ｘ，ｙ，ｚ）∈∂Ω１
を図４の四面体要素に分割することにより求める。この四面体分割の節点の座標と連結情報は図２および図３で与えるものとする。
図５は、本発明の解析方法を図４に適用した場合の行列ＫおよびＭの各成分の一覧である。
図６は、本発明の解析方法を図４に適用した場合の固有値および固有ベクトルの数値解である。
図７は、本発明の解析方法を図４に適用した場合の第１〜第５固有関数を可視化したものである。右図は左図の内部を取り出したものである。
図４は、領域Ω１
における立方体領域を四面体要素に分割したときの節点および要素の状態を示した図である。
領域Ω１
＝｛０＜ｘ＜１，０＜ｙ＜１，０＜ｚ＜１｝における
固有値問題：Δｕ＝λｕ ｕ＝ｕ（ｘ，ｙ，ｚ）；（ｘ，ｙ，ｚ）∈Ω１ 、
ディリクレ境界条件：ｕ（ｘ，ｙ，ｚ）＝０ （ｘ，ｙ，ｚ）∈∂Ω１
を図４の四面体要素に分割することにより求める。この四面体分割の節点の座標と連結情報は図２および図３で与えるものとする。
図５は、本発明の解析方法を図４に適用した場合の行列ＫおよびＭの各成分の一覧である。
図６は、本発明の解析方法を図４に適用した場合の固有値および固有ベクトルの数値解である。
図７は、本発明の解析方法を図４に適用した場合の第１〜第５固有関数を可視化したものである。右図は左図の内部を取り出したものである。
［実施例２］は、２次元平面領域の例である。
図１１は、領域Ω１の節点および要素の状態を示した図である。
領域Ω１
＝｛０＜ｘ＜１，０＜ｙ＜１｝における
固有値問題：Δｕ＝λｕ ｕ＝ｕ（ｘ，ｙ）；（ｘ，ｙ）∈Ω１ 、
ディリクレ境界条件：ｕ（ｘ，ｙ）＝０ （ｘ，ｙ）∈∂Ω１
を図１１の三角形要素に分割することにより求める。
図１１は、領域Ω１の節点および要素の状態を示した図である。
領域Ω１
＝｛０＜ｘ＜１，０＜ｙ＜１｝における
固有値問題：Δｕ＝λｕ ｕ＝ｕ（ｘ，ｙ）；（ｘ，ｙ）∈Ω１ 、
ディリクレ境界条件：ｕ（ｘ，ｙ）＝０ （ｘ，ｙ）∈∂Ω１
を図１１の三角形要素に分割することにより求める。
［実施例３］は、２次元平面領域の例である。
領域Ω２ ついて三角形分割を行い、Ω２ の固有値問題およびディリクレ境界条件及びノイマン境界条件に対する固有値問題に対する第１９の固有関数を計算し、可視化した結果を図１５と１６に示す。
領域Ω２ ついて三角形分割を行い、Ω２ の固有値問題およびディリクレ境界条件及びノイマン境界条件に対する固有値問題に対する第１９の固有関数を計算し、可視化した結果を図１５と１６に示す。
［実施例４］は、３次元の中空体の例である。
３次元中空体の曲面における固有関数の結果を図１７〜１９に示す。
固有値については、Ωを上面Ω＋ と下面Ω− の２重に三角形分割して得られる、潰れた曲面とみなし、固有値を出力する。固有関数については、出力された固有関数の上面Ω＋と下面Ω−
の節点の値の差を出力する。
３次元中空体の曲面における固有関数の結果を図１７〜１９に示す。
固有値については、Ωを上面Ω＋ と下面Ω− の２重に三角形分割して得られる、潰れた曲面とみなし、固有値を出力する。固有関数については、出力された固有関数の上面Ω＋と下面Ω−
の節点の値の差を出力する。
熱方程式に対して、温度分布を解析した実施例を示す。１ｍ四方の領域の頂点を、それぞれ（０，０）、（１，０）、（１，１）、（０，１）とおく。初期条件として、ｕ（ｔ，ｘ，ｙ）＝ｕ（０，０．５，０．２５）＝１０４を与える。この領域の時刻ｔ＝０．０２秒のディリクレ境界条件とノイマン境界条件の温度分布を図２０〜２３にそれぞれ示す。
波動方程式に対して、波の動きを解析した実施例を示す。１ｍ四方の領域の頂点をそれぞれ（０，０）、（１，０）、（１，１）、（０，１）とおく。初期条件として、ｕ（ｔ，ｘ，ｙ）＝ｕ（０，０．０１，０．０１）＝ｕ（０，０．９９，０．９９）＝１０４ を与える。
図２４は、この領域の時刻ｔ＝０．０４秒のディリクレ境界条件の時の波の動きを表したものである。
図２５は、この領域の時刻ｔ＝０．０３秒のノイマン境界条件の時の波の動きを表したものである。
図２４は、この領域の時刻ｔ＝０．０４秒のディリクレ境界条件の時の波の動きを表したものである。
図２５は、この領域の時刻ｔ＝０．０３秒のノイマン境界条件の時の波の動きを表したものである。
Claims (11)
- 有限要素法により解析対象物の数値モデルを作成し、数値解析を行い、解析結果を可視化するプログラムにおいて、
前記プログラムは有限要素法の基底関数から決まる剛性行列Ｋと質量行列Ｍを直接計算することにより、与えられた領域のラプラシアンの固有値問題の固有値と固有関数を求める
ことを特徴とするプログラム。 - 前記解析対象物が１次元モデルである場合に、解析対象物を１次元モデルによる解析方法より解析し、
前記解析対象物が２次元モデルである場合に、解析対象物を１次元モデルによる解析方法、２次元モデルによる解析方法又はこれらの組み合わせにより解析し、
前記解析対象物が３次元モデルである場合に、前記解析対象物を、１次元モデルによる解析方法、２次元モデルによる解析方法、３次元モデルによる解析方法又はこれらの任意の組み合わせにより解析する
ことを特徴とする請求項１に記載のプログラム。 - 有限要素法により解析対象物の数値モデルを作成し、数値解析を行い、解析結果を可視化する解析方法において、
前記解析方法は有限要素法の基底関数から決まる剛性行列Ｋと質量行列Ｍを直接計算することにより、与えられた領域のラプラシアンの固有値問題の固有値と固有関数を求める
ことを特徴とする解析方法。 - 前記解析対象物が１次元モデルである場合に、解析対象物を１次元モデルによる解析方法より解析し、
前記解析対象物が２次元モデルである場合に、解析対象物を１次元モデルによる解析方法、２次元モデルによる解析方法又はこれらの組み合わせにより解析し、
前記解析対象物が３次元モデルである場合に、前記解析対象物を、１次元モデルによる解析方法、２次元モデルによる解析方法、３次元モデルによる解析方法又はこれらの任意の組み合わせにより解析する
ことを特徴とする請求項３に記載の解析方法。 - 前記解析方法は、剛性行列Ｋと質量行列Ｍを直接計算する場合に、
前記解析対象物が２次元の平面のモデルであるとき、
剛性行列の要素Ｋｉｊについて、
の定義により、
剛性行列の要素Ｋｉｊを決定する式が、
ｉ≠ｊのとき、
ｉ＝ｊのとき、
のように得られる
ことを特徴とする請求項３又は４に記載の解析方法。 - 前記解析方法は、剛性行列Ｋと質量行列Ｍを直接計算する場合に、
前記解析対象物が２次元の平面のモデルであるとき、
質量行列の要素Ｍｉｊについて、
の定義により、
質量行列の要素Ｍｉｊを決定する式が、
ｉ≠ｊのとき、
ｉ＝ｊのとき、
のように得られる
ことを特徴とする請求項３又は４に記載の解析方法。 - 前記解析方法は、剛性行列Ｋと質量行列Ｍを直接計算する場合に、
前記解析対象物が３次元の中空体のモデルであるとき、
剛性行列の要素Ｋｉｊについて、
の定義により、
剛性行列の要素Ｋｉｊを決定する式が、
ｉ≠ｊのとき、
ｉ＝ｊのとき、
のように得られる
ことを特徴とする請求項３又は４に記載の解析方法。 - 前記解析方法は、剛性行列Ｋと質量行列Ｍを直接計算する場合に、
前記解析対象物が３次元の中空体のモデルであるとき、
質量行列の要素Ｍｉｊについて、
の定義により、
質量行列の要素Ｍｉｊを決定する式が、
ｉ≠ｊのとき、
ｉ＝ｊのとき、
のように得られる
ことを特徴とする請求項３又は４に記載の解析方法。 - 前記解析方法は、剛性行列Ｋと質量行列Ｍを直接計算する場合に、
前記解析対象物が３次元の中実体のモデルであるとき、
剛性行列の要素Ｋｉｊについて、
の定義により、
剛性行列の要素Ｋｉｊを決定する式が、
ｉ≠ｊのとき、
ｉ＝ｊのとき、
のように得られる
ことを特徴とする請求項３又は４に記載の解析方法。 - 前記解析方法は、剛性行列Ｋと質量行列Ｍを直接計算する場合に、
前記解析対象物が３次元の中実体のモデルであるとき、
質量行列の要素Ｍｉｊについて、
の定義により、
質量行列の要素Ｍｉｊを決定する式が、
ｉ≠ｊのとき、
ｉ＝ｊのとき、
のように得られる
ことを特徴とする請求項３又は４に記載の解析方法。 - 前記解析方法は、剛性行列Ｋと質量行列Ｍを直接計算することにより、与えられた領域のラプラシアンの固有値問題の固有値と固有関数を求める場合に、
固有値問題を、
で与えたとき、
第ｋ番目の固有値λｋ
に対応する固有ベクトルｕｋの第ｉ成分を、第ｋ番目の固有関数φｋ の第ｉ番目の節点の値とし、隣接する節点間の値は線形補間によって定めることにより第ｋ番目の固有関数φｋ
の値を求め、
その固有関数φｋ
の値を用いて解析対象物を可視化する
ことを特徴とする請求項３ないし１０のいずれかに記載の解析方法。
Priority Applications (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| JP2005134797A JP4830094B2 (ja) | 2005-05-06 | 2005-05-06 | 計算機 |
Applications Claiming Priority (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| JP2005134797A JP4830094B2 (ja) | 2005-05-06 | 2005-05-06 | 計算機 |
Publications (3)
| Publication Number | Publication Date |
|---|---|
| JP2006313400A true JP2006313400A (ja) | 2006-11-16 |
| JP2006313400A5 JP2006313400A5 (ja) | 2007-02-22 |
| JP4830094B2 JP4830094B2 (ja) | 2011-12-07 |
Family
ID=37534879
Family Applications (1)
| Application Number | Title | Priority Date | Filing Date |
|---|---|---|---|
| JP2005134797A Active JP4830094B2 (ja) | 2005-05-06 | 2005-05-06 | 計算機 |
Country Status (1)
| Country | Link |
|---|---|
| JP (1) | JP4830094B2 (ja) |
Cited By (3)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN103106182A (zh) * | 2011-11-14 | 2013-05-15 | 达索系统西姆利亚公司 | 使用自动多级子结构化生成子结构 |
| CN112364586A (zh) * | 2020-11-11 | 2021-02-12 | 南方电网科学研究院有限责任公司 | 直流换流器的直流回路阻抗计算方法、装置和计算机设备 |
| CN115345047A (zh) * | 2022-08-11 | 2022-11-15 | 南京理工大学 | 一种考虑宽箱梁剪力滞效应的空间梁单元构件方法 |
-
2005
- 2005-05-06 JP JP2005134797A patent/JP4830094B2/ja active Active
Cited By (5)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN103106182A (zh) * | 2011-11-14 | 2013-05-15 | 达索系统西姆利亚公司 | 使用自动多级子结构化生成子结构 |
| US10255392B2 (en) | 2011-11-14 | 2019-04-09 | Dassault Systemes Simulia Corp. | Substructure generation using automated multilevel substructuring |
| CN112364586A (zh) * | 2020-11-11 | 2021-02-12 | 南方电网科学研究院有限责任公司 | 直流换流器的直流回路阻抗计算方法、装置和计算机设备 |
| CN112364586B (zh) * | 2020-11-11 | 2024-05-14 | 南方电网科学研究院有限责任公司 | 直流换流器的直流回路阻抗计算方法、装置和计算机设备 |
| CN115345047A (zh) * | 2022-08-11 | 2022-11-15 | 南京理工大学 | 一种考虑宽箱梁剪力滞效应的空间梁单元构件方法 |
Also Published As
| Publication number | Publication date |
|---|---|
| JP4830094B2 (ja) | 2011-12-07 |
Similar Documents
| Publication | Publication Date | Title |
|---|---|---|
| Pilgerstorfer et al. | Bounding the influence of domain parameterization and knot spacing on numerical stability in isogeometric analysis | |
| Freytag et al. | Finite element analysis in situ | |
| Schaback | A computational tool for comparing all linear PDE solvers: Error-optimal methods are meshless | |
| Cervone et al. | An optimal constrained approach for divergence-free velocity interpolation and multilevel VOF method | |
| Liu et al. | Parameterization for molecular Gaussian surface and a comparison study of surface mesh generation | |
| JP4830094B2 (ja) | 計算機 | |
| US20150100289A1 (en) | Method and system for shapewise comparison | |
| US11003816B2 (en) | Structure analysis device and structure analysis method | |
| US11295050B2 (en) | Structural analysis method and structural analysis apparatus | |
| JP2006313400A5 (ja) | ||
| Khimich et al. | Numerical study of the stability of composite materials on computers of hybrid architecture | |
| Jahanshahloo et al. | Reconstruction of 3D shapes with B-spline surface using diagonal approximation BFGS methods | |
| Partimbene et al. | Asynchronous multi-splitting method for linear and pseudo-linear problems | |
| Aigner et al. | Approximate implicitization of space curves | |
| Dasgupta | Locking-free compressible quadrilateral finite elements: Poisson’s ratio-dependent vector interpolants | |
| Wang et al. | A consistent projection integration for Galerkin meshfree methods | |
| Weng et al. | Substructuring method for eigensolutions | |
| Dachsel et al. | Shape matching by time integration of partial differential equations | |
| Schuster et al. | A mesh-free parallel moving least-squares-based interpolation method for the application in aeroelastic simulations with the flow simulator | |
| Dachsel et al. | The classic wave equation can do shape correspondence | |
| Roque et al. | RBF-FD meshless optimization using direct search (GLODS) in the analysis of composite plates | |
| Kim et al. | A new hybrid interpolation method using surface tracking, fitting and smoothing function applied for aeroelasticity | |
| Wilczewska et al. | Displacement Field Estimation for Echocardiography Strain Imaging Using B-Spline Based Elastic Image Registration—Synthetic Data Study | |
| JP7096218B2 (ja) | 情報処理装置、情報処理方法、及びプログラム | |
| Dachsel et al. | A Study of Spectral Expansion for Shape Correspondence |
Legal Events
| Date | Code | Title | Description |
|---|---|---|---|
| A521 | Written amendment |
Free format text: JAPANESE INTERMEDIATE CODE: A523 Effective date: 20061227 |
|
| A621 | Written request for application examination |
Free format text: JAPANESE INTERMEDIATE CODE: A621 Effective date: 20080122 |
|
| A977 | Report on retrieval |
Free format text: JAPANESE INTERMEDIATE CODE: A971007 Effective date: 20090129 |
|
| A131 | Notification of reasons for refusal |
Free format text: JAPANESE INTERMEDIATE CODE: A131 Effective date: 20090616 |
|
| RD03 | Notification of appointment of power of attorney |
Free format text: JAPANESE INTERMEDIATE CODE: A7423 Effective date: 20090806 |
|
| A521 | Written amendment |
Free format text: JAPANESE INTERMEDIATE CODE: A523 Effective date: 20090813 |
|
| A521 | Written amendment |
Free format text: JAPANESE INTERMEDIATE CODE: A821 Effective date: 20090806 |
|
| A02 | Decision of refusal |
Free format text: JAPANESE INTERMEDIATE CODE: A02 Effective date: 20091124 |
|
| A521 | Written amendment |
Free format text: JAPANESE INTERMEDIATE CODE: A523 Effective date: 20100223 |
|
| A911 | Transfer to examiner for re-examination before appeal (zenchi) |
Free format text: JAPANESE INTERMEDIATE CODE: A911 Effective date: 20100303 |
|
| A912 | Re-examination (zenchi) completed and case transferred to appeal board |
Free format text: JAPANESE INTERMEDIATE CODE: A912 Effective date: 20100416 |
|
| A521 | Written amendment |
Free format text: JAPANESE INTERMEDIATE CODE: A523 Effective date: 20110728 |
|
| A01 | Written decision to grant a patent or to grant a registration (utility model) |
Free format text: JAPANESE INTERMEDIATE CODE: A01 |
|
| R150 | Certificate of patent or registration of utility model |
Free format text: JAPANESE INTERMEDIATE CODE: R150 |