【０００１】
【発明の属する技術分野】
本発明は、有限要素法の高性能化に関し、無節点（Ｎｏｄｅｌｅｓｓ）有限要素法モデルによる有限要素法の高性能化に関する。
【０００２】
【従来の技術】
有限要素法（ＦＥＭ（Ｆｉｎｉｔｅ Ｅｌｅｍｅｎｔ Ｍｅｔｈｏｄ））は構造物の変位や応力の解析や、電気或いは流体等の非構造物の特性解析に用いられ、解析計算装置としてのコンピュータ技術の発達と相俟って近年の応用力学工学分野、より具体的にはダムや堤防の構築やトンネル掘削、高速道路建設等の土木技術分野や、高層ビル建築、橋梁の架橋、溶鉱炉設計といった建築分野、圧延或いはプレスといった工業技術分野、電気回路設計や石油化学や食品プラントにおける流路設計等、・・・といったさまざまな産業分野等のＣＡＥ（Ｃｏｍｐｕｔｅｒ Ａｓｓｉｓｔｅｄ Ｅｎｇｉｎｅｅｒｉｎｇ（計算機援用エンジニアリング））の分野で実用性のある解析手段技術として多大な貢献をなしている。
【０００３】
ここで、従来の有限要素法（以下、節点有限要素法（ＮＯＤ ＦＥＭ）と記す）は解析すべき点を有限個の小領域（要素）に分割し、要素同士は節点で力学的に連結しており、要素の集合は連続領域を構成するものと仮定して解析を行うようにしている。
【０００４】
換言すれば、従来技術である節点有限要素法ではまず図１９に示すように連続体３００の内部に有限の離散的な接点３００−１、３００−２、・・、３００−６を定義する。また、いくつかの接点で囲まれた小領域３０１、３０２、・・、３０５を要素と呼ぶ。そして要素の力学的特性は仮想仕事の原理などの変分原理を用いて要素自身を囲む節点値（未知数）により数式表現する。これはつぎのマトリックス方程式となる。
【０００５】
【数１】
なお、図１９は膜のような平面領域問題のイメージを示す。また、図２０は節点に作用する変位及び集中力の説明図である。従来技術である節点有限要素法では図１９の例では小領域の膜（この例では三角形の小領域３０１、３０２、・・３０５）は節点間を結ぶバネに置き換わる。
【０００６】
図１９に示すような平面領域問題では、式１の左辺の｛ｕ｝は要素（小領域）を構成する節点の変位を表すベクトルであり、右辺の｛ｆ｝は節点に作用する集中力をあらわすベクトルであり、左辺の係数［Ｋ］ は複数のバネの働きを表すマトリックスである。図１９からも明らかなように１つの節点は複数の要素に共有されている。また、節点値は未知数であるが複数要素に対して共通に扱うものとする。また、１つの節点に及ぼす力は複数の節点からの寄与がある。
【０００７】
平面応力を求める場合、図２０に示すように節点ベクトルｆｉと集中共通ベクトル［Ｐｘ，Ｐｙ］Ｔ が節点３１に作用するものとする。この場合、節点３２０における値からの釣り合い方程式（平衡方程式）は次式となる。
【０００８】
【数２】
上記式２に式１を代入することにより節点力∫ｉｘ，∫ｉｙを消去する。このとき式２は節点を共有するすべての要素からの寄与がある。式２のΣｆｉは要素剛性（バネ）［Ｋ］の重ね合わせを意味する。式２の重ねあわせを全領域の節点に関して実行すると次式を得る。
【０００９】
【数３】
ここに、
［Ｋ］：全体剛性マトリックス、｛Ｕ｝：全体節点変位ベクトル
｛Ｐ｝：全体荷重ベクトル
である。
【００１０】
図２１は直列バネ問題（荷重Ｐを与えたときの直列バネの各バネ部分の変位の解析）の一実施例を示す図であり、直列バネ２１０は直列に接続するバネ２１１、２１２、２１３からなる。また、図２１で符号２１−１〜２１−４は節点有限要素法でいう節点であり、Ｐは荷重、Ｋ１〜Ｋ３はバネ２１１〜２１３の定数（バネ定数）を示す。また、図２１の直列バネ問題について式３に相当する式は次式となる。
【００１１】
【数４】
【００１２】
ここで、上記式４は全体節点変位ベクトル｛Ｕ｝を未知数とする連立一次方程式であり、マトリックス［ｋ］の２番目及び３番目の対角項がバネ定数の和（ｋ１＋ｋ２、ｋ２＋ｋ３）になっている。このような剛性要素の重ね合わせ（＝和）が従来技術である節点有限要素法、つまり、節点のある有限要素法の特徴である。
【００１３】
要約すれば、従来技術である節点有限要素法の基本式は、節点における平衡不定式である式２であり、連立一次方程式（式４）の各行は節点に関する平衡方程式である。
【００１４】
なお、荷重は、通常、物体の表面あるいは内部に分布して加わるが、節点有限要素法では分布荷重を力学的に等価な集中荷重に変換して節点に作用するものとして扱わなければならず、このため（つまり、荷重処理のために）力学的境界条件という概念を導入する。
【００１５】
図２２は力学的境界条件及び幾何学的境界の説明図であり、図２２で符号Ｐは荷重ベクトル、符号Ｓσは力学的境界である。また、符号ＳＵの表面（斜線で示されている部分）のように変位が指定されている部分（変位０を含む）を幾何学的境界（＝変位境界）という。ＳＵの節点値の指定は、変位境界条件を上記式４の左辺の｛Ｕ｝に代入することにより連立一次方程式（図４のマトリックス［ｋ］）の数値計算が実行可能となる。つまり、荷重ベクトルＰ１、Ｐ４は既知データであるから変位境界条件を変位ベクトル｛Ｕ｝に代入すれば、連立一次方程式解析プログラムを用いてコンピュータで数値計算を行い、結果を得ることができる。
【００１６】
以上により、全体節点変位ベクトル｛Ｕ｝を求め、このベクトルの中から特定の要素の節点変位成分｛ｕ｝を取り出して式１に代入すれば節点力ベクトル｛ｆ｝を求めることができる。ここで、節点力ベクトル｛ｆ｝はバネの内力に相当する。また、内力を変換した物理量が応力である。
【００１７】
【発明が解決しようとする課題】
通常構造物の強度設計においては、応力が許容値以下の値になるように構造物の形状あるいは材質を変更する作業が伴うが、従来技術である節点有限要素法では２次元または３次元の要素構成マトリックスを作る際には数値積分を実行する必要がある。
【００１８】
しかしながら、一般に数値積分には計算時間がかかるという問題点と、応力分布が要素と要素の境界（例えば図１９のＬ３１３６で示す部分）が不連続になるといった問題点があり、精度のよい結果を得るためには大規模解析計算用のポストプロセッサが不可欠になるといった問題点があった。
【００１９】
また、前述したように節点有限要素法では分布荷重を力学的に等価な集中荷重に変換して節点に作用するものとして扱うので、図３に示すような要素の辺の中間に他の要素のコーナー（ｃｏｒｎｅｒ）が配置されるイレギュラーなメッシュ分割については解析計算を実行することができないといった問題点があった。このことは、例えば、原子炉内部の熱交換用のパイプの歪みや破壊シミュレーションが困難であることを意味している。
【００２０】
つまり、節点有限要素法は節点法に基いているので要素の破壊により新たな節点が生じた場合に説店の平行方程式が新たに生じることとなる。しかし、破壊することによって生じた節点が隣接要素の返上に存在してはならないという点と、要素剛性を重ね合わせなければならないという点がコンピュータにより計算を行う際のコーディングの障壁となるという問題点がある。このことは従来技術である節点有限要素法では要素の破壊現象の追跡やシミュレーションが極めて困難であり、実用性がないということを意味している。
【００２１】
本発明は上記節点有限要素法の問題点を解決するために創案されたものであり、解析対象領域を有限個の要素に分割するが、節点を用いることなく高精度の解析を行い得ると共に、節点有限要素法では困難であった応用力学の諸問題にもコンピュータによる数値解析の手段を与える無節点有限要素法（ＮＯＤＥＬＥＳ ＦＥＭ）の提供を目的とする。また、本発明は大規模解析計算用のポストプロセッサでなくても、例えば、中規模のコンピュータ、或いはパソコン等でも短時間に高精度の解析を行い得る手法の提供を目的とする。
【００２２】
【課題を解決するための手段】
上記課題を解決するために、第１の発明の無節点有限要素法は、負荷の対象領域を小領域に分割し、分割された各小領域について負荷により生じた変位を計算し、計算により取得した各小領域の変位を隣接する各小領域の内部境界又は外部境界との相関関係により座標変換し、座標変換された各小領域の変位を基に対象領域全体の変位を算出することを特徴とする。
【００２３】
また、第２の発明の無節点有限要素法は、負荷の対象領域をサイズの異なる構造格子を重ねるようにして小領域に分割し、分割された各小領域について負荷により生じた変位を計算し、計算により取得した各小領域の変位を隣接する各小領域との内部境界又は外部境界との相関関係により座標変換し、座標変換された各小領域の変位を基に対象領域全体の変位を算出することを特徴とする。
【００２４】
また、第３の発明の無節点有限要素法は、上記第１又は第２の発明の無節点有限要素法により算出された変位を所定の変位・応力変換式により変換して各小領域又は対象領域全体の応力を得ることを特徴とする。
【００２５】
また、第４の発明の無節点有限要素法は、上記第１又は第２の発明の無節点有限要素法により算出された変位を所定の変位・境界力変換式により変換して分割された対象領域内の内部境界力及び外部境界力を得ることを特徴とする。
【００２６】
また、第５の発明は上記第１乃至第３のいずれかの無節点有限要素法において、対象領域は２次元又は３次元の固体上に画定された領域であることを特徴とする。
【００２７】
【発明の実施の形態】
各実施例の説明に先立ち、まず、無節点有限要素法により要素間の境界辺上の変位及び応力を節点を用いずに求める際に用いる各数式の導出方法について説明する。
要素間の境界辺上の変位及び応力を節点を用いずに求めるため、無節点有限要素法では次式１のエネルギー原理式で表される無節点有限要素の概念を導入する。
【００２８】
【数５】
ここに、ｕｉ：変位成分、ｔｉ：境界力成分、σｉｊ：応力成分、
バーｕｉ：規定変位 バーｔｉ：表面荷重成分又は隣接要素から外力として作用する境界力
バーｐｉ：物体力成分
である。
【００２９】
より具体化するため、平面応力について図２の三角形要素を例として無節点有限要素法による定式化について説明する。
図１は無節点有限要素法における無節点有限要素（小領域）の一実施例を示す図であり、符号１０は無節点有限要素としての三角形要素、符号１１〜１３は隣接要素（無節点有限要素）の例である。また、符号ｓ１、ｓ２、ｓ３は無節点有限要素１０の辺、符号ｓ３、ｓ４、ｓ５は隣接要素１１の辺、符号ｓ６、ｓ７、ｓ１は隣接要素１２の辺、符号ｓ８、ｓ９、ｓ２は隣接要素１３の辺であり、無節点有限要素１０と隣接要素１１〜１３が共有する辺は内部境界、辺ｓ４〜ｓ９は外部境界である。この場合、図１の三角形要素１０の内部の点Ｐの変位ｕを次式（式６）のように有限べき級数に展開する。
【００３０】
【数６】
ここに、ａｉ、ｂｊ ・・・は未知数である。
【００３１】
式６を座標と未知数に分離し、マトリックス−ベクトル形式で表すと次式（式７）のようになる。
【００３２】
【数７】
【００３３】
上記式７からひずみεが与えられる。また、ひずみ−応力関係式（σ＝Ｄ（ε−ε０）＋σ０；Ｄは材料定数を含んだ弾性マトリックス）から応力σが与えられる。また、三角形要素１０の点Ｐの応力σは次式８に示すひずみ−応力関係で与えられ、境界ｓ１〜ｓ３の境界力ｔｉはＣａｕｓｈｙ（人名）の式（次式８）から与えられる。
【００３４】
【数８】
上記式７、８、９を式５に代入して数学的操作を施すと次式（式１０）を得ることができる。
【００３５】
【数９】
【００３６】
図２は４要素に分割された平面板の例である。この例では要素１１の左辺で変位ｕ，ｖが固定され、要素２４の右辺にｘ方向の分布荷重バーｔが分布する。この例に上記式１０を摘要すると次式（式１１）が得られる。なお、要素２１〜２４は無節点有限要素である。
【００３７】
【数１０】
ここに、Ｍ１１、Ｍ１２はマトリックス、ａｉは要素１１の未定係数ベクトルである。
ここで、式１１をマトリックス−ベクトル表示すると次式（式１２）になる。
【００３８】
【数１１】
【００３９】
この式１２が計算すべき最終的な連立一次方程式である。環境条件は導入済みであるからコンピュータでマトリックス演算を行って連立一次方程式の係数値｛ａ１、ａ２、ａ３、ａ４ ｝を算出し、要素２４の右辺に分布するｘ方向の分布荷重バーｔを得ることができる。
【００４０】
図３は要素の辺上に別の要素のコーナーがあるイレギュラーメッシュの一例を示す図である。
従来技術である節点有限要素法の基本式は節点における平衡方程式である式２であり、連立一次方程式（式４）の各行は節点に関する平衡方程式である。図３に示したようなイレギュラーメッシュは塑性変形のシミュレーションの際に行うメッシュ分割で望ましいものであり、他の要素の辺の中間に他の要素のコーナーを配置するようにして分割するわけであるが、従来技術である節点有限要素法では節点における平衡方程式を順次計算し、その結果を次の要素のコーナーの各節点における平衡方程式を計算するといった操作を繰り返す必要があり、マトリックス演算の繰り返しが必要となり、計算時間が掛かるので大規模ポストプロセッサを用いる必要がある。また、計算結果の精度がマトリックス演算を繰り返す度に低下し、最終計算結果の信頼性が保証されないこととなる。
【００４１】
これに対し、本発明の無節点有限要素法は全領域を有限個の要素に分割するが、節点を用いないので前述した式２の平衡方程式を立てる必要がない（つまり、図２０の例では節点３２０に節点ベクトルｆ）と集中共通ベクトル［Ｐｘ，Ｐｙ］Ｔ が作用したが、図１、図２の例では節点を用いないので平衡方程式を立てる必要がない）。従って、図３に示したようなイレギュラーメッシュについても節点有限要素法のようにマトリックス演算を何度も繰り返す必要がないので大規模なポストプロセッサを必要とせず、また、計算結果の精度の信頼性を保証できる。また、要素毎に未定係数ベクトルａｉの数が異なっていても計算を実行することができる。このことは、無節点有限要素法ではメッシュ分割のプリプロセッサの構成が容易であること及び応力集中に関する計算結果の精度がよいことを意味する。
【００４２】
また、従来技術である節点有限要素法は離散集中系であるが、本発明の無節点有限要素法は全領域を有限個の要素に離散化するが力学的には分布形であることから、２次以上の変位関数を用いた無節点有限要素法に基く解析計算による要素間の応力の分布は図４に示すように連続的となる。
【００４３】
また、節点有限要素法における全体節点変位ベクトルを表す連立一次方程式（式４）と、無節点有限要素法における応力分布を表す連立一次方程式（式１２）の対比から無節点有限要素法には全体マトリックスを作るときに重ね合わせの過程（プロセス）がない。つまり、無節点有限要素法では式１２の連立一次方程式の順序は並び替えても計算結果に影響しない。
【００４４】
＜実施例＞
１．ハードウエア構成
図５は無節点有限要素法による解析計算を実行可能な数値解析システムの一実施例の構成を示すブロック図であり、この例は数値解析システムとしてパーソナルコンピュータ（以下、パソコン）を例としている。
図５で、パソコン５は制御部５１、マウス等のポインティングデバイス５２、入力装置５３、作業用メモリ５４、保存記憶用メモリ５５、表示装置５６、プリンタ５７を有している。
制御部５１はＣＰＵ（又はＭＰＵ）およびプログラム格納用ＲＯＭのほか周辺制御回路を備えた構成をなし、装置全体の動作を制御すると共に、保存記憶用メモリ２３に保存記憶されている無節点解析計算プログラム群（無節点解析手段）により本発明の無節点有限要素法に基く無節点解析処理の実行を制御する。
【００４５】
ポインティングデバイス５２は解析担当者による領域指定時の位置や描画パラメータの指定や各種操作指示等の手段であり、操作がなされると位置情報を制御部５１に送出する。また、入力装置５３はキーボードやマウス等のポインティングデバイス等を有しており、ユーザはデータやパラメータを入力したりすることができる。
【００４６】
また、作業用メモリ５４はＲＡＭ等の一時記憶メモリからなり、パソコンの起動時にプログラム格納用ＲＯＭに記憶されている制御プログラムが読み出されて記憶されると共に、固体力学の解析処理が選択されると具体的なモデル毎に対応の解析用プログラムが保存記憶メモリ５５から適時読み出されて記憶される。
【００４７】
また、保存記憶用メモリ５５は磁気ディスクやＦＤ等のリムーバブルな記録媒体からなっており、本発明の無節点有限要素解析手順を実現するために作成されている無節点有限要素解析用プログラムモジュール群及び解析支援プログラムや、入力データ、パラメータおよび解析結果（出力データ）、精度評価用データのほか、描画プログラムを含む対象領域設定プログラムや、領域分割プログラム、中抜き演算によるマトリックス方程式計算プログラム等のサブプログラム群やサポートプログラム群を保存記録している。
【００４８】
２．モデル
以下、上記パソコン５を用いた無節点有限要素法に基く数値解析の一実施例として従来の有限要素法（節点有限要素法）では計算が困難であった応力集中モデルを取り上げるが、無節点有限要素法の適用範囲は応力集中モデルに限定されない。つまり、従来の有限要素法で解決可能な力学モデルについては高速でより精度が高い計算結果を得ることが可能であり、それに加えて応力集中モデルのように従来困難であった分野のアプリケーションについても高速で且つ精度の高い結果を提供できる。
【００４９】
剛体変位の一例としての応力集中モデル（応力がどの部分にどのように集中し、伝播するかを考慮すること）は、例えば、原子炉のように冷却水還流用のパイプ孔を有する密閉容器、プラントの隔壁、ネジ止め溶接された板からなる構造物（例えば、航空機、船、橋等々）の設計時に必要となるものであり、図７に示すような有孔板の引張荷重モデルをその典型的な例としてあげることができる。
【００５０】
３．無節点有限要素法による手順
図６は無節点有限要素法によるモデル解析手順を示すフローチャートである。
まず、解析対象とする材料や構造物等（以下、対象領域と記す）を特定すると共に解析対象とする部位（以下、対象領域と記す）と対象領域中の特異領域（例えば、平板に穿たれた穴の周辺、屈曲面の周辺、溶接個所及びその周辺等）を特定し、そのサイズや重量等、特性、荷重等の負荷量、値からの方向等の条件値（パラメータとして変化させる場合は、そのピッチも含む）を設定する。なお、対象領域は通常は２次元（平面）又は３次元（立体）であるがこれに限定されない（ステップＳ１）。
【００５１】
次に、対象領域（特異領域がある場合はそれも含む）を図７（ａ）の例に示すようにモデル化する。モデル化にあたってはモデル全体を１とした正規化を行うことが望ましく、更に、図７（ｂ）の例に示すように解析計算が効率よくできるように解析対象領域（図７（ｂ）の例で符号７１’で示す部分）を決定することが計算上望ましい。モデルの作成及び正規化は後述するようにパソコン５を用いて容易に行うことができる（ステップＳ２）。
【００５２】
次に、モデルの解析対象領域全体をｎ×ｎ（ｎは２以上の整数）のメッシュに分割し（ステップＳ３）、更に、上記ステップＳ１で特定した特異領域の周辺についてステップＳ３で分割した各メッシュを細分して小メッシュに分割し、各小メッシュ、つまり各要素に番号をつけると共に各辺の長さ（＝分割数でもよい）を対応付ける。この際、図８に示すように対象領域全体から特異領域を絞り込むように何段階かに分けてメッシュの細分を行うようにするとよい。このようなメッシュ分割は後述するようにパソコン５を用いて簡単な操作で行うことができる（ステップＳ４）。
【００５３】
まず、第１段階（つまり、ステップＳ３）でメッシュ（２次元の場合は正方格子、３次元の場合は立方格子）に分割した各小領域（要素）毎にその座標値及び材料等の条件値（大きさ、重量、特性値、負荷、力の方向等（モデルが正規化されている場合は条件値も正規化したもの））を与え、前述した式１２の連立一次方程式を演算して各小領域毎の変位（未定係数ベクトル）を算出する。つまり、無節点有限要素法では小領域毎に結果が算出される（ステップＳ５）。
【００５４】
次に、第１段階で得たメッシュのうち細分されたメッシュについて段階ごとに細分された各小領域について前述した式１２の多元連立一次方程式（正方マトリックス方程式）を演算して各領域毎の変位（未定係数ベクトル）を算出する。この場合、各小領域に与えられる変位の総和は上記ステップＳ５で算出した小領域（細分前の小領域）の変位である。更に小領域が細分されている場合は同様の操作を繰り返す。つまり、サイズの異なるメッシュを重ね合わせることと同じである。なお、無節点有限要素法では各小領域を独立的に扱えるので、本願出願人が１９９９年５月２５日に出願の発明（マトリックスの中抜き圧縮処理（特開平２０００−３３９２９５））に基くマトリックスの中抜き圧縮処理を行ないつつ連立一次方程式の計算を行い、未定係数ベクトル（図７の例では各小領域ごとの変位）の近似解を得て、それを基に解析対象領域の変位の近似解を得ることができる（ステップＳ６）。
【００５５】
次に近似精度の検証を行い、精度が許容誤差内になった場合は解析計算を終了する。また、許容誤差に至らない場合はステップＳ３に戻って各要素を更に分割するか、ステップＳ１に戻って条件値の設定を変更する。なお、本発明における近似精度の検証は本願出願人によって同日に出願された「無節点有限要素法における近似解の収束手段」を用いて行う（ステップＳ７）。
【００５６】
上記ステップＳ７で得た変位を基にひずみ−応力関係式（式８）から応力を得ることができる。また、Ｃａｕｓｈｙ（人名）の式から境界力を得ることができる。
【００５７】
４．応力集中モデル
図７は、中央に円孔を有する正方形板の両端に一様な引張荷重を与えた有孔板の引張荷重モデルの説明図であり、図７（ａ）で符号７１は正方形の板、符号７２は板７１の中央に開口した円孔、符号７３、７３’は板の両端に一様に加えられた引張荷重Ｓを意味し、引張荷重７３、７３’はそれぞれ反対方向に加えられる。また、図７（ｂ）は板７１が二軸対称であることから解析対象領域をその１／４の板７１’として取り扱うことを示す図である。なお、この例では、板の長さを２Ｌとしたとき、Ｌ＝５０ｍｍ（ミリ）、孔の半径ｒ＝１０ｍｍ、引張荷重Ｓ＝１００ｋｇｆ／ｍｍ、Ｅ＝２００００ｋｇｆ／ｍｍ、ν＝０．３とした。また、符号７４〜７５は解析対象領域７１’の左側の隣接領域との境界辺（ＡＢ）で固定される変位、符号７６〜７９は解析対象領域７１’の下側の隣接領域との境界辺（ＣＤ）ので固定される変位を表す。なお、変位７５、７６は円孔７２との境界辺（弧ＢＣ（１／４周））に働く変位でもある。
【００５８】
図８は解析対象領域（図７（ｂ）の１／４板７１’）を図６のフローチャートのステップＳ３、Ｓ４で述べたように領域分割した例を示し、図８（ａ）はメッシュ（構造格子）による分割を重ねて小領域（要素）を得た例（つまり、１／４板７１’を一番大きな構造格子で分割した１６個の小領域８１を得てから、そのうちの円孔に近い小領域の分割を１／２の縮尺で更に３回繰り返して、円孔（内部境界領域）７２との境界に近い小領域８２、８３の分割を順次進めるようにして計４回の分割で、１６個の小領域８１、２０個の小領域８２、２７個の小領域８３及び７個小領域８４、つまり、合計７３個の小領域（要素）を得た例である。
【００５９】
また、図８（ｂ）は更に円孔７２の近傍の小領域８３、８４（円孔との境界領域を含む）を３分割して１３７個の小領域（要素）を得た例である。
【００６０】
このような場合、従来の有限要素法（節点有限要素法）の分割方法では、最初１６分割し、更に１／２分割を３回繰り返すと１２８要素となり、しかも、再分割を行う毎に全ての要素について計算を行わなければならないので大規模計算機による演算を繰り返さなければならず誤差が蓄積して精度がおちるが、本発明の分割方法では図８（ａ）に示したように７３要素で済み、しかも再分割した範囲の要素の解析だけでよいので、計算が簡単になり、計算結果（近似解）の精度がよい。また，従来の有限要素法の分割方法では、最初１６分割し、更に２分割を３回繰り返してから３分割を１回行うと４８４要素となり、再分割を行う毎に全ての要素について計算を行わなければならないが、本発明の分割方法では図８（ｂ）に示したように１３７要素で済み、しかも再分割した範囲の要素について計算を行うだけでよい。
【００６１】
３．パソコンによる計算手順
図９は無節点有限要素解析システムの起動時に表示されるジョブ選択メニューの一実施例を示す図であり、図１０は図５に示したパソコンで図７（ａ）の有孔板の引張荷重モデルに働く変位及び応力の解析計算を行なう場合の制御部５１による制御動作例を示すフローチャートである。以下の説明で制御部５１は無節点有限要素解析システムを構成する制御プログラム及び解析用プログラムに従って制御及び演算を行う。
【００６２】
（動作例）
制御プログラム及び各種解析用プログラム群及び画面表示データ等からなる無節点有限要素解析システムがインストールされているパソコン５を起動すると、制御部５１は表示部５６の画面上に図９に示すようなジョブ選択メニューを表示する（ステップＴ１）。解析担当者が表示されているジョブ選択メニューの中から２次元無節点有限要素解析アプリケーションを選択すると（ステップＴ２）、制御部５１は描画プログラムを含む対象領域設定プログラムを保存記憶用メモリ５５から読み出して作業用メモリ５４のアプリケーションプログラム領域に記憶すると共にその対象領域設定プログラムを起動し、解析担当者に対象領域（平面）の描画を促す（ステップＴ３）。
【００６３】
解析担当者がマウス５２で表示装置５６の画面上に対象領域（図７の例では有効板）を描画し、サイズ、及び単位等をキー入力すると制御部５１は描画された対象領域を図形化（図７（ａ）参照）して表示し、解析担当者に特異領域と荷重等の負荷が働く部分及びその方向等の指定を促す（ステップＴ４）。
【００６４】
解析担当者が図形化された対象領域上の特異領域（図７の例では円孔７２）と負荷が働く部分又はその方向（図７の例では右辺及び左辺に働く引張荷重７３、７３’）等をマウス５２で指定すると（ステップＴ５）、制御部５１は解析対象領域図形化プログラムにより図形の対称性等に基いた解析対象領域（図７（ｂ）参照）を抽出して表示装置５６の画面に表示する。解析対象領域の表示はその部分を破線等で囲んだり、対象領域の表示色とは異なる色で表示したりして明示する。また、ウインドウ表示するようにしてもよい。また、解析担当者の指示により解析対象領域の修正や変更、或いは代替案の表示を行うこともできる。なお、対象領域内に特異領域が複数ある場合は、解析担当者はそれぞれ複数個の特異領域を指定することができる。また、複数個の特異領域が指定された場合は、特異領域を１つづつ含むように特異領域の数分の解析対象領域を表示し、解析担当者に解析計算の順位指定を促す。また、解析担当者は、この際、図形の辺や点又は面にマウスカーソルを固定してキーボードから文字或いは記号を入力して図形の辺や点又は面に名称或いは記号を付することができる（ステップＴ６）。
【００６５】
次に、制御部５１は、領域分割プログラムを保存記憶用メモリ５５から読み出して作業用メモリ５４のアプリケーションプログラム領域に上書き記憶してその領域分割プログラムを起動し、ｋ＝１として、解析担当者によって指定された解析対象領域の第ｋ次領域分割数ｎの入力を促す（ステップＴ７）。
【００６６】
解析担当者が第ｋ次領域分割数ｎ（図７の例では５）を入力すると、制御部５１は領域分割プログラムにより解析対象領域を図８に示したようにｎ×ｎ（図７の例では５）のメッシュに分割し、画面上に再分割要否選択ボタン（又は選択メニュー）等を表示して分割された小領域（要素）の再分割要否の選択を促す（ステップＴ８）。
【００６７】
解析担当者が「再分割否」を選択した場合はステップＴ１２に遷移し、「再分割要」を選択した場合は、制御部５１はステップＴ７でｋ次領域分割された解析対象領域のうち更に分割対象とする領域範囲の指定を促すので（ステップＴ９）、解析担当者が特異領域を含む範囲（１次領域分割されたメッシュのうちの特定の範囲）をマウスで指定し、更に、分割数ｎｋ＋１を入力すると、第ｋ＋１次分割が行われ、指定範囲内の各メッシュを（ｎｋ＋１ ）×（ｎｋ＋１ ）分割し、更に、画面上に再分割要否選択ボタン（又は選択メニュー）等を表示して分割された小領域（要素）の更なる再分割要否の選択を促す。つまり、ｋ＝１の場合は、図７（ａ）の例で１１番目のメッシュ（小領域）８２−１１と２３番目のメッシュ８２−２３を指定し、分割数２を入力すると、特異領域（孔７２）側の９個のメッシュ８２が更に２×２の３６の小領域８３に分割され、更に、孔７２側の小領域８３の再々分割を行うか否かの選択要求が行われる（ステップＴ１０）。
【００６８】
更に、解析担当者が「再々分割否」を選択した場合はステップＴ１２に遷移し、「再々分割要」を選択した場合はｋに１を加えてステップＴ１０に戻る（ステップＴ１１）。
【００６９】
図８に示したような無節点有限要素法による領域分割が終了すると、制御部５１は解析支援プログラムを読み出して作業用メモリ５４のアプリケーションプログラム領域に上書き記憶してその解析支援プログラムを起動し、分割された各小領域に各要素番号を付与し、表示する。要素番号は通常自動付与されるが、解析担当者が小領域を指定して番号を付与することもできる（ステップＴ１２）。
【００７０】
要素番号の付与が終わると、制御部５１は表示装置５６に案内メッセージを表示して解析担当者にデータの入力を促すので、解析担当者が分割した各小領域（サイズは異なっていてもよい）の変位を未定係数ベクトルとする多元連立一次方程式（式１２参照）をたて（ステップＴ１３）、初期値（図７、図８の例では、サイズ、右辺の各要素に働く引張荷重Ｓ）を導入すると、制御部５１は特開平２０００−３３９２９５に開示のマトリックスの中抜き圧縮処理を行いつつ未定係数の近似値を求め、求めた未定係数の近似値から解析対象領域の変位を算出する（ステップＴ１４）。
次に、算出した近似値が所定の誤差範囲（許容誤差）内か否かを判定し、許容誤差内の場合にはステップＴ１９に遷移し、そうでない場合にはＴ１６に遷移する（ステップＴ１５）。
【００７１】
制御部５１は、対象とする要素（小領域）及びその周辺要素をパソコン５の表示装置２４の画面に表示し、更に、ステップＴ８で用いたものと同じ領域分割プログラムを保存記憶用メモリ５５から読み出して作業用メモリ５４のアプリケーションプログラム領域に上書き記憶してその領域分割プログラムを起動し、解析担当者によって対象とする要素（小領域）の分割数ｊの入力を促す（ステップＴ１６）。
【００７２】
解析担当者が分割数ｊを入力すると、制御部５１は領域分割プログラムにより対象とする要素をメッシュで分割する（ステップＴ１７）。
【００７３】
対象とする要素の分割が終了すると、制御部５１は解析支援プログラムを読み出して作業用メモリ５４のアプリケーションプログラム領域に上書き記憶してその解析支援プログラムを起動し、分割された各小領域に各要素番号を付け直してステップＴ１４に戻る。これにより、再分割された要素を加えた要素全体について更に詳しい変位の近似解を得ることができる（ステップＴ１８）。
【００７４】
計算結果（近似解）が許容誤差内の場合はそれを解として表示装置５３に表示する。また、指定によりプリンタ５７で印刷出力したり保存記憶用メモリ５５に書き込むこともできる（ステップＴ１９）。
【００７５】
また、上記ステップＴ１８で得た変位の近似解を基にひずみ−応力関係式（式８）から応力を得ることができる。また、Ｃａｕｓｈｙ（人名）の式から境界力を得ることができる。
【００７６】
図１１は無節点有限要素法による有孔板の引張荷重モデルのＹ−Ｙ’間の変位の分布を示す図であり、上記図１０のフローチャートに基いてパソコン５で得た結果を図示したものであり、図７（ａ）のＹ−Ｙ’間、つまり、図７（ｂ）の１／４板７１’でＹ軸方向に１０ｍｍ（有孔部分の上部）から５０ｍｍ（端部）の間に生ずる応力の分布を示したものである。
【００７７】
図１１から明らかなように、Ｔｒｅｆｆｔｚ（人名）による計算結果と無節点有限要素法による計算結果は一致している。つまり、図７の有孔板モデルは図８に示すように領域の重ね合わせ分割によるイレギュラーメッシュであり、従来の節点有限要素法では計算時間及び精度の点から不可能に等しかったイレギュラーメッシュについて計算が容易に、しかも高精度でできることを示している。図１１で孔の傍（１０ｍｍ位置）の変位が最も大きく、距離が大きくなるにつれて変位が少なくなる様子が明らかとなっている。
【００７８】
図１２は図７（ｂ）のように領域分割した１／４板７１’に引張荷重Ｓを加えた場合の板７１’の変形の様子を示す図、すなわち、無節点有限要素法による有孔板の引張荷重モデルのＹ−Ｙ’間の変位の分布のシミュレーションの一実施例を示す図であり、図１２から円孔７２の近傍の変形が著しいことがわかる。シミュレーションは図１０のフローチャートのステップＴ２でジョブ選択メニューのうちの２次元解析シミュレーション又は３次元解析シミュレーションを選択し、変数（パラメータ）の初期値及びピッチを指定して図１０のステップＴ１４〜Ｔ１８に相当する動作を繰り返すことにより行うことができる。
【００７９】
なお、上記ステップＴ３で描画プログラムを用いたがＣＡＤプログラムでもよい。また、上記ステップＴ８、Ｔ１０で領域をメッシュ分割する際のメッシュの形状は図８では正方格子の例を示したが正方格子に限定されない。
【００８０】
また、上記ステップＴ２で２次元無節点有限要素解析アプリケーションを選択した場合について述べたが３次元無節点有限要素解析の場合も、ステップＴ３以下の動作と同様の動作で近似解を得ることができる。この場合、解析対象領域はステップＴ８でｎ×ｎ×ｎの立体格子（正方立体格子に限定されない）で分割される。また、必要に応じてステップＴ１０で再分割できる。
【００８１】
なお、上記例ではパソコンを用いた例について説明したが本発明による解析計算装置はパソコンに限定されない。つまり、中型、大型、もしくは大規模電子計算機についても適用できることはいうまでもない。即ち、解決を要する固体力学問題の規模に応じて解法モデルを作り適宜コンピュータの規模を選べばよい。以下、解析装置の他の実施例について述べる。
【００８２】
図１３は、本発明の無節点有限要素法による計算を実行可能な解析計算装置の一実施例の構成を示すブロック図であり、多重プロセッサ（ｍｕｌｔｉｐｒｏｓｅｓｓｏｒ）による並列処理を行なう例である。図１３で、多重プロセッサ１３０は複数の処理装置（プロセッサ）１３１−１、１３１−２、・・・、１３１−ｎ（実施例ではｎ＝１６）が記憶装置１３２と入出力装置１３３を共有している。また、入出力装置１３３には外部装置（この例ではパーソナルコンピュータ（以下、パソコン）等の電子計算機）１３５が接続している。なお、入出力装置１３３はキーボード等の入力装置、プリンタやモニター装置等の出力装置および磁気ディスクや光ディスク装置からなっている。
【００８３】
上記ハードウエア構成で、外部装置１３５側では図１０のフローチャートのステップＴ１〜Ｔ１３、Ｔ１５〜Ｔ１９を実行し、多重プロセッサ１３０側ではステップＴ１４を実行する。なお、図１３の例では多重プロセッサ１３０と外部装置１３５をケーブル等で接続した構成を示したが、多重プロセッサ１３５と外部装置１３５は分離していてもよく、データの授受を記録媒体を介して行なうようにしてもよい。
【００８４】
図１４は本発明の無節点有限要素法による計算を実行可能な解析計算装置の他の実施例の構成を示すブロック図であり、通信ネットワークを用いた並列的処理によって解析計算を行なう場合の構成例を示す。
【００８５】
図１４で、ネットワーク１４０は通信回線１４５に接続するコンピュータ１４２−１、１４２−２、・・・、１４２−ｎ（実施例ではｎ＝６）と、これらコンピュータ装置１４２−１、１４２−２、・・・、１４２−ｎとの間で通信回線１４５を介してデータを授受すると共に、コンピュータ１４２−１、１４２−２、・・・、１４２−ｎ間でのデータの授受を制御するサーバ１４１を有している。なお、ネットワーク１４５はインターネット等の通信ネットワークを用いてもよく、或いは各研究所等を結ぶ専用回線からなるネットワークでもよく、また、これらの混合型であってもよい。
【００８６】
サーバ１４１は処理装置、通信制御部、メモリおよび入出力装置を備えたコンピュータからなり、データベースを接続できる。
【００８７】
コンピュータ１４２−１、１４２−２、・・・、１４２−ｎは大型コンピュータに限られず、所定次数の正方マトリックス演算を行ない得る程度の処理速度とメモリー容量を有し、通信制御機能を備えた装置であればよく、例えば、パソコンやワークステーション等であってもよい。
【００８８】
なお、正方マトリックス配分の要請上、数値解析を行なう場合にはそれぞれのコンピュータの処理速度およびメモリー容量はあまり掛離れていると、最小能力のコンピュータにより効率が制限される（例えば、ネットワーク１４５に接続するコンピュータ１４２−１、１４２−２、・・・、１４２−６が大型計算機、コンピュータ１４２−７、・・・、１４２−ｎがパソコンといった構成はシステム全体の処理効率がパソコンの処理効率で限定されてしまう）ので、互いの能力差が一定範囲内にあることが望ましい。
【００８９】
通信ネットワークを用いた場合には、図１０のフローチャートに準じた手順で解析計算を行なうことができる。
この場合、ステップＴ１からＴ１２をサーバ１４１側で実行し、ステップＴ１３の初期値及びマトリックスをサーバーから各コンピュータ１４２−１、１４２−２、・・・、１４２−ｎに配分して通信回線を介して送信し、各コンピュータ１４２−１、１４２−２、・・・、１４２−ｎでステップＶ１４を実行するようにして、実行終了後サーバ１４１に近似解を送信してサーバ１４１でステップＴ１５〜Ｔ１９を実行するようにできる。
【００９０】
なお、図１４の例ではサーバ１４１と端末を別の装置としたが、各端末（コンピュータ１４２−１〜１４２−ｎ）のうちのある一つ（例えば、コンピュータ１４２−１）をサーバと兼用するようにしてもよい。
【００９１】
（剛体変位モデルへの無節点有限要素法の応用）
無節点有限要素法は下記に述べるように剛体変位を明快に扱うことができる。
例証として平面応力モデルを挙げて説明する。まず、下記式１３の変位関数を仮定する。
【００９２】
【数１２】
【００９３】
ここに、ｕ：ｘ方向の変位 ν：ｙ方向の変位
ｕ０：ｘ方向の剛体変位 ν０：ｙ方向の剛体変位
χ０：平面内の剛体回転 εｘ０、εｙ０、γｘｙ０ ：ひずみ成分
を表す。
【００９４】
式１３は要素（小領域）内の任意点の変位をひずみと剛体変位による生じる項に分離して表している。但し、ｕ０ 、ν０ 、εｘ０、εｙ０、γｘｙ０ は各要素に対する未知数であり、全体の未知数は６×要素数となる。式１３を新しエネルギー原理を表す次式１４と組み合わせてせん断荷重が作用する片持平面板のモデルを解析すると図１５、図１６を得る。
【００９５】
図１５はせん断荷重が作用する片持平面板の変形（変位）を示す図であり、図１６はせん断荷重が作用する片持平面板の応力分布を示す図である。図１５、図１６から４０×２０要素分割の場合、変位も応力も十分な精度を持っていることがわかる。
【００９６】
【数１３】
【００９７】
次に、式１３から剛体回転χ０ の項をわざわざ落とした下記変位関数（式１５）を用いて同じモデルを解析すると、変位も応力も極端に小さくなることが分かる。この数値実験により剛体変位は全て陽に取り扱うべきことがわかる。
【００９８】
【数１４】
【００９９】
次に、剛体変位を陽に含めた板曲げモーメントのたわみ関数を次式１６に示す。
【０１００】
【数１５】
【０１０１】
式１６を用いた解析計算結果を図１７、図１８に示す。図１７、図１８は等分布荷重が左右する単純支持板のたわみの計算結果を図示したものであり、メッシュ分割数は２５×２５要素である。また、図１８は数値計算するときに要素間の結合力を緩めるテクニックを使って得られた図である。
【０１０２】
本発明の無節点有限要素法では、変形後、個々の要素（小領域）がアイソレート（ｉｓｏｌａｔｅ）していることが分かる。すなわち、要素と要素がくっついていなくても計算が可能である。また、無節点有限要素法の部分対角マトリックスＭｉｊ（式１２参照）は対称で逆行列が存在する。従って、連立一次方程式のベクトル消去法又はベクトル反復法が成立する。これは、要素毎に消去することもできるし、要素毎に反復計算が可能であることも示している。このような高計算効率の性質は非線型モデルの解析計算に威力を発揮する。
【０１０３】
以上、本発明の一実施例について説明したが本発明は上記実施例に限定されるものではなく、種々の変形実施が可能であることはいうまでもない。
【０１０４】
【発明の効果】
上記説明したように、本発明の無節点有限要素法によれば対象領域を有限弧の要素に離散化するが、力学的には分布形であるため、２次以上の変位関数を用いた無節点有限要素解析による応力は要素（小領域）間の分布が連続となる。また、無節点有限要素法には全体マトリックスを作るときに節点有限要素法のように重ね合わせの過程がないので、サイズの異なる構造格子を重ねるようにして分割された小領域を含むイレギュラーメッシュについても容易に数値解析を行うことができる。
【０１０５】
また、本発明の無節点有限要素法では、変形後、個々の要素（小領域）がアイソレートしているので、すなわち、要素と要素がくっついていなくても計算が可能である。また、無節点有限要素法の部分対角マトリックスＭｉｊ（式１２参照）は対称で逆行列が存在する。従って、連立一次方程式のベクトル消去法又はベクトル反復法が成立する。これは、要素毎に消去することもできるし、要素毎に反復計算が可能であることも示している。このような高計算効率の性質は非線型モデルの解析計算に威力を発揮する。
【０１０６】
また、各小領域が隣接していなくても計算が可能であることから連立一次方程式（ベクトル方程式）の順序は並び替えても計算上何ら影響がない。この性質と重ねあわせの必要がないことを考慮してプログラムを作成すれば、小領域（要素）が力学的な意味で破壊し、全体の要素数が動的に増加しても計算が可能となる。
【図面の簡単な説明】
【図１】無節点有限要素法における無節点要素（小領域）の一実施例を示す図である。
【図２】４要素に分割された平面板の例である。
【図３】要素の辺上に別の要素のコーナーがあるイレギュラーメッシュの一例を示す図である。
【図４】無節点有限要素法に基く解析計算による要素間の応力の分布の一例を示す図である。
【図５】無節点有限要素法による解析計算を実行可能な数値解析システムの一実施例の構成を示すブロック図である。
【図６】無節点有限要素法によるモデル解析手順を示すフローチャートである。
【図７】有孔板の引張荷重モデルの説明図である。
【図８】領域分割の一実施例を示す図である。
【図９】無節点有限要素解析システムの起動時に表示されるジョブ選択メニューの一実施例を示す図である。
【図１０】パソコンで有孔板の引張荷重モデルに働く変位及び応力の解析計算を行なう場合の制御部による制御動作例を示すフローチャートである。
【図１１】無節点有限要素法による有孔板の引張荷重モデルのＹ−Ｙ’間の変位の分布を示す図である。
【図１２】無節点有限要素法による有孔板の引張荷重モデルのＹ−Ｙ’間の変位の分布のシミュレーションの一実施例を示す図である。
【図１３】本発明の無節点有限要素法による計算を実行可能な解析計算装置の一実施例の構成を示すブロック図である。
【図１４】本発明の無節点有限要素法による計算を実行可能な解析計算装置の他の実施例の構成を示すブロック図である。
【図１５】せん断荷重が作用する片持平面板の変形（変位）を示す図である。
【図１６】せん断荷重が作用する片持平面板の応力分布を示す図である。
【図１７】等分布荷重が左右する単純支持板のたわみの計算結果の一実施例を示す図である。
【図１８】等分布荷重が左右する単純支持板のたわみの計算結果の一実施例を示す図である。
【図１９】膜のような平面領域問題のイメージを示す。
【図２０】節点に作用する変位及び集中力の説明図である。
【図２１】直列バネ問題の一実施例を示す図であり。
【図２２】力学的境界条件及び幾何学的境界の説明図である。
【符号の説明】
１０、１１、１２、１３ 要素（小領域）
７１ 有孔板（対象領域）
７１’ 解析対象領域
８１、８２、８３、８４ メッシュ（小領域、要素）
ｓ１、ｓ２、ｓ３ 内部境界
ｓ４〜ｓ９ 外部境界[0001]
TECHNICAL FIELD OF THE INVENTION
The present invention relates to improving the performance of a finite element method, and more particularly to improving the performance of a finite element method using a nodeless finite element model.
[0002]
[Prior art]
The finite element method (FEM) is used for analyzing the displacement and stress of a structure and for analyzing the characteristics of a non-structure such as electricity or fluid, and is coupled with the development of computer technology as an analysis and calculation device. In recent years, applied mechanical engineering fields, more specifically, civil engineering fields such as dam and embankment construction, tunnel excavation, and highway construction, high-rise buildings, bridge construction, blast furnace design, rolling and pressing, etc. Analysis means technology that is practical in the field of CAE (Computer Assisted Engineering) in various industrial fields such as the industrial technology field, electric circuit design, flow path design in petrochemical and food plants, etc. Has made a great contribution.
[0003]
Here, the conventional finite element method (hereinafter referred to as a nodal finite element method (NOD FEM)) divides a point to be analyzed into a finite number of small regions (elements), and the elements are dynamically connected at nodes. The analysis is performed assuming that a set of elements constitutes a continuous area.
[0004]
In other words, in the nodal finite element method of the prior art, first, finite discrete contacts 300-1, 300-2,..., 300-6 are defined inside the continuum 300 as shown in FIG. Further, the small areas 301, 302,... 305 surrounded by some contact points are called elements. The mechanical characteristics of the element are expressed mathematically by node values (unknowns) surrounding the element itself using variational principles such as the principle of virtual work. This results in the following matrix equation:
[0005]
(Equation 1)
[0006]
In the planar area problem as shown in FIG. 19, {u} on the left side of Expression 1 is a vector representing displacement of a node constituting an element (small area), and {f} on the right side represents a concentration force acting on the node. The coefficient [K] on the left side is a matrix representing the function of a plurality of springs. As is clear from FIG. 19, one node is shared by a plurality of elements. Although the node value is an unknown value, it is assumed that the node value is commonly handled for a plurality of elements. Further, the force acting on one node has contributions from a plurality of nodes.
[0007]
When obtaining the plane stress, as shown in FIG. i And the central common vector [P x , P y ] T Acts on the node 31. In this case, the balance equation (balance equation) from the value at the node 320 is as follows.
[0008]
(Equation 2)
By substituting Equation 1 into Equation 2 above, the nodal force ∫ ix , ∫ iy To delete. At this time, Equation 2 has contributions from all elements sharing the node. Σf in Equation 2 i Means superposition of element rigidity (spring) [K]. When the superposition of Expression 2 is performed on the nodes in all regions, the following expression is obtained.
[0009]
[Equation 3]
here,
[K]: overall rigidity matrix, {U}: overall nodal displacement vector
{P}: Overall load vector
It is.
[0010]
FIG. 21 is a diagram showing an embodiment of the series spring problem (analysis of displacement of each spring portion of the series spring when a load P is applied). The series spring 210 is composed of springs 211, 212, and 213 connected in series. Become. In FIG. 21, reference numerals 21-1 to 21-4 denote nodes in the node finite element method, P denotes a load, and K1 to K3 denote constants (spring constants) of the springs 211 to 213. Further, the equation corresponding to equation 3 for the series spring problem in FIG. 21 is as follows.
[0011]
(Equation 4)
[0012]
Here, the above equation 4 is a simultaneous linear equation in which the total nodal displacement vector {U} is an unknown number, and the second and third diagonal terms of the matrix [k] are the sum of spring constants (k1 + k2, k2 + k3). ing. The superposition (= sum) of such rigid elements is a feature of the nodal finite element method of the prior art, that is, the finite element method with nodal points.
[0013]
In summary, the basic formula of the nodal finite element method of the prior art is Formula 2 which is a balanced indefinite formula at the node, and each line of the simultaneous linear equation (Formula 4) is a balanced equation relating to the node.
[0014]
In addition, the load is usually distributed and applied to the surface or inside of the object, but in the nodal finite element method, the distributed load must be converted into a mechanically equivalent concentrated load and treated as acting on the node, For this reason (ie, for load handling) the concept of mechanical boundary conditions is introduced.
[0015]
FIG. 22 is an explanatory diagram of a mechanical boundary condition and a geometric boundary. In FIG. 22, reference symbol P is a load vector, and reference symbol Sσ is a mechanical boundary. Also, the symbol S U A portion (including a displacement of 0) whose displacement is designated, such as the surface (portion indicated by hatching), is called a geometric boundary (= displacement boundary). S U The numerical calculation of the simultaneous linear equation (matrix [k] in FIG. 4) can be executed by substituting the displacement boundary condition into {U} on the left side of the above equation (4). In other words, since the load vectors P1 and P4 are known data, if the displacement boundary condition is substituted into the displacement vector {U}, a numerical calculation can be performed by a computer using a simultaneous linear equation analysis program to obtain a result.
[0016]
As described above, the nodal force vector {f} can be obtained by obtaining the entire nodal displacement vector {U}, extracting the nodal displacement component {u} of a specific element from this vector, and substituting it into Equation 1. Here, the nodal force vector {f} corresponds to the internal force of the spring. The physical quantity obtained by converting the internal force is the stress.
[0017]
[Problems to be solved by the invention]
Normally, in the strength design of a structure, it is necessary to change the shape or material of the structure so that the stress becomes a value equal to or less than an allowable value. However, the two-dimensional or three-dimensional Numerical integration must be performed when constructing the configuration matrix.
[0018]
However, numerical integration generally requires a long calculation time and stress distribution has a problem that the boundary between elements (for example, a portion indicated by L3136 in FIG. 19) is discontinuous. There is a problem that a post-processor for large-scale analysis calculation is indispensable to obtain.
[0019]
Further, as described above, in the nodal point finite element method, the distributed load is converted into a mechanically equivalent concentrated load and is treated as acting on the nodal point. Therefore, as shown in FIG. There has been a problem that analysis calculation cannot be performed on irregular mesh division in which corners are arranged. This means, for example, that it is difficult to simulate the distortion or fracture of the heat exchange pipe inside the reactor.
[0020]
That is, since the nodal finite element method is based on the nodal method, when a new node is generated by destruction of an element, a parallel equation of a hypothesis is newly generated. However, the problem that nodes generated by destruction must not be on the back of adjacent elements and that element rigidity must be overlapped is a barrier to coding when performing calculations by computer. There is. This means that it is extremely difficult to track and simulate the destruction phenomenon of the element by the nodal finite element method, which is the conventional technique, and it is not practical.
[0021]
The present invention has been devised to solve the problem of the node finite element method, and divides the analysis target region into a finite number of elements, but can perform high-precision analysis without using nodes, An object of the present invention is to provide a nodeless finite element method (NODELES FEM) that provides a means of numerical analysis by a computer to various problems of applied mechanics, which was difficult with the node finite element method. It is another object of the present invention to provide a method capable of performing high-accuracy analysis in a short time even with a medium-scale computer or a personal computer, for example, without using a post-processor for large-scale analysis calculation.
[0022]
[Means for Solving the Problems]
In order to solve the above problem, a nodeless finite element method according to a first aspect of the present invention divides a load target area into small areas, calculates a displacement caused by a load for each of the divided small areas, and obtains the displacement by calculation. The displacement of each of the small regions is subjected to coordinate transformation based on the correlation with the internal boundary or the external boundary of each of the adjacent small regions, and the displacement of the entire target region is calculated based on the displacement of each of the coordinate transformed small regions. And
[0023]
Further, the nodeless finite element method of the second invention divides a load target area into small areas by overlapping structural grids having different sizes, and calculates a displacement caused by a load for each of the divided small areas. The coordinates of the displacement of each small region obtained by the calculation are converted by the correlation with the inner boundary or the outer boundary with the adjacent small regions, and the displacement of the entire target region is calculated based on the displacement of each coordinated small region. It is characterized in that it is calculated.
[0024]
Further, the nodeless finite element method of the third invention is a method of converting the displacement calculated by the nodeless finite element method of the first or second invention according to a predetermined displacement / stress conversion formula to each small area or object. It is characterized in that the stress of the entire region is obtained.
[0025]
Further, the nodeless finite element method of the fourth invention is an object obtained by converting the displacement calculated by the nodeless finite element method of the first or second invention by a predetermined displacement / boundary force conversion formula and dividing the object. It is characterized in that an internal boundary force and an external boundary force within a region are obtained.
[0026]
According to a fifth aspect of the present invention, in any one of the first to third nodeless finite element methods, the target region is a region defined on a two-dimensional or three-dimensional solid.
[0027]
BEST MODE FOR CARRYING OUT THE INVENTION
Prior to the description of each embodiment, first, a method of deriving each mathematical expression used for obtaining displacement and stress on a boundary side between elements without using a node will be described.
In order to obtain the displacement and stress on the boundary side between elements without using nodes, the nodeless finite element method introduces the concept of a nodeless finite element expressed by the following energy principle equation (1).
[0028]
(Equation 5)
Where u i : Displacement component, t i : Boundary force component, σ ij : Stress component,
Bar u i : Specified displacement bar t i : Boundary force acting as external force from surface load component or adjacent element
Bar p i : Body force component
It is.
[0029]
For more specificity, the formulation of the plane stress by the no-node finite element method using the triangular element of FIG. 2 as an example will be described.
FIG. 1 is a diagram showing an embodiment of a nodeless finite element (small area) in the nodeless finite element method. Reference numeral 10 denotes a triangle element as a nodeless finite element, and reference numerals 11 to 13 denote adjacent elements (nodeless finite element). Element). Reference numerals s1, s2, and s3 are sides of the no-node finite element 10, reference numerals s3, s4, and s5 are sides of the adjacent element 11, reference numerals s6, s7, and s1 are sides of the adjacent element 12, and reference numerals s8, s9, and s2 are Sides of the adjacent element 13 that are shared by the nodeless finite element 10 and the adjacent elements 11 to 13 are internal boundaries, and sides s4 to s9 are external boundaries. In this case, the displacement u of the point P inside the triangular element 10 in FIG. 1 is developed into a finite power series as in the following equation (Equation 6).
[0030]
(Equation 6)
Where a i , B j .. Are unknowns.
[0031]
Expression 6 is separated into coordinates and unknowns, and is represented by the following expression (Expression 7) in a matrix-vector format.
[0032]
(Equation 7)
[0033]
From Equation 7, the strain ε is given. Further, a strain-stress relational expression (σ = D (ε−ε 0 ) + Σ 0 D is an elastic matrix containing material constants), and stress σ is given. The stress σ at the point P of the triangular element 10 is given by the strain-stress relationship shown in the following equation 8, and the boundary force t at the boundaries s1 to s3 i Is given by the equation of Caushy (personal name) (formula 8 below).
[0034]
(Equation 8)
By substituting Equations 7, 8, and 9 into Equation 5, and performing a mathematical operation, the following Equation (Equation 10) can be obtained.
[0035]
(Equation 9)
[0036]
FIG. 2 is an example of a plane plate divided into four elements. In this example, the displacements u and v are fixed on the left side of the element 11, and the distributed load bar t in the x direction is distributed on the right side of the element 24. When the above equation 10 is used in this example, the following equation (equation 11) is obtained. Elements 21 to 24 are nodeless finite elements.
[0037]
(Equation 10)
Where M 11 , M 12 Is the matrix, a i Is an undetermined coefficient vector of the element 11.
Here, the following equation (Equation 12) is obtained by expressing Equation 11 in a matrix-vector representation.
[0038]
[Equation 11]
[0039]
Equation 12 is the final simultaneous linear equation to be calculated. Since the environmental conditions have already been introduced, a matrix operation is performed by a computer, and the coefficient value of the simultaneous linear equation ｛a 1 , A 2 , A 3 , A 4 By calculating｝, a distributed load bar t in the x direction distributed on the right side of the element 24 can be obtained.
[0040]
FIG. 3 is a diagram illustrating an example of an irregular mesh having a corner of another element on the side of the element.
The basic formula of the nodal finite element method according to the prior art is Formula 2 which is a balance equation at the node, and each line of the simultaneous linear equation (Formula 4) is a balance equation relating to the node. The irregular mesh as shown in FIG. 3 is desirable for mesh division performed in the simulation of plastic deformation, and is divided such that a corner of another element is arranged in the middle of the side of another element. However, in the conventional node finite element method, it is necessary to repeatedly calculate the equilibrium equation at the node and calculate the equilibrium equation at each node at the corner of the next element. Is required, and it takes a long calculation time, so that a large-scale post-processor must be used. In addition, the accuracy of the calculation result decreases each time the matrix operation is repeated, and the reliability of the final calculation result cannot be guaranteed.
[0041]
On the other hand, the no-node finite element method of the present invention divides the entire region into a finite number of elements, but does not need to establish the above-mentioned equilibrium equation of Equation 2 since no nodes are used (that is, in the example of FIG. 20). Node 320 has node vector f ) And the central common vector [P x , P y ] T However, in the examples of FIGS. 1 and 2, there is no need to establish an equilibrium equation because no nodes are used.) Therefore, even for an irregular mesh as shown in FIG. 3, it is not necessary to repeat the matrix operation many times unlike the nodal finite element method, so that a large-scale post-processor is not required, and the accuracy of the calculation result is not dependable. Quality can be guaranteed. In addition, an undetermined coefficient vector a i The calculation can be performed even if the numbers of are different. This means that the nodeless finite element method facilitates the configuration of a mesh-divided preprocessor and that the accuracy of the calculation result regarding stress concentration is high.
[0042]
Further, the nodal finite element method of the prior art is a discrete lumped system, but the no-node finite element method of the present invention discretizes the entire region into a finite number of elements. As shown in FIG. 4, the distribution of stress between elements by an analytical calculation based on a no-node finite element method using a second-order or higher displacement function is continuous.
[0043]
From the contrast between the simultaneous linear equation (Equation 4) representing the total nodal displacement vector in the nodal finite element method and the simultaneous linear equation (Equation 12) representing the stress distribution in the no-node finite element method, When creating a matrix, there is no superposition process. That is, in the nodeless finite element method, even if the order of the simultaneous linear equations of Expression 12 is rearranged, the order does not affect the calculation result.
[0044]
<Example>
1. Hardware configuration
FIG. 5 is a block diagram showing a configuration of an embodiment of a numerical analysis system capable of executing an analysis calculation by the nodeless finite element method. In this example, a personal computer (hereinafter, a personal computer) is used as the numerical analysis system.
In FIG. 5, the personal computer 5 includes a control unit 51, a pointing device 52 such as a mouse, an input device 53, a working memory 54, a storage memory 55, a display device 56, and a printer 57.
The control unit 51 has a configuration including a CPU (or MPU) and a ROM for storing programs, and a peripheral control circuit. The control unit 51 controls the operation of the entire apparatus, and performs the no-node analysis calculation stored in the storage memory 23. A program group (nodeless analysis means) controls execution of a nodeless analysis process based on the nodeless finite element method of the present invention.
[0045]
The pointing device 52 is a unit for designating a position when the area is designated by the person in charge of analysis, a drawing parameter, various operation instructions, and the like, and sends position information to the control unit 51 when an operation is performed. The input device 53 has a keyboard and a pointing device such as a mouse, and the user can input data and parameters.
[0046]
The work memory 54 is a temporary storage memory such as a RAM. When the personal computer is started, the control program stored in the program storage ROM is read and stored, and the solid-state mechanics analysis process is selected. The analysis program corresponding to each specific model is read out from the storage memory 55 as appropriate and stored.
[0047]
The storage memory 55 is made of a removable recording medium such as a magnetic disk or an FD, and is a program module group for nodeless finite element analysis created to realize the nodeless finite element analysis procedure of the present invention. In addition to analysis support programs, input data, parameters and analysis results (output data), accuracy evaluation data, target area setting programs including drawing programs, area division programs, matrix equation calculation programs by hollowing out operations, etc. Stores and records programs and support programs.
[0048]
2. model
Hereinafter, as an example of the numerical analysis based on the no-node finite element method using the personal computer 5, a stress concentration model which is difficult to calculate with the conventional finite element method (nodal finite element method) will be described. The applicable range of the element method is not limited to the stress concentration model. In other words, it is possible to obtain faster and more accurate calculation results for dynamic models that can be solved by the conventional finite element method, and also for applications in fields that were conventionally difficult such as stress concentration models. Fast and accurate results can be provided.
[0049]
A stress concentration model (considering how and where stress is concentrated and propagated) as an example of a rigid body displacement is, for example, a closed vessel having a pipe hole for cooling water recirculation such as a nuclear reactor, This is necessary when designing a structure (for example, an aircraft, a ship, a bridge, etc.) made of a partition wall and a screw-welded plate of a plant, and a tensile load model of a perforated plate as shown in FIG. A typical example.
[0050]
3. Procedure using no-node finite element method
FIG. 6 is a flowchart showing a model analysis procedure using the nodeless finite element method.
First, a material, a structure, and the like to be analyzed (hereinafter, referred to as a target region) are specified, and a part to be analyzed (hereinafter, referred to as a target region) and a unique region in the target region (for example, a flat plate). The perimeter of the perforated hole, the perimeter of the bent surface, the welding location and its surroundings) are specified, and the size, weight, etc., characteristics, the amount of load such as load, and the condition value such as the direction from the value (when changing as a parameter, , Including its pitch). The target area is usually two-dimensional (plane) or three-dimensional (three-dimensional), but is not limited to this (step S1).
[0051]
Next, the target region (including the unique region, if any) is modeled as shown in the example of FIG. In modeling, it is desirable to perform normalization with the entire model being 1, and furthermore, as shown in the example of FIG. 7B, the analysis target region (the example of FIG. It is computationally desirable to determine the portion indicated by reference numeral 71 '). Model creation and normalization can be easily performed using the personal computer 5 as described later (step S2).
[0052]
Next, the entire analysis target region of the model is divided into n × n (n is an integer of 2 or more) meshes (step S3), and the periphery of the singular region specified in step S1 is divided in step S3. The mesh is subdivided and divided into small meshes, each small mesh, that is, each element is numbered, and the length of each side (= the number of divisions) is associated. At this time, as shown in FIG. 8, the mesh may be subdivided into several stages so as to narrow down the unique region from the entire target region. Such mesh division can be performed by a simple operation using the personal computer 5 as described later (step S4).
[0053]
First, in the first stage (that is, step S3), for each small area (element) divided into meshes (square lattice in the case of two dimensions, cubic lattice in the case of three dimensions), condition values such as coordinate values and materials (Size, weight, characteristic value, load, direction of force, etc. (when the model is normalized, the condition values are also normalized)), and the above-described simultaneous linear equations of Expression 12 are calculated to obtain The displacement (undetermined coefficient vector) for each small area is calculated. That is, in the no-node finite element method, a result is calculated for each small area (step S5).
[0054]
Next, for each of the sub-regions subdivided for each stage among the subdivided meshes of the mesh obtained in the first stage, the above-described multiple simultaneous linear equation (square matrix equation) of Expression 12 is calculated to calculate the displacement for each region. (Undetermined coefficient vector) is calculated. In this case, the sum of the displacements given to each small area is the displacement of the small area (the small area before subdivision) calculated in step S5. If the small area is further subdivided, the same operation is repeated. That is, this is the same as overlapping meshes having different sizes. Since the no-node finite element method can treat each small area independently, the applicant of the present invention applied a matrix based on the invention filed on May 25, 1999 (the matrix compression processing (Japanese Patent Laid-Open No. 2000-339295)). Calculates the simultaneous linear equations while performing the compression and compression processing, and obtains an approximate solution of the undetermined coefficient vector (the displacement for each small area in the example of FIG. 7), and approximates the displacement of the analysis target area based on the approximate solution. A solution can be obtained (step S6).
[0055]
Next, the approximation accuracy is verified, and when the accuracy is within the allowable error, the analysis calculation is terminated. If the allowable error is not reached, the process returns to step S3 to further divide each element, or returns to step S1 to change the setting of the condition value. Note that the verification of the approximation accuracy in the present invention is performed by using “a means for converging an approximate solution in the no-node finite element method” filed on the same day by the present applicant (step S7).
[0056]
The stress can be obtained from the strain-stress relational expression (Equation 8) based on the displacement obtained in step S7. In addition, the boundary force can be obtained from the equation of Cauchy (person name).
[0057]
4. Stress concentration model
FIG. 7 is an explanatory diagram of a tensile load model of a perforated plate in which a uniform tensile load is applied to both ends of a square plate having a circular hole in the center. In FIG. Reference numeral 72 denotes a circular hole opened at the center of the plate 71, and reference numerals 73 and 73 'denote tensile loads S uniformly applied to both ends of the plate, and the tensile loads 73 and 73' are respectively applied in opposite directions. FIG. 7B is a diagram showing that the analysis target area is treated as a quarter of the plate 71 ′ because the plate 71 is biaxially symmetric. In this example, when the length of the plate is 2L, L = 50 mm (mm), the radius of the hole r = 10 mm, the tensile load S = 100 kgf / mm, E = 20,000 kgf / mm, and ν = 0.3. did. Reference numerals 74 to 75 denote displacements fixed at a boundary side (AB) with the adjacent region on the left side of the analysis target region 71 ′, and reference numerals 76 to 79 denote boundary sides with the lower adjacent region of the analysis target region 71 ′. (CD) represents a fixed displacement. Note that the displacements 75 and 76 are also displacements acting on the boundary side (arc BC (1/4 circumference)) with the circular hole 72.
[0058]
FIG. 8 shows an example in which the analysis target area (1/4 plate 71 'of FIG. 7B) is divided into areas as described in steps S3 and S4 of the flowchart of FIG. 6, and FIG. Example in which small areas (elements) are obtained by superimposing divisions based on a structural lattice (that is, 16 small areas 81 obtained by dividing a quarter plate 71 ′ by the largest structural lattice, The division of the small area close to the circle is repeated three more times at 1/2 scale, and the division of the small areas 82 and 83 close to the boundary with the circular hole (internal boundary area) 72 is sequentially advanced to a total of four times. This is an example in which 16 small areas 81, 20 small areas 82, 27 small areas 83, and 7 small areas 84, that is, a total of 73 small areas (elements) are obtained.
[0059]
FIG. 8B shows an example in which small regions 83 and 84 (including a boundary region with the circular hole) near the circular hole 72 are further divided into three to obtain 137 small regions (elements).
[0060]
In such a case, in the conventional finite element method (node finite element method), the division method is divided into 16 at first, and the 分割 division is repeated three times, resulting in 128 elements. Since the calculation must be performed on the elements, the operation by the large-scale computer must be repeated, and errors accumulate and the accuracy is reduced. However, in the dividing method of the present invention, only 73 elements are required as shown in FIG. In addition, since it is only necessary to analyze the elements in the subdivided range, the calculation is simplified, and the accuracy of the calculation result (approximate solution) is good. Further, in the conventional division method of the finite element method, first, 16 divisions are performed, and further two divisions are repeated three times, and then three divisions are performed once, resulting in 484 elements. However, the division method of the present invention requires only 137 elements as shown in FIG. 8B, and furthermore, only needs to perform calculations on the elements in the subdivided range.
[0061]
3. Calculation procedure by PC
FIG. 9 is a diagram showing an embodiment of a job selection menu displayed when the nodeless finite element analysis system is started, and FIG. 10 is a personal computer shown in FIG. 9 is a flowchart illustrating an example of a control operation performed by a control unit 51 when performing an analysis calculation of displacement and stress acting on a model. In the following description, the control unit 51 performs control and calculation according to a control program and an analysis program constituting the nodeless finite element analysis system.
[0062]
(Operation example)
When the personal computer 5 on which the no-node finite element analysis system including the control program, various analysis programs, screen display data, and the like is installed is started, the control unit 51 displays a job as shown in FIG. A selection menu is displayed (step T1). When the analysis person selects the two-dimensional no-node finite element analysis application from the displayed job selection menu (step T2), the control unit 51 reads out the target area setting program including the drawing program from the storage memory 55. Then, the program is stored in the application program area of the work memory 54, and the target area setting program is started, so that the person in charge of analysis is prompted to draw the target area (plane) (step T3).
[0063]
When the analyst draws a target area (effective plate in the example of FIG. 7) on the screen of the display device 56 with the mouse 52 and inputs a size, a unit, and the like by a key, the control unit 51 converts the drawn target area into a graphic. (See FIG. 7 (a)) and displayed, prompting the person in charge of analysis to specify the unique region, the portion where the load such as the load is applied, and the direction thereof (step T4).
[0064]
The singular area (circular hole 72 in the example of FIG. 7) on the target area, which is plotted by the analyst, and the portion or direction in which the load acts (tensile loads 73 and 73 ′ acting on the right and left sides in the example of FIG. 7). Is designated by the mouse 52 (step T5), the control unit 51 extracts an analysis target area (see FIG. 7B) based on the symmetry of the graphic by the analysis target area graphicization program, and Display on the screen. The display of the analysis target region is clearly indicated by surrounding the portion with a broken line or the like or displaying the target region in a color different from the display color of the target region. Alternatively, a window may be displayed. In addition, the analysis target area can be corrected or changed, or an alternative can be displayed according to an instruction of an analysis person. When there are a plurality of specific regions in the target region, the person in charge of analysis can specify a plurality of specific regions. When a plurality of unique regions are specified, the analysis target regions corresponding to the number of the unique regions are displayed so as to include the unique regions one by one, and the analysis person is urged to specify the order of the analysis calculation. At this time, the person in charge of analysis can fix a mouse cursor on a side, a point, or a surface of the figure, input a character or a symbol from a keyboard, and attach a name or a symbol to the side, point, or surface of the figure. (Step T6).
[0065]
Next, the control unit 51 reads out the area division program from the storage memory 55, overwrites and stores it in the application program area of the working memory 54, starts the area division program, sets k = 1, and The user is prompted to input the k-th area division number n of the designated analysis target area (step T7).
[0066]
When the person in charge of analysis inputs the k-th area division number n (5 in the example of FIG. 7), the control unit 51 uses the area division program to specify the analysis target area as n × n (see FIG. 7). Then, the mesh is divided into the mesh of 5), and a subdivision necessity selection button (or selection menu) or the like is displayed on the screen to urge the selection of the necessity of subdivision of the divided small area (element) (step T8).
[0067]
If the analysis technician selects “no re-division”, the process proceeds to step T12. If “re-division required” is selected, the control unit 51 further selects the analysis target area divided in k-th order in step T7. Since the user is prompted to specify a region range to be divided (step T9), the person in charge of analysis specifies the range including the singular region (a specific range of the mesh divided into primary regions) with the mouse, and furthermore, n k + 1 Is input, the (k + 1) th-order division is performed, and each mesh within the specified range is (n k + 1 ) × (n k + 1 ) Divide and further display a re-division necessity selection button (or selection menu) or the like on the screen to urge the user to select further necessity of re-division of the divided small area (element). That is, when k = 1, the 11th mesh (small area) 82-11 and the 23rd mesh 82-23 are specified in the example of FIG. 7A, and when the number of divisions 2 is input, the singular area ( The nine meshes 82 on the hole 72) side are further divided into 36 2 × 2 small regions 83, and a request is made to determine whether or not to re-divide the small region 83 on the hole 72 side (step). T10).
[0068]
Further, if the analysis staff selects “No re-division”, the process proceeds to step T12, and if “re-division is required”, 1 is added to k and the process returns to step T10 (step T11).
[0069]
When the area division by the no-node finite element method as shown in FIG. 8 is completed, the control unit 51 reads out the analysis support program, overwrites and stores the analysis support program in the application program area of the working memory 54, and starts the analysis support program. Each element number is assigned to each divided small area and displayed. Although the element number is usually automatically assigned, the person in charge of analysis can also assign a number by designating a small area (step T12).
[0070]
When the assignment of the element numbers is completed, the control unit 51 displays a guide message on the display device 56 to prompt the analyst to input data. Therefore, each small area (the size may be different) divided by the analyst. ) Is defined as a system of linear equations (see equation 12) using the undetermined coefficient vector as an undetermined coefficient vector (step T13), and initial values (in the examples of FIGS. 7 and 8, the size and the tensile load S acting on each element on the right side) Is introduced, the control unit 51 obtains the approximate value of the undetermined coefficient while performing the centering compression process disclosed in Japanese Patent Laid-Open No. 2000-339295, and calculates the displacement of the analysis target area from the obtained approximate value of the undetermined coefficient ( Step T14).
Next, it is determined whether or not the calculated approximate value is within a predetermined error range (allowable error). If the calculated approximate value is within the allowable error, the process proceeds to step T19, and if not, the process proceeds to T16 (step T15). .
[0071]
The control unit 51 displays the target element (small area) and its peripheral elements on the screen of the display device 24 of the personal computer 5, and further executes the same area division program used in step T8 from the storage memory 55. The data is read out, overwritten and stored in the application program area of the working memory 54, the area division program is started, and the person in charge of analysis prompts the input of the division number j of the target element (small area) (step T16).
[0072]
When the person in charge of the analysis inputs the number of divisions j, the control unit 51 divides the target element by a mesh using the area division program (step T17).
[0073]
When the division of the target element is completed, the control unit 51 reads the analysis support program, overwrites the analysis support program in the application program area of the working memory 54, starts the analysis support program, and stores each element in each divided small area. The number is reset, and the process returns to step T14. Thereby, a more detailed approximate solution of displacement can be obtained for the entire element including the subdivided element (step T18).
[0074]
When the calculation result (approximate solution) is within the allowable error, the result is displayed on the display device 53 as a solution. Further, it can be printed out by the printer 57 or written in the storage memory 55 as specified (step T19).
[0075]
Further, the stress can be obtained from the strain-stress relational expression (Expression 8) based on the approximate solution of the displacement obtained in Step T18. In addition, the boundary force can be obtained from the equation of Cauchy (person name).
[0076]
FIG. 11 is a diagram showing the distribution of the displacement between YY ′ of the tensile load model of the perforated plate by the no-node finite element method, and illustrates the result obtained by the personal computer 5 based on the flowchart of FIG. 7A, that is, between 10 mm (upper part of the perforated portion) and 50 mm (end) in the Y-axis direction in the 1/4 plate 71 ′ of FIG. 7B. 3 shows the distribution of the stresses generated in FIG.
[0077]
As is clear from FIG. 11, the calculation result by Treftz (person's name) matches the calculation result by the no-node finite element method. In other words, the perforated plate model of FIG. 7 is an irregular mesh obtained by overlapping and dividing regions as shown in FIG. 8, and the irregular mesh which was impossible in the conventional nodal finite element method in terms of calculation time and accuracy. This indicates that the calculation can be easily performed with high accuracy. In FIG. 11, it is clear that the displacement near the hole (at a position of 10 mm) is the largest, and the displacement decreases as the distance increases.
[0078]
FIG. 12 is a view showing a state of deformation of a plate 71 'when a tensile load S is applied to a quarter plate 71' divided into regions as shown in FIG. 7B, that is, a perforated hole by the no-node finite element method. FIG. 13 is a diagram illustrating an example of a simulation of the distribution of displacement between YY ′ of the tensile load model of the plate, and FIG. 12 shows that deformation near the circular hole 72 is remarkable. In the simulation, a two-dimensional analysis simulation or a three-dimensional analysis simulation is selected from the job selection menu in step T2 of the flowchart of FIG. 10, and an initial value and a pitch of a variable (parameter) are designated, and the process proceeds to steps T14 to T18 in FIG. It can be performed by repeating the corresponding operation.
[0079]
Although the drawing program is used in step T3, it may be a CAD program. Further, in FIG. 8, the shape of the mesh when the area is divided into meshes in steps T8 and T10 is shown as an example of a square lattice, but is not limited to the square lattice.
[0080]
Further, the case where the two-dimensional no-node finite element analysis application is selected in step T2 has been described. However, in the case of three-dimensional no-node finite element analysis, an approximate solution can be obtained by the same operation as the operation after step T3. . In this case, the analysis target area is divided by an n × n × n solid grid (not limited to a square solid grid) in step T8. Further, if necessary, it can be re-divided in step T10.
[0081]
In the above example, an example using a personal computer has been described, but the analytical calculation device according to the present invention is not limited to a personal computer. That is, it goes without saying that the present invention can be applied to a medium-sized, large-sized, or large-scale computer. That is, a solution model may be created according to the scale of the solid mechanics problem that needs to be solved, and the scale of the computer may be appropriately selected. Hereinafter, another embodiment of the analyzer will be described.
[0082]
FIG. 13 is a block diagram showing the configuration of an embodiment of an analysis calculation apparatus capable of executing calculations by the no-node finite element method according to the present invention, in which parallel processing is performed by a multiprocessor. In FIG. 13, a multiprocessor 130 includes a plurality of processing devices (processors) 131-1, 131-2,..., 131-n (n = 16 in the embodiment) sharing a storage device 132 and an input / output device 133. ing. In addition, an external device (in this example, an electronic computer such as a personal computer (hereinafter, a personal computer)) 135 is connected to the input / output device 133. The input / output device 133 includes an input device such as a keyboard, an output device such as a printer and a monitor device, and a magnetic disk and an optical disk device.
[0083]
With the above hardware configuration, the external device 135 executes steps T1 to T13 and T15 to T19 in the flowchart of FIG. 10, and the multiprocessor 130 executes step T14. In the example shown in FIG. 13, the multiprocessor 130 and the external device 135 are connected by a cable or the like. However, the multiprocessor 135 and the external device 135 may be separated from each other. It may be performed.
[0084]
FIG. 14 is a block diagram showing the configuration of another embodiment of the analysis calculation apparatus capable of executing the calculation by the nodeless finite element method according to the present invention. The configuration in the case of performing the analysis calculation by the parallel processing using the communication network. Here is an example.
[0085]
In FIG. 14, a network 140 includes computers 142-1, 142-2,..., 142-n (n = 6 in the embodiment) connected to the communication line 145, and computer devices 142-1 and 142-2, ,..., 142-n via the communication line 145, and a server 141 for controlling data exchange between the computers 142-1, 142-2,. have. Note that the network 145 may be a communication network such as the Internet, or may be a network including a dedicated line connecting the research laboratories, etc., or may be a mixed type of these.
[0086]
The server 141 includes a computer having a processing device, a communication control unit, a memory, and an input / output device, and can connect to a database.
[0087]
The computers 142-1, 142-2,..., 142-n are not limited to large computers, and have a processing speed and a memory capacity sufficient to perform a square matrix operation of a predetermined order, and have a communication control function. For example, a personal computer or a workstation may be used.
[0088]
In the case of performing a numerical analysis on the request of the square matrix distribution, if the processing speed and the memory capacity of each computer are far apart, the efficiency of the computer is limited by the computer having the minimum capacity (for example, when the computer is connected to the network 145). .., 142-n are personal computers, the processing efficiency of the entire system is limited by the processing efficiency of the personal computers. Therefore, it is desirable that the difference between the abilities is within a certain range.
[0089]
When a communication network is used, the analysis calculation can be performed in a procedure according to the flowchart of FIG.
In this case, steps T1 to T12 are executed on the server 141 side, and the initial value and matrix of step T13 are distributed from the server to the computers 142-1, 142-2,... , 142-n are executed by each of the computers 142-1, 142-2,..., 142-n. After the execution is completed, the approximate solution is transmitted to the server 141, and the server 141 executes steps T15 to T19. Can be executed.
[0090]
In the example of FIG. 14, the server 141 and the terminal are different devices, but one of the terminals (computers 142-1 to 142-n) (for example, the computer 142-1) is also used as the server. You may do so.
[0091]
(Application of nodeless finite element method to rigid body displacement model)
The nodeless finite element method can handle rigid body displacement clearly as described below.
A plane stress model will be described as an example. First, a displacement function of the following Expression 13 is assumed.
[0092]
(Equation 12)
[0093]
Where u: displacement in x direction ν: displacement in y direction
u 0 : Rigid displacement in the x direction ν 0 : Rigid body displacement in y direction
χ 0 ： Rigid body rotation in plane ε x0 , Ε y0 , Γ xy0 : Strain component
Represents
[0094]
Equation 13 expresses the displacement of an arbitrary point in the element (small area) separately into terms caused by strain and rigid displacement. Where u 0 , Ν 0 , Ε x0 , Ε y0 , Γ xy0 Is the unknown for each element, and the total unknown is 6 × the number of elements. When the model of a cantilever flat plate on which a shear load is applied is analyzed by combining Expression 13 with the following Expression 14, which represents a new energy principle, FIGS. 15 and 16 are obtained.
[0095]
FIG. 15 is a diagram showing the deformation (displacement) of the cantilever flat plate to which a shear load acts, and FIG. 16 is a diagram showing the stress distribution of the cantilever flat plate to which a shear load acts. FIGS. 15 and 16 show that in the case of 40 × 20 element division, both displacement and stress have sufficient accuracy.
[0096]
(Equation 13)
[0097]
Next, from Equation 13, the rigid body rotation χ 0 When the same model is analyzed by using the following displacement function (Equation 15) in which the term is intentionally dropped, it can be seen that both the displacement and the stress become extremely small. This numerical experiment shows that all rigid body displacements should be handled explicitly.
[0098]
[Equation 14]
[0099]
Next, the following equation 16 shows the bending function of the plate bending moment including the rigid body displacement explicitly.
[0100]
[Equation 15]
[0101]
FIGS. 17 and 18 show the results of the analysis calculation using Expression 16. FIG. 17 and FIG. 18 show the calculation results of the deflection of the simple support plate that is affected by the uniformly distributed load, and the number of mesh divisions is 25 × 25 elements. FIG. 18 is a diagram obtained by using a technique for loosening the coupling force between elements when performing numerical calculation.
[0102]
In the nodeless finite element method of the present invention, it can be seen that, after deformation, individual elements (small areas) are isolated. That is, the calculation can be performed even if the elements are not attached to each other. Also, the partial diagonal matrix M of the no-node finite element method ij (See Equation 12) is symmetric and has an inverse matrix. Therefore, a vector elimination method or a vector repetition method of simultaneous linear equations is established. This indicates that the data can be deleted for each element, and iterative calculation can be performed for each element. Such a property of high computational efficiency exerts its power in the analytical calculation of a nonlinear model.
[0103]
As mentioned above, although one Example of this invention was described, this invention is not limited to the said Example, It cannot be overemphasized that various deformation | transformation implementation is possible.
[0104]
【The invention's effect】
As described above, according to the no-node finite element method of the present invention, the target area is discretized into finite arc elements. The stress obtained by the node finite element analysis has a continuous distribution between elements (small areas). In addition, the no-node finite element method does not have a superposition process like the nodal finite element method when creating the whole matrix, so an irregular mesh including small areas divided by overlapping structural grids of different sizes Numerical analysis can also be easily performed for.
[0105]
In the nodeless finite element method of the present invention, since individual elements (small areas) are isolated after deformation, that is, calculation can be performed even if elements do not adhere to each other. Also, the partial diagonal matrix M of the no-node finite element method ij (See Equation 12) is symmetric and has an inverse matrix. Therefore, a vector elimination method or a vector repetition method of simultaneous linear equations is established. This indicates that the data can be deleted for each element, and iterative calculation can be performed for each element. Such a property of high computational efficiency exerts its power in the analytical calculation of a nonlinear model.
[0106]
In addition, since the calculation can be performed even when the small areas are not adjacent to each other, the order of the simultaneous linear equations (vector equations) is not changed even if the order is rearranged. Considering this property and the fact that superposition is not necessary, if a program is created, it is possible to calculate even if the small area (element) is destroyed in a dynamic sense and the total number of elements increases dynamically. Become.
[Brief description of the drawings]
FIG. 1 is a diagram showing one embodiment of a nodeless element (small area) in a nodeless finite element method.
FIG. 2 is an example of a plane plate divided into four elements.
FIG. 3 is a diagram illustrating an example of an irregular mesh having a corner of another element on a side of the element;
FIG. 4 is a diagram showing an example of a stress distribution between elements by an analytical calculation based on a nodeless finite element method.
FIG. 5 is a block diagram showing a configuration of an embodiment of a numerical analysis system capable of executing an analysis calculation by a nodeless finite element method.
FIG. 6 is a flowchart illustrating a model analysis procedure using a nodeless finite element method.
FIG. 7 is an explanatory diagram of a tensile load model of a perforated plate.
FIG. 8 is a diagram showing an example of area division.
FIG. 9 is a diagram illustrating an example of a job selection menu displayed when the nodeless finite element analysis system is activated.
FIG. 10 is a flowchart illustrating an example of a control operation performed by a control unit when an analysis calculation of displacement and stress acting on a tensile load model of a perforated plate is performed by a personal computer.
FIG. 11 is a diagram showing a distribution of displacement between YY ′ of a tensile load model of a perforated plate by the no-node finite element method.
FIG. 12 is a diagram showing an example of a simulation of a displacement distribution between YY ′ of a tensile load model of a perforated plate by a no-node finite element method.
FIG. 13 is a block diagram showing a configuration of an embodiment of an analysis calculation apparatus capable of executing calculation by the nodeless finite element method of the present invention.
FIG. 14 is a block diagram showing a configuration of another embodiment of the analysis calculation device capable of executing calculation by the nodeless finite element method of the present invention.
FIG. 15 is a diagram showing deformation (displacement) of a cantilever flat plate to which a shear load is applied.
FIG. 16 is a diagram showing a stress distribution of a cantilever flat plate to which a shear load is applied.
FIG. 17 is a diagram illustrating an example of a calculation result of a deflection of a simple support plate that is affected by an evenly distributed load.
FIG. 18 is a diagram illustrating an example of a calculation result of a deflection of a simple support plate that is affected by an evenly distributed load.
FIG. 19 shows an image of a planar area problem such as a film.
FIG. 20 is an explanatory diagram of displacement and concentration force acting on a node.
FIG. 21 is a diagram showing one embodiment of a series spring problem.
FIG. 22 is an explanatory diagram of a mechanical boundary condition and a geometric boundary.
[Explanation of symbols]
10, 11, 12, 13 elements (small area)
71 Perforated plate (target area)
71 'Analysis target area
81, 82, 83, 84 mesh (small area, element)
s1, s2, s3 inner boundary
s4-s9 outer boundary
