図面 (/)

技術 有限要素法直接計算プログラムおよび解析方法

出願人 国立大学法人東北大学
発明者 浦川肇十文字正樹
出願日 2005年5月6日 (15年7ヶ月経過) 出願番号 2005-134797
公開日 2006年11月16日 (14年1ヶ月経過) 公開番号 2006-313400
状態 特許登録済
技術分野 複合演算 CAD
主要キーワード 近似方程式 有限要素法解析プログラム 集中質量 境界値問題 基準要素 荷重ベクトル 連結ファイル 三角形要素
関連する未来課題
重要な関連分野

この項目の情報は公開日時点(2006年11月16日)のものです。
また、この項目は機械的に抽出しているため、正しく解析できていない場合があります

図面 (20)

課題

要素分割データ入力に際し、データ入力作業から各要素の節点局所番号付けする作業をなくし、高次要素を用いた計算であっても、データを与えるとほとんど誤差なく短時間に、剛性行列Kと質量行列Mを直接計算して求め、固有関数を精度良く可視化して出力することができる有限要素法直接計算プログラムおよび解析方法を提供することを目的とする。

解決手段

有限要素法により解析対象物数値モデルを作成し、数値解析を行い、解析結果を可視化するプログラムの解析方法において、有限要素法の基底関数から決まる剛性行列Kと質量行列Mを直接計算することにより、与えられた領域のラプラシアン固有値問題固有値と固有関数を求める。

概要

背景

有限要素法では、連続体である物体有限の形状と寸法を持つ要素の集合により近似して扱う。従来の有限要素法解析プログラムでは、解析対象物に適した数値解析を行うために、要素の形状、節点座標、および各要素の局所的な節点番号を定め、これらの要素分割のデータを入力してから、剛性行列Kと質量行列Mを計算する必要があった。

例えば、解析対象物が2次元の平面のモデルである場合には、剛性行列Kと質量行列Mは、[数1]と[数2]の定義に従って計算させるプログラムになっている。

2次元の平面領域三角形分割εに対する剛性行列K=(Kij)は、φi を点Pi における基底関数とするとき、次のように定義される。





ただし、G(ε)は要素の和集合である。

2次元の平面領域の三角形分割εに対する質量行列M=(Mij)は、φi を点Pi における基底関数とするとき、次のように定義される。




ただし、G(ε)は要素の和集合である。

また、これらの式は、3次元の中空体においても同様に定義され、解析対象物が3次元の中空体のモデルである場合には、剛性行列Kと質量行列Mは、[数3]と[数4]の定義に従って計算させるプログラムになっている。

3次元の中空体の三角形分割εに対する剛性行列K=(Kij)は、φi を点Pi における基底関数とするとき、次のように定義される。




ただし、G(ε)は要素の和集合であり、〈∇φi ,∇φj 〉g は中空体に誘導される計量ベクトルgに関する内積であり、vg
体積要素を表す。

3次元の中空体の三角形分割εに対する質量行列M=(Mij)は、φi を点Pi における基底関数とするとき、次のように定義される。




ただし、G(ε)は要素の和集合であり、vgは体積要素を表す。

さらに、解析対象物が3次元の中実体のモデルである場合には、剛性行列Kと質量行列Mは、[数5]と[数6]の定義に従って計算させるプログラムになっている。

3次元の中実体の四面体分割εに対する剛性行列K=(Kij)は、φi を点Pi における基底関数とするとき、次のように定義される。




ただし、G(ε)は要素の和集合である。

3次元の中実体の四面体分割εに対する質量行列M=(Mij)は、φi を点Pi における基底関数とするとき、次のように定義される。




ただし、G(ε)は要素の和集合である。

上記従来技術では、計算速度を早くするために要素分割を工夫し、その他有限要素法の使用に当たって留意した場合であっても、計算結果を出力させるまでに、多くのステップの計算を必要とした。このため、多くの時間がかかる上に誤差が累積し、出力された結果は不正確なものになるという問題があった。
また、剛性行列Kと質量行列Mの計算に関し、直接剛性法により、要素ごとに積分計算することを全要素について行い、全体的な近似方程式を組み立てるという従来技術もあるが、データ入力に際し、各要素の節点番号付け順序や、局所的な節点番号と全体的な節点番号の対応付けが必要になるなどの制約がある(例えば、非特許文献1参照。)。
さらに、有限要素法解析の精度を高めるために要素を小さくすると、要素数や節点数が増加し、データ入力とステップの計算に多くの作業と時間を要することになる。とくに、解析対象物が3次元中実体の複雑なモデルである場合には、要素数が数千万以上の規模になることがあり、剛性行列Kと質量行列Mの計算量が増加する。このため多くの時間がかかる上に誤差が累積し、出力された結果は不正確なものになる。さらに、剛性行列Kと質量行列Mの計算量が膨大になると、大きな記憶容量と高速計算処理を行うことができる高性能計算機が必要になり、一般に使用されている計算機では対応できなくなるという問題が生じる。

このため、従来から、剛性行列Kと質量行列Mの計算に関心が集まり、解析対象物の要素分割を工夫し、計算機における計算処理を分散化させるなどの改良が進められている。例えば、複数のプロセッサからなり、各プロセッサにメモリを持つようなマルチプロセッサ上での並列処理方式を採用し、各プロセッサが可能な限り、格納する剛性行列の成分、荷重ベクトルの成分を計算するのに必要となる要素剛性行列、応力を計算するように各プロセッサへの割り当てを行うことにより、バッファの容量を削減する有限要素法の処理方式(例えば、特許文献1参照。)が開示されている。また、要素分割の対象となる構造物亀裂先端部分の作成位置分割数を指定するだけで自動的に要素の頂点とその中間点の節点の座標データを作成し、前記亀裂先端部の基準要素に属する節点番号の集合を有する要素データを求め、前記基準要素の節点番号をもとにして他の要素の節点番号を求めて、その要素の節点番号の集合を有する要素データを次々に求める有限要素法解析用モデルの作成方法(例えば、特許文献2参照。)、有限要素法による解析対象物体の二次元又は三次元の有限要素データの作成方法において、要素の分割数の少ない方向から優先的に節点番号を付ける工程を備えた有限要素データの作成方法(例えば、特許文献3参照。)が開示されている。さらに、有限要素法の定型的手順により構造物の解析モデルから剛性行列および集中質量行列を生成し、これらの行列で表される一般固有値問題標準的固有値問題Ay=λyに変換し、集中質量行列が負の成分を含む場合、アーノルディ反復法によりAを射影し行列T={hij}を得るが、このときの初期ベクトルを、対応する集中質量行列の成分が正であるような成分hに対しては実数、負であるような成分hに対しては純虚数となるように設定することにより、行列Aを実ヘッセンベルグ行列Tに変換し、このヘッセンベルグ行列Tが与える標準的固有値問題の固有値QR法により求める方法およびシステム(例えば、特許文献4参照。)が開示されている。

地文雄著「有限要素法概説」1980年、P.58−67、株式 会社サイエンス
特開平05−108696号公報
特開平11−272649号公報
特開平11−283050号公報
特開2001−256216号公報

概要

要素分割のデータ入力に際し、データ入力作業から各要素の節点の局所番号付けする作業をなくし、高次要素を用いた計算であっても、データを与えるとほとんど誤差なく短時間に、剛性行列Kと質量行列Mを直接計算して求め、固有関数を精度良く可視化して出力することができる有限要素法直接計算プログラムおよび解析方法を提供することを目的とする。 有限要素法により解析対象物の数値モデルを作成し、数値解析を行い、解析結果を可視化するプログラムの解析方法において、有限要素法の基底関数から決まる剛性行列Kと質量行列Mを直接計算することにより、与えられた領域のラプラシアンの固有値問題の固有値と固有関数を求める。

目的

しかしながら、特許文献1の開示技術では、解析対象物が3次元の中実体の複雑なモデルになると、要素数が数千万以上の規模になることがあるため、計算処理方法を改良しても、計算機では処理できない膨大な計算量になる可能性がある。また、特許文献2ないし4のいずれの開示技術でも、本発明のように剛性行列の要素Kijを決定する式そのものを改良することにより、計算時間を大幅に短縮し、解析精度を向上させるという新たな提案はなされていない。
本発明は上記事情に鑑みてなされたものであり、要素分割のデータ入力に際し、データ入力作業から各要素の節点の局所番号付けする作業をなくし、高次要素を用いた計算であっても、データを与えるとほとんど誤差なく短時間に、剛性行列Kと質量行列Mを決定する式により計算して求め(以下、「直接計算」と呼ぶ。)、固有関数を精度良く可視化して出力することができる有限要素法直接計算プログラムおよび解析方法を提供することを目的とする。

効果

実績

技術文献被引用数
0件
牽制数
0件

この技術が所属する分野

ライセンス契約や譲渡などの可能性がある特許掲載中! 開放特許随時追加・更新中 詳しくはこちら

請求項1

有限要素法により解析対象物数値モデルを作成し、数値解析を行い、解析結果を可視化するプログラムにおいて、前記プログラムは有限要素法の基底関数から決まる剛性行列Kと質量行列Mを直接計算することにより、与えられた領域のラプラシアン固有値問題固有値固有関数を求めることを特徴とするプログラム。

請求項2

前記解析対象物が1次元モデルである場合に、解析対象物を1次元モデルによる解析方法より解析し、前記解析対象物が2次元モデルである場合に、解析対象物を1次元モデルによる解析方法、2次元モデルによる解析方法又はこれらの組み合わせにより解析し、前記解析対象物が3次元モデルである場合に、前記解析対象物を、1次元モデルによる解析方法、2次元モデルによる解析方法、3次元モデルによる解析方法又はこれらの任意の組み合わせにより解析することを特徴とする請求項1に記載のプログラム。

請求項3

有限要素法により解析対象物の数値モデルを作成し、数値解析を行い、解析結果を可視化する解析方法において、前記解析方法は有限要素法の基底関数から決まる剛性行列Kと質量行列Mを直接計算することにより、与えられた領域のラプラシアンの固有値問題の固有値と固有関数を求めることを特徴とする解析方法。

請求項4

前記解析対象物が1次元モデルである場合に、解析対象物を1次元モデルによる解析方法より解析し、前記解析対象物が2次元モデルである場合に、解析対象物を1次元モデルによる解析方法、2次元モデルによる解析方法又はこれらの組み合わせにより解析し、前記解析対象物が3次元モデルである場合に、前記解析対象物を、1次元モデルによる解析方法、2次元モデルによる解析方法、3次元モデルによる解析方法又はこれらの任意の組み合わせにより解析することを特徴とする請求項3に記載の解析方法。

請求項5

前記解析方法は、剛性行列Kと質量行列Mを直接計算する場合に、前記解析対象物が2次元の平面のモデルであるとき、剛性行列の要素Kijについて、の定義により、剛性行列の要素Kijを決定する式が、i≠jのとき、 i=jのとき、 のように得られることを特徴とする請求項3又は4に記載の解析方法。

請求項6

前記解析方法は、剛性行列Kと質量行列Mを直接計算する場合に、前記解析対象物が2次元の平面のモデルであるとき、質量行列の要素Mijについて、の定義により、質量行列の要素Mijを決定する式が、i≠jのとき、i=jのとき、 のように得られることを特徴とする請求項3又は4に記載の解析方法。

請求項7

前記解析方法は、剛性行列Kと質量行列Mを直接計算する場合に、前記解析対象物が3次元の中空体のモデルであるとき、剛性行列の要素Kijについて、の定義により、剛性行列の要素Kijを決定する式が、i≠jのとき、 i=jのとき、 のように得られることを特徴とする請求項3又は4に記載の解析方法。

請求項8

前記解析方法は、剛性行列Kと質量行列Mを直接計算する場合に、前記解析対象物が3次元の中空体のモデルであるとき、質量行列の要素Mijについて、 の定義により、質量行列の要素Mijを決定する式が、i≠jのとき、i=jのとき、 のように得られることを特徴とする請求項3又は4に記載の解析方法。

請求項9

前記解析方法は、剛性行列Kと質量行列Mを直接計算する場合に、前記解析対象物が3次元の中実体のモデルであるとき、剛性行列の要素Kijについて、の定義により、剛性行列の要素Kijを決定する式が、i≠jのとき、 i=jのとき、 のように得られることを特徴とする請求項3又は4に記載の解析方法。

請求項10

前記解析方法は、剛性行列Kと質量行列Mを直接計算する場合に、前記解析対象物が3次元の中実体のモデルであるとき、質量行列の要素Mijについて、の定義により、質量行列の要素Mijを決定する式が、i≠jのとき、 i=jのとき、 のように得られることを特徴とする請求項3又は4に記載の解析方法。

請求項11

前記解析方法は、剛性行列Kと質量行列Mを直接計算することにより、与えられた領域のラプラシアンの固有値問題の固有値と固有関数を求める場合に、固有値問題を、 で与えたとき、第k番目固有値λkに対応する固有ベクトルukの第i成分を、第k番目の固有関数φk の第i番目節点の値とし、隣接する節点間の値は線形補間によって定めることにより第k番目の固有関数φkの値を求め、その固有関数φkの値を用いて解析対象物を可視化することを特徴とする請求項3ないし10のいずれかに記載の解析方法。

技術分野

0001

本発明は、有限要素法により解析対象物解析する場合に、剛性行列Kと質量行列Mの計算時間を短縮し、精度よく固有関数可視化することができる有限要素法直接計算プログラムおよび解析方法に関する。

背景技術

0002

有限要素法では、連続体である物体有限の形状と寸法を持つ要素の集合により近似して扱う。従来の有限要素法解析プログラムでは、解析対象物に適した数値解析を行うために、要素の形状、節点座標、および各要素の局所的な節点番号を定め、これらの要素分割のデータを入力してから、剛性行列Kと質量行列Mを計算する必要があった。

0003

例えば、解析対象物が2次元の平面のモデルである場合には、剛性行列Kと質量行列Mは、[数1]と[数2]の定義に従って計算させるプログラムになっている。

0004

2次元の平面領域三角形分割εに対する剛性行列K=(Kij)は、φi を点Pi における基底関数とするとき、次のように定義される。





ただし、G(ε)は要素の和集合である。

0005

2次元の平面領域の三角形分割εに対する質量行列M=(Mij)は、φi を点Pi における基底関数とするとき、次のように定義される。




ただし、G(ε)は要素の和集合である。

0006

また、これらの式は、3次元の中空体においても同様に定義され、解析対象物が3次元の中空体のモデルである場合には、剛性行列Kと質量行列Mは、[数3]と[数4]の定義に従って計算させるプログラムになっている。

0007

3次元の中空体の三角形分割εに対する剛性行列K=(Kij)は、φi を点Pi における基底関数とするとき、次のように定義される。




ただし、G(ε)は要素の和集合であり、〈∇φi ,∇φj 〉g は中空体に誘導される計量ベクトルgに関する内積であり、vg
体積要素を表す。

0008

3次元の中空体の三角形分割εに対する質量行列M=(Mij)は、φi を点Pi における基底関数とするとき、次のように定義される。




ただし、G(ε)は要素の和集合であり、vgは体積要素を表す。

0009

さらに、解析対象物が3次元の中実体のモデルである場合には、剛性行列Kと質量行列Mは、[数5]と[数6]の定義に従って計算させるプログラムになっている。

0010

3次元の中実体の四面体分割εに対する剛性行列K=(Kij)は、φi を点Pi における基底関数とするとき、次のように定義される。




ただし、G(ε)は要素の和集合である。

0011

3次元の中実体の四面体分割εに対する質量行列M=(Mij)は、φi を点Pi における基底関数とするとき、次のように定義される。




ただし、G(ε)は要素の和集合である。

0012

上記従来技術では、計算速度を早くするために要素分割を工夫し、その他有限要素法の使用に当たって留意した場合であっても、計算結果を出力させるまでに、多くのステップの計算を必要とした。このため、多くの時間がかかる上に誤差が累積し、出力された結果は不正確なものになるという問題があった。
また、剛性行列Kと質量行列Mの計算に関し、直接剛性法により、要素ごとに積分計算することを全要素について行い、全体的な近似方程式を組み立てるという従来技術もあるが、データ入力に際し、各要素の節点番号付け順序や、局所的な節点番号と全体的な節点番号の対応付けが必要になるなどの制約がある(例えば、非特許文献1参照。)。
さらに、有限要素法解析の精度を高めるために要素を小さくすると、要素数や節点数が増加し、データ入力とステップの計算に多くの作業と時間を要することになる。とくに、解析対象物が3次元中実体の複雑なモデルである場合には、要素数が数千万以上の規模になることがあり、剛性行列Kと質量行列Mの計算量が増加する。このため多くの時間がかかる上に誤差が累積し、出力された結果は不正確なものになる。さらに、剛性行列Kと質量行列Mの計算量が膨大になると、大きな記憶容量と高速計算処理を行うことができる高性能計算機が必要になり、一般に使用されている計算機では対応できなくなるという問題が生じる。

0013

このため、従来から、剛性行列Kと質量行列Mの計算に関心が集まり、解析対象物の要素分割を工夫し、計算機における計算処理を分散化させるなどの改良が進められている。例えば、複数のプロセッサからなり、各プロセッサにメモリを持つようなマルチプロセッサ上での並列処理方式を採用し、各プロセッサが可能な限り、格納する剛性行列の成分、荷重ベクトルの成分を計算するのに必要となる要素剛性行列、応力を計算するように各プロセッサへの割り当てを行うことにより、バッファの容量を削減する有限要素法の処理方式(例えば、特許文献1参照。)が開示されている。また、要素分割の対象となる構造物亀裂先端部分の作成位置分割数を指定するだけで自動的に要素の頂点とその中間点の節点の座標データを作成し、前記亀裂先端部の基準要素に属する節点番号の集合を有する要素データを求め、前記基準要素の節点番号をもとにして他の要素の節点番号を求めて、その要素の節点番号の集合を有する要素データを次々に求める有限要素法解析用モデルの作成方法(例えば、特許文献2参照。)、有限要素法による解析対象物体の二次元又は三次元の有限要素データの作成方法において、要素の分割数の少ない方向から優先的に節点番号を付ける工程を備えた有限要素データの作成方法(例えば、特許文献3参照。)が開示されている。さらに、有限要素法の定型的手順により構造物の解析モデルから剛性行列および集中質量行列を生成し、これらの行列で表される一般固有値問題標準的固有値問題Ay=λyに変換し、集中質量行列が負の成分を含む場合、アーノルディ反復法によりAを射影し行列T={hij}を得るが、このときの初期ベクトルを、対応する集中質量行列の成分が正であるような成分hに対しては実数、負であるような成分hに対しては純虚数となるように設定することにより、行列Aを実ヘッセンベルグ行列Tに変換し、このヘッセンベルグ行列Tが与える標準的固有値問題の固有値QR法により求める方法およびシステム(例えば、特許文献4参照。)が開示されている。

0014

地文雄著「有限要素法概説」1980年、P.58−67、株式 会社サイエンス
特開平05−108696号公報
特開平11−272649号公報
特開平11−283050号公報
特開2001−256216号公報

発明が解決しようとする課題

0015

しかしながら、特許文献1の開示技術では、解析対象物が3次元の中実体の複雑なモデルになると、要素数が数千万以上の規模になることがあるため、計算処理方法を改良しても、計算機では処理できない膨大な計算量になる可能性がある。また、特許文献2ないし4のいずれの開示技術でも、本発明のように剛性行列の要素Kijを決定する式そのものを改良することにより、計算時間を大幅に短縮し、解析精度を向上させるという新たな提案はなされていない。
本発明は上記事情に鑑みてなされたものであり、要素分割のデータ入力に際し、データ入力作業から各要素の節点の局所番号付けする作業をなくし、高次要素を用いた計算であっても、データを与えるとほとんど誤差なく短時間に、剛性行列Kと質量行列Mを決定する式により計算して求め(以下、「直接計算」と呼ぶ。)、固有関数を精度良く可視化して出力することができる有限要素法直接計算プログラムおよび解析方法を提供することを目的とする。

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

0016

上記課題を解決するための手段として、本発明は以下の特徴を有している。
本発明のプログラムは、有限要素法により解析対象物の数値モデルを作成し、数値解析を行い、解析結果を可視化するプログラムにおいて、有限要素法の基底関数から決まる剛性行列Kと質量行列Mを直接計算することにより、与えられた領域のラプラシアンの固有値問題の固有値と固有関数を求めることを特徴とする。
本発明によれば、剛性行列の要素Kijを決定する式そのものを改良することにより、要素分割のデータ入力に際し、データ入力作業から各要素の節点の局所番号付けを省くことができ、データを与えるとほとんど誤差なく短時間に剛性行列Kと質量行列Mを直接計算して求め、固有値と固有関数を出力させることができる。本発明では、各要素の節点のつながりが把握できればよいため、従来必要とされていた各要素の節点番号付けの順序や、局所的な節点番号と全体的な節点番号の対応付けは不要になる。設定できる要素数と節点数の制限をなくしたことにより、1次元の線分、2次元の平面はもちろん、3次元の中空体、又は3次元の中実体の複雑な解析対象物の高次要素を用いた計算にも対応できるようにしている。

0017

本発明のプログラムは、前記解析対象物が1次元モデルである場合に、解析対象物を1次元モデルによる解析方法より解析し、前記解析対象物が2次元モデルである場合に、解析対象物を1次元モデルによる解析方法、2次元モデルによる解析方法又はこれらの組み合わせにより解析し、前記解析対象物が3次元モデルである場合に、前記解析対象物を、1次元モデルによる解析方法、2次元モデルによる解析方法、3次元モデルによる解析方法又はこれらの任意の組み合わせにより解析することを特徴とする。

0018

本発明の解析方法は、有限要素法により解析対象物の数値モデルを作成し、数値解析を行い、解析結果を可視化する解析方法において、有限要素法の基底関数から決まる剛性行列Kと質量行列Mを直接計算することにより、与えられた領域のラプラシアンの固有値問題の固有値と固有関数を求めることを特徴とする。

0019

本発明の解析方法は、前記解析対象物が1次元モデルである場合に、解析対象物を1次元モデルによる解析方法より解析し、前記解析対象物が2次元モデルである場合に解析対象物を1次元モデルによる解析方法、2次元モデルによる解析方法又はこれらの組み合わせにより解析し、前記解析対象物が3次元モデルである場合に、前記解析対象物を、1次元モデルによる解析方法、2次元モデルによる解析方法、3次元モデルによる解析方法又はこれらの任意の組み合わせにより解析することを特徴とする。

0020

本発明の解析方法は、剛性行列Kと質量行列Mを直接計算する場合に、前記解析対象物が2次元の平面のモデルであるとき、剛性行列の要素Kijについて、


の定義により、剛性行列の要素Kijを決定する式が、
i≠jのとき、


i=jのとき、


のように得られることを特徴とする。
本発明の剛性行列の要素Kijを決定する式は、式そのものを改良したものであり、2次元の平面領域の解析対象物を三角形要素εに分割し、三角形の各線分の長さを分子に、それに対応する三角形の面積分母にした量で定式化したものである。これらの式は、3次元の中空体においても、同様に適用することができる。
図8は、2次元平面領域および3次元中空体における、i≠jのときの剛性行列Kijを決定する式を説明するための概念図である。節点Pi と節点Pj が異なるときは、図8のようになる。上記i≠jのときのKijの式の分母は対応する三角形の面積を表し、area(eμ)は三つの節点Pi ,Pj ,Pkからなる三角形を、area(ev) は三つの節点Pi ,Pj , Pl からなる三角形を表す。分子はそれぞれの辺の長さの2乗などを足したり引いたりしたものである。
図9は、2次元平面領域および3次元中空体における、i­=jのときの剛性行列Kiiを決定する式を説明するための概念図である。 節点Pi と節点Pj が同じときは、図9のようになる。節点Pi
に接続する節点と三角形が書かれている。上記i­=jのときのKiiの式は、各三角形のPi を含まない辺の長さを2乗し、三角形の面積で割った量の和で与えられる。
従って、局所的な節点番号と全体的な節点番号を対応づけは不要となる。3次元の中空体の場合も同様である。

0021

本発明の解析方法は、剛性行列Kと質量行列Mを直接計算する場合に、前記解析対象物が2次元の平面のモデルであるとき、質量行列の要素Mijについて、


の定義により、質量行列の要素Mijを決定する式が、
i≠jのとき、


i=jのとき、


のように得られることを特徴とする。
本発明の質量行列の要素Mijを決定する式は、3次元の中空体においても、同様に適用することができる。

0022

本発明の解析方法は、剛性行列Kと質量行列Mを直接計算する場合に、前記解析対象物が3次元の中空体のモデルであるとき、剛性行列の要素Kijについて、


の定義により、剛性行列の要素Kijを決定する式が、
i≠jのとき、


i=jのとき、


のように得られることを特徴とする。

0023

本発明の解析方法は、剛性行列Kと質量行列Mを直接計算する場合に、前記解析対象物が3次元の中空体のモデルであるとき、質量行列の要素Mijについて、


の定義により、質量行列の要素Mijを決定する式が、
i≠jのとき、


i=jのとき、


のように得られることを特徴とする。

0024

本発明の解析方法は、剛性行列Kと質量行列Mを直接計算する場合に、前記解析対象物が3次元の中実体のモデルであるとき、剛性行列の要素Kijについて、


の定義により、剛性行列の要素Kijを決定する式が、
i≠jのとき、


i=jのとき、


のように得られることを特徴とする。
本発明の剛性行列の要素Kijを決定する式は、式そのものを改良したものであり、3次元の中実体の解析対象物を四面体要素eに分割し、四面体の各線分の長さを分子に、それに対応する四面体の体積を分母にし、四面体要素eすべてにわたる和として定式化したものである。

0025

本発明の解析方法は、剛性行列Kと質量行列Mを直接計算する場合に、前記解析対象物が3次元の中実体のモデルであるとき、質量行列の要素Mijについて、


の定義により、質量行列の要素Mijを決定する式が、
i≠jのとき、


i=jのとき、


のように得られることを特徴とする。
本発明の質量行列の要素Mijを決定する式は、四面体要素eすべてにわたる和として定式化したものであり、Vol(e)は四面体の体積を表す。

0026

本発明の解析方法は、剛性行列Kと質量行列Mを直接計算することにより、与えられた領域のラプラシアンの固有値問題の固有値と固有関数を求める場合に、固有値問題を、


で与えたとき、第k番目固有値λ
に対応する固有ベクトルukの第i成分を、第k番目の固有関数φk の第i番目の節点の値とし、隣接する節点間の値は線形補間によって定めることにより第k番目の固有関数φk
の値を求め、その固有関数φk の値を用いて解析対象物を可視化することを特徴とする。
本発明の可視化方法は、剛性行列Kと質量行列Mを直接計算することにより、固有空間所属する第k番目の固有関数φkの値を用いて解析対象物を可視化して短時間に精度良く出力することができる。
図10は、解析対象物が2次元平面および3次元中空体のモデルであるときの、隣接する節点間の値を線形補間によって定める方法の説明をするための概念図である。

発明の効果

0027

上記課題を解決するための手段により、本発明によれば、有限要素法により対象物を解析する場合に、要素分割のデータ入力に際し、データ入力作業から各要素の節点の局所番号付けを省くことができた。また、従来必要とされていた各要素の節点番号付けの順序や局所的な節点番号と全体的な節点番号の対応付けを不要にしたことによって、設定できる要素数と節点数の制限をなくすことができた。本発明によれば、1次元の線分、2次元の平面はもちろん、3次元の中空体、又は3次元の中実体の複雑な解析対象物を解析する場合であっても、データを与えると、従来技術に比べて、ほとんど誤差なく短時間に、剛性行列Kと質量行列Mを直接計算して求めることができた。
また、剛性行列Kと質量行列Mを直接計算した結果を用いて、ディリクレ境界値問題ノイマン境界値固有問題の固有値と固有関数、および高次の固有値と固有関数を容易かつ精度良く可視化して出力することができることから、定量的で具体的な解析をすることができ、取り扱いが容易で実用的な有限要素法直接計算プログラムおよび解析方法を提供することができた。

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

0028

本発明を実施するための最良の形態を、図1から図3に基づいて説明する。図1から図3では、連続体である3次元の中実体の解析対象物を、有限の形状と寸法を持つ四面体要素eに分割し、各要素の集合により近似して扱う。
ただし、これらは一実施形態にすぎず、本発明の特許請求の範囲を限定するものではない。

0029

図1は、有限要素法直接計算プログラムの解析方法の手順を示したフローチャートである。ステップS01〜S14により剛性行列Kと質量行列Mを決定して出力し、ステップS15〜S16により固有値と固有関数を求めて出力する手順により、解析する。

0030

図2は、図4における領域Ω1 の四面体要素の35個の節点の座標を与える、節点の座標ファイルの例である。節点とは各四面体要素の頂点をいう。図2の節点の座標ファイルは、前もって準備しておく。節点数を1行目に入力し、2行目以降は節点の座標を1行ずつ入力する。この順番が節点番号となるので、節点番号を与える必要はない。座標は、左からx座標、y座標およびz座標とする。例えば、1行目は節点数が35であり、2行目は節点番号1のx座標が0.5、y座標が0.5、z座標が0.5であることを意味する。

0031

図3は、図4における領域Ω1
の96個の四面体要素の節点の連結関係を与える、連結ファイルの例である。連結ファイルとは、各要素をなす節点の番号のファイルをいう。節点番号は、35個の節点P1
〜P35 の添字を用いて1〜35のように表す。図3の連結ファイルは、前もって準備しておく。図3の連結ファイルは、連結関係にある4点の節点番号の組の総数を1行目に入力し、2行目以降は連結関係にある4点の節点番号の組を1行ずつ入力する。この場合、四面体をなす節点番号が分かればよく、時計、反時計まわりのどちらの順番で入力してもよい。図3の2行目は節点番号16,35,23,2からなる四面体の要素を表している。

0032

本発明の剛性行列の要素Kijと質量行列の要素Mijを決定する式の計算方法について述べる。

0033

i≠jのとき、剛性行列Kと質量行列Mの(i,j)成分は、[数7]、[数8]および[数9]により計算される。ここで、i、j、k、lは自然数である。

0034

節点Pの座標を(x(P),y(P),z(P))で表したとき、始点P、終点QのベクトルPQと始点R、終点SのベクトルRSの内積は次のように表される。

0035

ここで、右辺の和は辺Pi
Pj を含む四面体eすべてにわたる和であり、四面体eの4つの節点Pi ,Pj ,Pk (e),Pl (e)で表す。また、PiPl (e)等は、始点をPi 、終点をPl (e)とするベクトルを表し、(PiPl (e),Pj Pl (e))は2つのベクトルPiPl (e)とPj Pl (e)との、[数7]における内積を表す。Vol(e)は四面体eの体積を表している。

0036

ここで、右辺の和は辺Pi
Pj を含む四面体eすべてにわたる和であり、Vol(e)は四面体eの体積を表している。

0037

i=jのとき、剛性行列Kと質量行列Mの(i,i)成分は、[数7]、[数10]および[数11]により計算される。ここで、i、j、k、lは自然数である。

0038

ここで、右辺の和は節点Pi
を含む四面体eすべてにわたる和である。Vol(e)は四面体eの体積を表す。

0039

ここで、右辺の和は節点Pi
を含む四面体eすべてにわたる和である。Vol(e)は四面体eの体積を表す。

0040

以下、図1に示すステップS01〜S14により剛性行列Kと質量行列Mを決定して出力する解析方法について説明する。
ステップS01では、節点の座標ファイルを読み込む。節点の座標ファイルは、図2のように前もって準備しておく。
ステップS02では、連結ファイルを読み込む。連結ファイルは、図3のように前もって準備しておく。
ステップS03では、L×L行列のすべての成分をとする初期化を行っている。L×L行列は、解決すべき問題に対し、必要な行列とする。ただし、Lは1≦L≦Nの定数とし、Lはプログラムの中で与える。ここで、Nは自然数であり、節点数を表す。

0041

ステップS04〜S07では、i≠jのときの、剛性行列の要素Kijと質量行列の要素Mijを決定する式の計算をしている。
ステップS04では、ステップS02で読み込んだ節点Pi
およびPj(i<j)をもつ四面体eの残りの節点がPk ,Pl ならば、ステップS01で読み込んだ各節点の座標(x(i),y(i),z(i))、(x(j),y(j),z(j))、(x(k),y(k),z(k))、(x(l),y(l),z(l))を用いて、内積(PiPl,PjPl)、(PiPk,PjPk)、(PiPl,PjPk)、(PiPk,PjPl)を計算する。
ステップS05では、四面体eの体積Vol(e)を求める。
ステップS06では、AK(i,j)、AM(i,j)を求める。
[数8]より


の値をAK(i,j)、
[数7]より


の値をAM(i,j)とする。
ステップS07では、ステップS06と同様に、四面体eの辺Pi
Pj 以外についても計算する。具体的には、
AK(i,k)、AK(i,l)、AK(j,k)、AK(j,l)、AK(k,l)、
AM(i,k)、AM(i,l)、AM(j,k)、AM(j,l)、AM(k,l)
である。

0042

ステップS08〜S11では、i=jのときの、剛性行列の要素Kijと質量行列の要素Mijを決定する式の計算をしている。
ステップS08では、節点Pi
を含む四面体eの残りの節点がPj ,Pk ,Pl ならばステップS01で読み込んだ各節点の座標(x(i),y(i),z(i))、(x(j),y(j),z(j))、(x(k),y(k),z(k))(x(l),y(l),z(l))を用いて、内積(Pj Pl,Pj
Pk )と‖PjPl‖2、‖PjPk‖2
を計算する。
ステップS09では、四面体eの体積Vol(e)を求める。
ステップS10では、AK(i,i)、AM(i,i)を求める。
[数8]より


の値をAK(i,i)、
[数9]より


の値をAM(i,i)とする。
ステップS11では、ステップS10と同様に、四面体eの点Pi
以外についても計算する。具体的には、
AK(j,j)、AK(k,k)、AK(l,l)、AM(j,j)、AM(k,k)、
AM(l,l)
である。

0043

ステップS12では、すべての四面体要素に対してS04〜S11で行ったことを繰り返す。そして、各要素を加えた和集合が、行列KおよびMの要素となる。
ステップS13では、ステップS12で求めた行列KおよびMの成分を用いて、(j,i)成分(i>j)を求める。連結関係のない成分はすべて零であり、ステップS03で初期化した結果をそのまま用いる。以上で、行列KおよびMの成分が決定できる。
ステップS14では、行列KおよびMのすべての成分を出力する。このステップは省略してもよい。

0044

次に、ステップS15〜S16により固有値と固有関数を求めて出力することにより、解析対象物を可視化する解析方法について説明する。
ステップS15では、[数12]で表される行列の一般固有値問題を解く。ここでは既に使用されている行列の一般固有値問題を解くプログラムを組み込んで解いているが、ステップS14で出力してある行列KおよびMの成分を用いて別のプログラムを使用しても構わない。

0045

ただし、Δはラプラシアンであって、φのラプラシアンは次のように表される。

0046

[数12]における第k番目の固有値を、[数13]における第k番目の固有値λkと定める。また、[数12]の第k番目の固有値λk
に対応する固有ベクトルukから、[数13]における第k番目の固有関数φkを得るには、固有ベクトルukの第i成分を固有関数φk
の第i番目の節点の値とし、隣接する節点間の値を線形補間によって固有関数φk の値を定める。
すなわち、(x,y,z)∈四面体Pi Pj Pk Pl において、
φk
(x,y,z)=a×(uの第i成分)+b×(uの第j成分)+c×(uの第k成分)+d×(uの第l成分)
によって与える。
ただし、0≦a,b,c,d≦1、a+b+c+d≦1である。

0047

ステップS16では、[数13]における第k番目の固有値λk
、固有関数φk の出力を行う。ただし、kは1≦k≦Lの定数とする。
以上で、図1に示すステップS01〜S19の一連の処理を終了する。

0048

次に、本発明の剛性行列Kの計算時間について、従来の計算方法と比較するために、2次元の場合について説明する。
従来の計算方法については、非特許文献1に記載の剛性行列の計算法に従って演算回数を調べる。三角形要素の1つをP1 (x1 ,y1 ),P2
(x2 ,y2 ),P3 (x3 ,y3 )とし、反時計回りに並んでいるものとする。従来の方法では、3点がP1 ,P2 ,P3 の順に反時計回りに並んでいる必要があるが、本発明の方法では向きを考慮する必要はない。

0049

従来の計算方法によると、
(1)面積Sの2倍に等しい行列式Dの計算
D=2×S=x1×(y2
−y3 )+x2×(y3 −y1 )+x3×(y1
−y2 )
と計算できるので、加算2回、減算3回、乗算3回の合計8回の演算が必要となる。
(2)面積座標での近似関数係数の計算


は減算1回、除算1回でのそれぞれで行うので、合計12回の演算が必要となる。
(3)S=D/2より、除算1回が必要となる。
(4)要素係数行列の計算
(2)で求めた係数を用いて、Aij=S×(bj×bi+cj×ci)を計算する。
すなわち、A11 ,A22,A33 ,A12(=A21),A23(=A32),A31(=A13)のそれぞれで乗算3回、加算1回を6回行うので、合計24回の演算が必要となる。
(5)9つの要素を加える
合計9回の演算が必要となる。
上より、要素がNならば、(8+12+1+24+9)N=54Nの演算回数が必要となる。

0050

本発明の計算方法によると、
(1)面積Sの2倍に等しい行列式Dの計算
従来の方法と同様に8回の演算が必要となる。
(2)線分の2乗の計算
PiPj2 =(xi
−xj )2+(yi −yj)2
と計算できるので、加算1回、減算2回、乗算2回となり、
PjPk2,PkPi2のそれぞれで行うので、合計15回の演算が必要となる。
(3)S1 =−0.25/D、S2 =0.5/Dより、除算2回が必要となる。
(4)要素係数行列の計算

Aij=S1×(PiPk2+PjPk2−PiPj2)より、加算1回、減算1回、乗算1回
Aii=S2×PjPk2より、乗算1回を用いて、
A11 ,A22,A33 ,A12(=A21),A23(=A32),A31(=A13)のそれぞれを計算するので、合計12回の演算が必要となる。
(5)9つの要素を加える
合計9回の演算が必要となる。
以上より、要素がNならば、(8+15+2+12+9)N=46Nの演算回数が必要となる。

0051

従って、従来の計算方法と本発明の計算方法とを比較すると、
(1−46N/54N)×100=14.8%
となり、およそ15%の計算時間の短縮が見込まれる。
15%の計算時間の短縮により、従来の計算方法では3ギガバイトの記憶容量が必要であったものが、本発明の計算方法では2.5ギガバイトの記憶容量で足りることになる。経費にすると、およそ5、6万円の削減が可能になる。
また、要素数5000(節点2601)の場合、従来の方法では0.00138秒という結果が得られている。このことから、要素の数が数千万以上の規模となった場合、本発明の計算方法を適用すれば、計算時間の短縮に寄与することが可能になる。

0052

また、3次元の場合、中空体、中実体のいずれも、剛性行列Kの計算時間について厳密に計算ステップを記載してある従来の計算方法は見られないようであるが、3次元の中実体の複雑な解析対象物になると、従来の計算方法では計算に1日ほど要することが見込まれる。
しかしながら、本発明の計算方法を適用すれば、3次元の中空体、又は3次元の中実体の複雑な解析対象物を解析する場合であっても、データを与えると、従来の計算方法に比べて、ほとんど誤差なく短時間に、剛性行列Kと質量行列Mを直接計算することにより解析することができる。
例えば、節点2379、要素数8871の直方体の場合には、
(1)
行列KおよびMの決定に、84秒を要し、
(2)
固有値、固有ベクトルの決定に、440秒を要した。
以上より、84秒+440秒=524秒=8分44秒の計算時間で処理できる。
ただし、(1)(2)いずれもデータの入出力等の時間が含まれており、計算のみに比べれば時間は余分にかかっている。

0053

[実施例1]は、3次元の中実体の例である。
図4は、領域Ω1
における立方体領域を四面体要素に分割したときの節点および要素の状態を示した図である。
領域Ω1
={0<x<1,0<y<1,0<z<1}における
固有値問題:Δu=λu u=u(x,y,z);(x,y,z)∈Ω1 、
ディリクレ境界条件:u(x,y,z)=0 (x,y,z)∈∂Ω1
図4の四面体要素に分割することにより求める。この四面体分割の節点の座標と連結情報図2および図3で与えるものとする。
図5は、本発明の解析方法を図4に適用した場合の行列KおよびMの各成分の一覧である。
図6は、本発明の解析方法を図4に適用した場合の固有値および固有ベクトルの数値解である。
図7は、本発明の解析方法を図4に適用した場合の第1〜第5固有関数を可視化したものである。右図は左図の内部を取り出したものである。

0054

[実施例2]は、2次元平面領域の例である。
図11は、領域Ω1の節点および要素の状態を示した図である。
領域Ω1
={0<x<1,0<y<1}における
固有値問題:Δu=λu u=u(x,y);(x,y)∈Ω1 、
ディリクレ境界条件:u(x,y)=0 (x,y)∈∂Ω1
図11の三角形要素に分割することにより求める。

0055

[実施例3]は、2次元平面領域の例である。
領域Ω2 ついて三角形分割を行い、Ω2 の固有値問題およびディリクレ境界条件及びノイマン境界条件に対する固有値問題に対する第19の固有関数を計算し、可視化した結果を図15と16に示す。

0056

[実施例4]は、3次元の中空体の例である。
3次元中空体の曲面における固有関数の結果を図17〜19に示す。
固有値については、Ωを上面Ω+ と下面Ω− の2重に三角形分割して得られる、潰れた曲面とみなし、固有値を出力する。固有関数については、出力された固有関数の上面Ω+と下面Ω−
の節点の値の差を出力する。

0057

方程式に対して、温度分布を解析した実施例を示す。1m四方の領域の頂点を、それぞれ(0,0)、(1,0)、(1,1)、(0,1)とおく。初期条件として、u(t,x,y)=u(0,0.5,0.25)=104を与える。この領域の時刻t=0.02秒のディリクレ境界条件とノイマン境界条件の温度分布を図20〜23にそれぞれ示す。

0058

波動方程式に対して、波の動きを解析した実施例を示す。1m四方の領域の頂点をそれぞれ(0,0)、(1,0)、(1,1)、(0,1)とおく。初期条件として、u(t,x,y)=u(0,0.01,0.01)=u(0,0.99,0.99)=104 を与える。
図24は、この領域の時刻t=0.04秒のディリクレ境界条件の時の波の動きを表したものである。
図25は、この領域の時刻t=0.03秒のノイマン境界条件の時の波の動きを表したものである。

図面の簡単な説明

0059

有限要素法直接計算プログラムの解析方法の手順を示したフローチャートで ある。
図4における領域Ω1の四面体要素の35個の節点の座標を与える、節点の座標ファイルの例である。
図4における領域Ω1の96個の四面体要素の節点の連結関係を与える、連 結ファイルの例である。
領域Ω1における立方体領域を四面体要素に分割したときの節点および要素 の状態を示した図である。
本発明の解析方法を図4に適用した場合の行列KおよびMの各成分の一覧である。
本発明の解析方法を図4に適用した場合の固有値および固有ベクトルの数値解である。
本発明の解析方法を図4に適用した場合の第1〜第5固有関数を可視化したものである。
2次元平面領域および3次元中空体における、i­=jのときの剛性行列Kijを決定する式を説明するための概念図である。
2次元平面領域および3次元中空体における、i≠jのときの剛性行列Kiiを決定する式を説明するための概念図である。
解析対象物が2次元平面および3次元中空体のモデルであるときの、隣接する節点間の値を線形補間によって定める方法の説明をするための概念図である。
領域Ω1 の節点および要素の状態を示した図である。
本発明の解析方法を図11に適用した場合の行列KおよびMの各成分の一覧である。
本発明の解析方法を図11に適用した場合の固有値および固有ベクトルの数値解である。
領域Ω2 の節点および要素の状態を示した図である。
本発明の解析方法を図14に適用した場合のディリクレ境界条件における第19番目の固有関数を可視化した図である。
本発明の解析方法を図14に適用した場合のノイマン境界条件における第19番目の固有関数を可視化した図である。
本発明の解析方法を球面に適用した場合の第24番目の固有関数を可視化した図である。
本発明の解析方法をトーラスに適用した場合の第27番目の固有関数を可視化した図である。
本発明の解析方法を楕円面に適用した場合の第26番目の固有関数を可視化した図である。
本発明の解析方法を正方形領域に適用した場合のディリクレ境界条件における時刻t=0.02での温度分布の正面図である。
本発明の解析方法を正方形領域に適用した場合のディリクレ境界条件における時刻t=0.02での温度分布の立体図である。
本発明の解析方法を正方形領域に適用した場合のノイマン境界条件における時刻t=0.02での温度分布の正面図である。
本発明の解析方法を正方形領域に適用した場合のノイマン境界条件における時刻t=0.02での温度分布の立体図である。
本発明の解析方法を正方形領域に適用した場合のディリクレ境界条件における時刻t=0.04での波の動きを表す正面図である。
本発明の解析方法を正方形領域に適用した場合のノイマン境界条件における時刻t=0.03での波の動きを表す立体図である。

ページトップへ

この技術を出願した法人

この技術を発明した人物

ページトップへ

関連する挑戦したい社会課題

関連する公募課題

ページトップへ

技術視点だけで見ていませんか?

この技術の活用可能性がある分野

分野別動向を把握したい方- 事業化視点で見る -

ページトップへ

おススメ サービス

おススメ astavisionコンテンツ

新着 最近 公開された関連が強い技術

この 技術と関連性が強い人物

関連性が強い人物一覧

この 技術と関連する社会課題

関連する挑戦したい社会課題一覧

この 技術と関連する公募課題

関連する公募課題一覧

astavision 新着記事

サイト情報について

本サービスは、国が公開している情報(公開特許公報、特許整理標準化データ等)を元に構成されています。出典元のデータには一部間違いやノイズがあり、情報の正確さについては保証致しかねます。また一時的に、各データの収録範囲や更新周期によって、一部の情報が正しく表示されないことがございます。当サイトの情報を元にした諸問題、不利益等について当方は何ら責任を負いかねることを予めご承知おきのほど宜しくお願い申し上げます。

主たる情報の出典

特許情報…特許整理標準化データ(XML編)、公開特許公報、特許公報、審決公報、Patent Map Guidance System データ