(57)【要約】【課題】対象物の形状を近似的に表現するとともにその対象物の変形を有限要素法によって解析するために用いられる要素の幾何学的なゆがみに起因して有限要素法による解析結果に生じる解析誤差を定量的に推定する。【解決手段】対象物の変形時にその要素に発生する変形モード（Ｓ２０１における変形モード値ｊ）を引張りと曲げとせん断とねじりとのうちの少なくとも一つとして想定するとともに、その想定された各変形モードにつき、要素のゆがみであるアスペクト比Ａ、反り度Ｗ、斜傾度Ｓおよび台形度Ｔに関連付けて、そのゆがみに起因する解析誤差をひずみエネルギ誤差Ｅｅとして推定する（Ｓ２０２ないしＳ２０８）。【選択図】  図４３
2004303227.tif 000002
【特許請求の範囲】
【請求項１】
対象物の形状を近似的に表現するとともにその対象物の変形を有限要素法によって解析するために用いられる要素の幾何学的なゆがみとそれに起因して有限要素法による解析結果に生じる解析誤差との関係を要素ゆがみ−解析誤差間関係としてコンピュータを用いて推定する要素ゆがみ−解析誤差間関係推定方法であって、
前記対象物の変形時に前記要素に発生する変形モードを引張りと曲げとせん断とねじりとのうちの少なくとも一つとして想定するとともに、その想定された各変形モードにつき、前記ゆがみに関連付けて前記解析誤差を推定する推定工程を含む要素ゆがみ−解析誤差間関係推定方法。
【請求項２】
前記要素が、面を有する面状要素であり、
前記推定工程が、
（ａ）前記面状要素をアイソパラメトリック要素に座標変換するとともに、前記面状要素の幾何学的なゆがみである斜傾性および台形性を互いに独立とするような複数の基本ゆがみ変数によって表現し、
（ｂ）そのアイソパラメトリック要素として前記面状要素を扱うことにより、その面状要素の変形モードを引張り、せん断および曲げを含む複数の基本変形モードに分離し、
（ｃ）前記アイソパラメトリック要素として前記面状要素を扱うことにより、その面状要素の垂直ひずみおよびせん断ひずみの各成分を前記複数の基本ゆがみ変数と前記複数の基本変形モードとを用いて定式化し、それにより、前記面状要素のひずみエネルギを計算し、
（ｄ）その計算されたひずみエネルギを用いることにより、各基本変形モードごとに、各基本ゆがみ変数の値に関連付けて前記解析誤差を計算する第１の計算工程を含む請求項１に記載の要素ゆがみ−解析誤差間関係推定方法。
【請求項３】
前記推定工程が、前記想定された各変形モードにつき、前記要素のひずみエネルギの解析解であって有限要素解析によることなく前記要素の領域内の応力およびひずみの厳密解を用いて計算されたものと、前記要素のひずみエネルギの近似解であって有限要素近似により計算されたものとの差を用いて前記解析誤差を計算する第２の計算工程を含む請求項１に記載の要素ゆがみ−解析誤差間関係推定方法。
【請求項４】
前記第２の計算工程が、前記対象物の形状を互いに共同して近似的に表現する複数の要素を１つの要素によって代表させるとともに、その代表要素について前記ひずみエネルギの解析解と近似解とを、前記代表要素は考慮するが他の要素は考慮しないで計算する解計算工程を含む請求項３に記載の要素ゆがみ−解析誤差間関係推定方法。
【請求項５】
前記ゆがみが、前記要素の幾何学的特徴を互いに独立して表現する複数の幾何学変数によって表現される請求項１ないし４のいずれかに記載の要素ゆがみ−解析誤差間関係推定方法。
【請求項６】
前記複数の幾何学変数が、少なくとも、前記要素の実形状のアスペクト比と、前記要素の実形状の、それと同じアスペクト比を有する基本長方形に対する斜傾度と、前記要素の実形状の、前記基本長方形に対する台形度とを含む請求項５に記載の要素ゆがみ−解析誤差間関係推定方法。
【請求項７】
前記斜傾度が、前記要素の実形状の、前記基本長方形に対するスキュー量を用いて定義され、
前記台形度が、前記要素の実形状の、前記基本長方形に対するテーパ量を用いて定義される請求項６に記載の要素ゆがみ−解析誤差間関係推定方法。
【請求項８】
前記複数の幾何学変数が、さらに、前記要素の実形状の、前記基本長方形に対する反り度を含む請求項６または７に記載の要素ゆがみ−解析誤差間関係推定方法。
【請求項９】
前記推定工程が、前記複数の幾何学変数の各々に関連付けて前記解析誤差を個別解析誤差として推定し、それら複数の幾何学変数に関連付けてそれぞれ推定された複数の個別解析誤差に基づき、前記ゆがみの全体に起因する前記解析誤差を総合解析誤差として推定するものである請求項５ないし８のいずれかに記載の要素ゆがみ−解析誤差間関係推定方法。
【請求項１０】
対象物の形状を近似的に表現する複数の要素の組合体を用いて前記対象物の変形を有限要素法によって解析した場合にその解析結果に生じる解析誤差をコンピュータを用いて推定する解析誤差推定方法であって、
前記コンピュータのユーザの変形モード特定指令に従い、各要素ごとに、前記対象物の変形時に各要素に発生することが予想される変形モードを引張りと曲げとせん断とねじりとのうちの少なくとも一つとして特定する変形モード特定工程と、
前記変形モードの特定後、前記複数の要素のうちの少なくとも一つにつき、各要素の幾何学的なゆがみに起因した前記解析誤差を定量的に表示する表示工程と
を含む解析誤差推定方法。
【請求項１１】
前記組合体に対して有限要素法による数値解析が実施されるのに先立って実行される請求項１０に記載の解析誤差推定方法。
【請求項１２】
さらに、各要素ごとに、前記特定された変形モードに基づき、かつ、要素の幾何学的なゆがみとそれに起因して有限要素法による解析結果に生じる解析誤差との間に予め推定された要素ゆがみ−解析誤差間関係に従って前記解析誤差を前記コンピュータを用いて推定する推定工程を含み、
前記表示工程が、その推定された解析誤差を表示するものである請求項１０または１１に記載の解析誤差推定方法。
【請求項１３】
前記要素ゆがみ−解析誤差間関係が、請求項１ないし９のいずれかに記載の要素ゆがみ−解析誤差間関係推定方法により推定されたものである請求項１２に記載の解析誤差推定方法。
【請求項１４】
前記ゆがみが、前記要素の幾何学的特徴を互いに独立して表現する複数の幾何学変数によって表現され、
前記推定工程が、それら複数の幾何学変数の各々に関連付けて前記解析誤差を個別解析誤差として推定する個別解析誤差推定工程を含み、
前記表示工程が、前記少なくとも一つの要素の各々に関連付けて、前記複数の幾何学変数の各々について、前記推定された個別解析誤差を表示するものである請求項１２または１３に記載の解析誤差推定方法。
【請求項１５】
前記推定工程が、前記少なくとも一つの要素の各々に関連付けて、前記ゆがみの全体に起因する前記解析誤差を総合解析誤差として推定する総合解析誤差推定工程を含み、
前記表示工程が、前記少なくとも一つの要素の各々に関連付けて、前記推定された総合解析誤差を表示するものである請求項１２ないし１４のいずれかに記載の解析誤差推定方法。
【請求項１６】
前記表示工程が、前記複数の要素を前記推定された総合解析誤差が許容値以上であるものとそうではないものとに区別して前記コンピュータのユーザに認識させるための表示を行う区別表示工程を含む請求項１５に記載の解析誤差推定方法。
【請求項１７】
対象物の形状を近似的に表現する複数の要素の組合体を用いて前記対象物の変形を有限要素法によってコンピュータを用いて解析する数値解析方法であって、
前記対象物の形状を前記複数の要素に前記コンピュータを用いて分割する要素分割工程と、
前記複数の要素への分割後、前記コンピュータのユーザの変形モード特定指令に従い、各要素ごとに、前記対象物の変形時に各要素に発生することが予想される変形モードを引張りと曲げとせん断とねじりとのうちの少なくとも一つとして特定する変形モード特定工程と、
前記変形モードの特定に応答し、前記複数の要素のうちの少なくとも一つにつき、各要素の幾何学的なゆがみに起因した解析誤差を定量的に表示する表示工程と、
前記組合体を用いて前記対象物の変形を有限要素法によって前記コンピュータによって解析する解析工程と
を含む数値解析方法。
【請求項１８】
請求項１ないし１７のいずれかに記載の方法を実施するためにコンピュータにより実行されるプログラム。
【請求項１９】
請求項１８に記載のプログラムをコンピュータ読み取り可能に記録した記録媒体。
【発明の詳細な説明】
【技術分野】
【０００１】
本発明は、対象物の形状を近似的に表現するとともにその対象物の変形を有限要素法によって解析するために用いられる要素の幾何学的なゆがみに起因して有限要素法による解析結果に生じる解析誤差を推定する技術に関するものである。
【背景技術】
【０００２】
対象物の変形を解析する手法として有限要素法が存在する。この有限要素法においては、対象物が複数の要素の組合体として近似的に表現され、それら複数の要素の組合体を用いて対象物の変形が解析される。そして、この有限要素法の一例として、対象物の形状を近似的に表現する要素を用いてその対象物の変形を解析する手法が存在する。要素は、例えば、膜要素として定義したり、薄肉を有するシェル要素として定義したり、四辺形を成す面を有する要素として定義したり、三角形を成す面を有する要素として定義することが可能である。
【０００３】
有限要素法による解析精度は対象物の要素分割に依存する。そのため、有限要素法の解析誤差を縮減するため、例えば、対象物を複数の要素に、対象物上の要素位置に応じて要素の大きさが適正となるように分割したり、要素の形状が適正となるように分割することが行われる。（例えば、特許文献１および２参照。）
【０００４】
特許文献１には、解析対象となる構造物を有限個のメッシュ（要素）に分割した後にメッシュ分割を最適化するために、応力、ひずみエネルギ等の解析誤差を考慮して、解析対象に設定された節点を移動する技術が記載されている。
【０００５】
特許文献２には、要求される解析精度を満たすサイズを有するメッシュを生成するとともに、解析対象上のメッシュ位置に応じてメッシュ分割の種類を、均等メッシュ分割、不均等メッシュ分割およびジャンクションメッシュ分割のうちのいずれかに設定する技術が記載されている。
【０００６】
要素形状が適正となるように要素分割を行う技術の一従来例によれば、同じ対象物の表面が分割されたすべての要素（四辺形を成す要素）ができる限り正方形すなわちゆがみのない四辺形を成すように要素分割が行われる。
【０００７】
さらに、その要素分割によって生成された各要素の形状が適正であるか否かを評価するために、各要素の形状の幾何学的なゆがみであって上記要素分割に起因して発生するものが、各要素のアスペクト比と反り度と斜傾度と台形度という４つの幾何学的特徴量によって表現される。
【０００８】
アスペクト比は、例えば、ゆがんだ要素の基準点（ゆがんだ要素の概して中央に位置する点）とそのゆがんだ要素の４つの要素辺のうちの第１要素辺の中点との第１距離と、基準点と、それら４つの要素辺のうち第１要素辺に隣接する２つの要素辺のいずれかである第２要素辺の中点との第２距離との比として定義される。
【０００９】
反り度（ワープ）は、例えば、ゆがんだ要素に対して投影平面を、そのゆがんだ要素の４つの節点から等距離にあるように想定し、その距離（投影高さ）の、ゆがんだ要素の４つの要素辺のうちのいずれかの長さの半分に対する比率の逆正弦として定義される。
【００１０】
斜傾度は、例えば、９０度から、前記基準点と前記第１要素辺の中点とを通過する第１直線と、基準点と前記第２要素辺の中点とを通過する第２直線との成す角度（小さい方）を引き算した値として定義される。
【００１１】
台形度は、例えば、ゆがんだ要素の各要素辺ごとに、前記基準点と各要素辺の２つの節点とを結ぶ２つの直線辺と対応する要素辺とによって三角形を作成し、最終的に作成された４つの三角形のうち面積が最小のものの面積を４倍したものをそれら４つの三角形の合計面積で割り算した値として定義される。
【００１２】
そして、前記従来例によれば、各要素ごと、各幾何学的特徴量について実際値と許容値とが互いに比較され、ある要素について少なくとも一つの幾何学的特徴量の実際値が許容値を超えると、その要素の形状が適正ではないと判定される。
【特許文献１】特開平４−２６８６７４号公報
【特許文献２】特開平１１−２５２９２号公報
【発明の開示】
【発明が解決しようとする課題】
【００１３】
しかし、この従来例は、幾何学的特徴量とそれに起因して解析結果に生ずる解析誤差との関係を明確にすることなく、経験則により、幾何学的特徴量の実際値から各要素の形状の適否を判定していた。
【００１４】
そのため、この従来例を採用する場合には、解析誤差が実際には小さい形状を有する要素についてその形状が適正ではないとの過った判定がなされ、その結果、不要な要素修正に時間が消費されてしまう可能性がある。
【００１５】
さらに、この従来例を採用する場合には、解析誤差が実際には大きい形状を有する要素についてその形状が適正であるとの誤った判定がなされ、その結果、必要な要素修正を経ないで得られた解析結果を解析者が信頼してしまう可能性もある。
【００１６】
本発明者らは、幾何学的特徴量と解析誤差との関係について研究を行い、その結果、解析対象の変形時に各要素に発生する変形モードが特定されれば、幾何学的特徴量と解析誤差との関係が数学的に誘導できることに気が付いた。
【００１７】
このような知見を背景にして、本発明は、各要素の幾何学的なゆがみに起因した解析誤差を定量的に推定することを課題としてなされたものである。
【課題を解決するための手段】
【００１８】
本発明によって下記の各態様が得られる。各態様は、項に区分し、各項には番号を付し、必要に応じて他の項の番号を引用する形式で記載する。これは、本発明が採用し得る技術的特徴の一部およびそれの組合せの理解を容易にするためであり、本発明が採用し得る技術的特徴およびそれの組合せが以下の態様に限定されると解釈すべきではない。すなわち、下記の態様には記載されていないが本明細書には記載されている技術的特徴を本発明の技術的特徴として適宜抽出して採用することは妨げられないと解釈すべきなのである。
【００１９】
さらに、各項を他の項の番号を引用する形式で記載することが必ずしも、各項に記載の技術的特徴を他の項に記載の技術的特徴から分離させて独立させることを妨げることを意味するわけではなく、各項に記載の技術的特徴をその性質に応じて適宜独立させることが可能であると解釈すべきである。
【００２０】
（１）  対象物の形状を近似的に表現するとともにその対象物の変形を有限要素法によって解析するために用いられる要素の幾何学的なゆがみとそれに起因して有限要素法による解析結果に生じる解析誤差との関係を要素ゆがみ−解析誤差間関係としてコンピュータを用いて推定する要素ゆがみ−解析誤差間関係推定方法であって、
前記対象物の変形時に前記要素に発生する変形モードを引張りと曲げとせん断とねじりとのうちの少なくとも一つとして想定するとともに、その想定された各変形モードにつき、前記ゆがみに関連付けて前記解析誤差を推定する推定工程を含む要素ゆがみ−解析誤差間関係推定方法。
【００２１】
後に理論的根拠を詳述するが、本発明者らは、任意の要素の幾何学的なゆがみは、複数の幾何学的特徴に着目することにより、完全に定義できるとともに、それら複数の幾何学的特徴は互いに独立しているという知見を得た。この知見によれば、要素のゆがみが正確に表現できることとなる。
【００２２】
力学的なつり合い条件から任意の仮想変位を考えると、外力のした仕事すなわち仮想変位と外力との積はひずみエネルギの増分に等しい。これは、仮想変位の原理と称される。線形弾性については、この仮想変位の原理は最小ポテンシャルエネルギの定理と等価なものとなる。この仮想変位の原理をエネルギによって表すと、正しい変位の解はポテンシャルエネルギと称される変位の関数を最小にするという最小ポテンシャルエネルギの定理が得られる。この定理によれば、真の変位ｕに対してポテンシャルエネルギＩが最小となる。ポテンシャルエネルギＩは、要素内に蓄えられるひずみエネルギＵと外力のポテンシャルエネルギとを含んでいる。この最小ポテンシャルエネルギの定理は有限要素法の理論的基礎である。
【００２３】
このように、ひずみエネルギと変位との間に依存関係が存在し、よって、ひずみエネルギの計算値に誤差が発生すれば、変位の計算値すなわち変形の解析結果にも誤差が発生することになる。
【００２４】
後に理論的根拠を詳述するが、要素のゆがみ程度が増加するにつれて、その要素に蓄えられるひずみエネルギも増加する。また、要素の変形モードが与えられれば、その要素のひずみエネルギの解析解を、有限要素解析によることなく、要素の応力およびひずみの理論解（厳密解）を用いて計算することができる。さらに、要素の変形モードが与えられれば、その要素のひずみエネルギの近似解を、有限要素近似により、その要素の節点変位の理論解（厳密解）を用いて計算することができる。それら解析解と近似解との差はひずみエネルギ誤差を意味する。このひずみエネルギ誤差は、上述の説明から明らかなように、有限要素法による解析結果の誤差につながる。
【００２５】
以上説明した知見に基づき、本項に係る方法においては、対象物の変形時に要素に発生する変形モードが引張りと曲げとせん断とねじりとのうちの少なくとも一つとして想定されるとともに、その想定された各変形モードにつき、ゆがみに関連付けて解析誤差が推定される。
【００２６】
ここに、要素の「ゆがみ」は、その要素の実形状の、基本形状（例えば、基本正方形）からのずれを意味する。
【００２７】
（２）  前記要素が、面を有する面状要素であり、
前記推定工程が、
（ａ）前記面状要素をアイソパラメトリック要素に座標変換するとともに、前記面状要素の幾何学的なゆがみである斜傾性および台形性を互いに独立とするような複数の基本ゆがみ変数によって表現し、
（ｂ）そのアイソパラメトリック要素として前記面状要素を扱うことにより、その面状要素の変形モードを引張り、せん断および曲げを含む複数の基本変形モードに分離し、
（ｃ）前記アイソパラメトリック要素として前記面状要素を扱うことにより、その面状要素の垂直ひずみおよびせん断ひずみの各成分を前記複数の基本ゆがみ変数と前記複数の基本変形モードとを用いて定式化し、それにより、前記面状要素のひずみエネルギを計算し、
（ｄ）その計算されたひずみエネルギを用いることにより、各基本変形モードごとに、各基本ゆがみ変数の値に関連付けて前記解析誤差を計算する第１の計算工程を含む（１）項に記載の要素ゆがみ−解析誤差間関係推定方法。
【００２８】
本項における「面状要素」は、狭義の面要素、板要素、膜要素、シェル要素等を意味するように解釈することが可能である。一般に、狭義の面要素は通常、面内変形を伴う二次元要素と称され、板要素は、横変形（面外変形）(transverse deformation)を受け、そして、シェル要素は、面内変形と横変形（面外変形）(transverse deformation)との双方を含む。
【００２９】
この定義に従い、本項における「面状要素」は、例えば、平面に投影されるシェル要素に分類される。
【００３０】
本項において「面状要素の垂直ひずみおよびせん断ひずみの各成分を定式化する」とは、例えば、長方形要素の幾何学変数を用いて面状要素の垂直ひずみとせん断ひずみとを展開することを意味するように解釈することが可能である。ひずみの成分には、例えば、垂直ひずみ、面内せん断ひずみ、面外ひずみ(transverse strain)等がある。
【００３１】
（３）  前記推定工程が、前記想定された各変形モードにつき、前記要素のひずみエネルギの解析解であって有限要素解析によることなく前記要素の領域内の応力およびひずみの厳密解を用いて計算されたものと、前記要素のひずみエネルギの近似解であって有限要素近似により計算されたものとの差を用いて前記解析誤差を計算する第２の計算工程を含む（１）項に記載の要素ゆがみ−解析誤差間関係推定方法。
【００３２】
この方法によれば、想定された各変形モードにつき、要素のひずみエネルギの解析解と近似解との差を用いて解析誤差が計算される。解析解は、有限要素解析によることなく、要素の応力状態の厳密解を用いて計算されたものである。これに対し、近似解は、有限要素近似により計算されたものである。
【００３３】
（４）  前記第２の計算工程が、前記対象物の形状を互いに共同して近似的に表現する複数の要素を１つの要素によって代表させるとともに、その代表要素について前記ひずみエネルギの解析解と近似解とを、前記代表要素は考慮するが他の要素は考慮しないで計算する解計算工程を含む（３）項に記載の要素ゆがみ−解析誤差間関係推定方法。
【００３４】
対象物の形状を互いに共同して近似的に表現する複数の要素の各々につき、その要素のみならず、有限要素近似についての別のスキームに基づく他の要素をも考慮してひずみエネルギの解析解と近似解とを計算することは可能である。しかし、それら複数の要素を代表する１つの要素は考慮するが他の要素は考慮しないでひずみエネルギの解析解と近似解とを計算することは可能であり、それらの計算値を用いて前記解析誤差を許容精度範囲内で推定することも可能である。
【００３５】
このような知見に基づき、本項に係る方法が提案されたのである。
【００３６】
（５）  前記ゆがみが、前記要素の幾何学的特徴を互いに独立して表現する複数の幾何学変数によって表現される（１）ないし（４）項のいずれかに記載の要素ゆがみ−解析誤差間関係推定方法。
【００３７】
（６）  前記複数の幾何学変数が、少なくとも、前記要素の実形状のアスペクト比と、前記要素の実形状の、それと同じアスペクト比を有する基本長方形に対する斜傾度と、前記要素の実形状の、前記基本長方形に対する台形度とを含む（５）項に記載の要素ゆがみ−解析誤差間関係推定方法。
【００３８】
この方法によれば、要素のゆがみが、少なくとも要素のアスペクト比と斜傾度と台形度とを含んで互いに独立した複数の幾何学変数によって表現される。
【００３９】
本項における「アスペクト比」は、要素の実形状の縦方向寸法と横方向寸法との比を意味し、種々の計算式を採用し得る。また、「斜傾度」は、要素の実形状の、基本形状に対する傾斜の度合いを意味し、種々の計算式を採用し得る。また、「台形度」は、要素の実形状の、基本形状に対する広がりの度合いを意味し、種々の定義を採用し得る。
【００４０】
本項における「要素の実形状」は、例えば、対象物の形状が分割されて取得された各要素の形状を意味する。
【００４１】
（７）  前記斜傾度が、前記要素の実形状の、前記基本長方形に対するスキュー量を用いて定義され、
前記台形度が、前記要素の実形状の、前記基本長方形に対するテーパ量を用いて定義される（６）項に記載の要素ゆがみ−解析誤差間関係推定方法。
【００４２】
後に理論的根拠を詳述するが、本発明者らは、少なくとも要素のアスペクト比と斜傾度と台形度とを含んで互いに線形独立した複数の幾何学変数によって要素のゆがみを表現するためには、斜傾度を、要素の実形状の、基本長方形に対するスキュー量を用いて定義するとともに、台形度を、要素の実形状の、基本長方形に対するテーパ量を用いて定義することが重要であることに気が付いた。
【００４３】
この知見に基づき、本項に係る方法が提案されたのである。
【００４４】
本項における「スキュー量」は、例えば、基本長方形をせん断変形させて要素の実形状を再現することを想定した場合のそのせん断変形量として定義することが可能である。
【００４５】
（８）  前記複数の幾何学変数が、さらに、前記要素の実形状の、前記基本長方形に対する反り度またはゆがみ度を含む（６）または（７）項に記載の要素ゆがみ−解析誤差間関係推定方法。
【００４６】
この方法によれば、要素のゆがみが、その要素の反りまたはゆがみをも考慮して表現することが可能となり、よって、反りまたはゆがみを考慮しないで表現する場合に比較し、その表現精度を容易に向上させることが可能となる。
【００４７】
（９）  前記推定工程が、前記複数の幾何学変数の各々に関連付けて前記解析誤差を個別解析誤差として推定し、それら複数の幾何学変数に関連付けてそれぞれ推定された複数の個別解析誤差に基づき、前記ゆがみの全体に起因する前記解析誤差を総合解析誤差として推定するものである（５）ないし（８）項のいずれかに記載の要素ゆがみ−解析誤差間関係推定方法。
【００４８】
本発明者らの研究により、互いに独立した複数の幾何学変数によって要素のゆがみを表現する場合には、ゆがみに起因した解析誤差を、各幾何学変数に関連付けて推定可能であることが判明した。
【００４９】
この知見に基づき、本項に係る方法においては、互いに独立した複数の幾何学変数の各々に関連付けて解析誤差が個別解析誤差として推定され、それら複数の幾何学変数に関連付けてそれぞれ推定された複数の個別解析誤差に基づき、要素のゆがみの全体に起因する解析誤差が総合解析誤差として推定される。
【００５０】
（１０）  前記推定工程が、前記推定された複数の個別解析誤差に基づき、前記解析結果に対する幾何学変数間の相乗効果を考慮して前記総合解析誤差を推定するものである（９）項に記載の要素ゆがみ−解析誤差間関係推定方法。
【００５１】
本発明者らの研究により、要素のゆがみのうち各幾何学変数により表現される各成分に起因した複数の個別解析誤差の単純な合計値より、ゆがみの全体に起因する総合解析誤差の方が大きい場合があることが判明した。この理由としては、いくつかの幾何学変数間の相乗効果が解析結果に対して与えられることが考えられる。
【００５２】
この知見に基づき、本項に係る方法が提案されたのである。
【００５３】
（１１）  対象物の形状を近似的に表現する複数の要素の組合体を用いて前記対象物の変形を有限要素法によって解析した場合にその解析結果に生じる解析誤差をコンピュータを用いて推定する解析誤差推定方法であって、
前記コンピュータのユーザの変形モード特定指令に従い、各要素ごとに、前記対象物の変形時に各要素に発生することが予想される変形モードを引張りと曲げとせん断とねじりとのうちの少なくとも一つとして特定する変形モード特定工程と、
前記変形モードの特定後、前記複数の要素のうちの少なくとも一つにつき、各要素の幾何学的なゆがみに起因した前記解析誤差を定量的に表示する表示工程と
を含む解析誤差推定方法。
【００５４】
この方法によれば、変形モードの特定後、複数の要素のうちの少なくとも一つにつき、各要素の幾何学的なゆがみに起因した解析誤差が定量的に表示される。
【００５５】
その定量的表示は、例えば、解析誤差の数値を直接に表示する方式としたり、図形、グラフ、色彩等を用いてグラフィカルに表示する方式とすることが可能である。
【００５６】
（１２）  前記組立体に対して有限要素法による数値解析が実施されるのに先立って実行される（１１）項に記載の解析誤差推定方法。
【００５７】
この方法は、対象物を互いに共同して表現する複数の要素の組立体に対して有限要素法による数値解析が実施されるのに先立って実行される。したがって、この方法によれば、要素ゆがみに起因した解析誤差を得ることを目的として有限要素法による数値解析を複数の要素の全体に対して実施せずに済み、解析誤差を取得するのに必要な時間を短縮することが容易となる。
【００５８】
（１３）  さらに、各要素ごとに、前記特定された変形モードに基づき、かつ、要素の幾何学的なゆがみとそれに起因して有限要素法による解析結果に生じる解析誤差との間に予め推定された要素ゆがみ−解析誤差間関係に従って前記解析誤差を前記コンピュータを用いて推定する推定工程を含み、
前記表示工程が、その推定された解析誤差を表示するものである（１１）または（１２）項に記載の解析誤差推定方法。
【００５９】
（１４）  前記要素ゆがみ−解析誤差間関係が、（１）ないし（１０）項のいずれかに記載の要素ゆがみ−解析誤差間関係推定方法により推定されたものである（１３）項に記載の解析誤差推定方法。
【００６０】
（１５）  前記ゆがみが、前記要素の幾何学的特徴を互いに独立して表現する複数の幾何学変数によって表現され、
前記推定工程が、それら複数の幾何学変数の各々に関連付けて前記解析誤差を個別解析誤差として推定する個別解析誤差推定工程を含み、
前記表示工程が、前記少なくとも一つの要素の各々に関連付けて、前記複数の幾何学変数の各々について、前記推定された個別解析誤差を表示するものである（１３）または（１４）項に記載の解析誤差推定方法。
【００６１】
（１６）  前記推定工程が、前記少なくとも一つの要素の各々に関連付けて、前記ゆがみの全体に起因する前記解析誤差を総合解析誤差として推定する総合解析誤差推定工程を含み、
前記表示工程が、前記少なくとも一つの要素の各々に関連付けて、前記推定された総合解析誤差を表示するものである（１３）ないし（１５）項のいずれかに記載の解析誤差推定方法。
【００６２】
（１７）  前記表示工程が、前記複数の要素を前記推定された総合解析誤差が許容値以上であるものとそうではないものとに区別して前記コンピュータのユーザに認識させるための表示を行う区別表示工程を含む（１６）項に記載の解析誤差推定方法。
【００６３】
（１８）  対象物の形状を近似的に表現する複数の要素の組合体を用いて前記対象物の変形を有限要素法によってコンピュータを用いて解析する数値解析方法であって、
前記対象物の形状を前記複数の要素に前記コンピュータを用いて分割する要素分割工程と、
前記複数の要素への分割後、前記コンピュータのユーザの変形モード特定指令に従い、各要素ごとに、前記対象物の変形時に各要素に発生することが予想される変形モードを引張りと曲げとせん断とねじりとのうちの少なくとも一つとして特定する変形モード特定工程と、
前記変形モードの特定に応答し、前記複数の要素のうちの少なくとも一つにつき、各要素の幾何学的なゆがみに起因した前記解析誤差を定量的に表示する表示工程と、
前記組合体を用いて前記対象物の変形を有限要素法によって前記コンピュータによって解析する解析工程と
を含む数値解析方法。
【００６４】
（１９）  （１）ないし（１８）項のいずれかに記載の方法を実施するためにコンピュータにより実行されるプログラム。
【００６５】
このプログラムがコンピュータにより実行されれば、前記（１）ないし（１８）項のいずれかに係る方法におけると基本的に同じ原理に従い、同様な効果が実現され得る。
【００６６】
このプログラムは、それの機能を果たすためにコンピュータにより実行される指令の組合せのみならず、各指令に従って処理されるファイルやデータをも含むように解釈することが可能である。
【００６７】
また、このプログラムは、それ単独でコンピュータに実行されることにより、所期の目的を達するものとしたり、他のプログラムと共にコンピュータに実行されることにより、所期の目的を達するものとすることができる。後者の場合、本項に係るプログラムは、データを主体とするものとすることができる。
【００６８】
（２０）  （１９）項に記載のプログラムをコンピュータ読み取り可能に記録した記録媒体。
【００６９】
この記録媒体に記録されているプログラムがコンピュータにより実行されれば、前記（１）ないし（１８）項のいずれかに係る方法と同じ作用効果が実現され得る。
【００７０】
この記録媒体は種々な形式を採用可能であり、例えば、フレキシブル・ディスク等の磁気記録媒体、ＣＤ、ＣＤ−ＲＯＭ等の光記録媒体、ＭＯ等の光磁気記録媒体、ＲＯＭ等のアンリムーバブル・ストレージ等のいずれかを採用し得る。
【００７１】
（２１）  （１）ないし（１８）項のいずれかに記載の方法をコンピュータによって実施する装置。
【発明を実施するための最良の形態】
【００７２】
以下、本発明のさらに具体的に実施の形態の一つを図面に基づいて詳細に説明する。
【００７３】
本発明の一実施形態に従う数値解析方法は、解析対象である構造物の形状を近似的に表現する複数の要素の組合体を用いることにより、その構造物の変形を有限要素法によって解析する方法である。
【００７４】
この数値解析方法においては、概略的に説明すれば、構造物が複数の要素に自動的に分割され、その分割された各要素ごとに、要素のゆがみに起因して将来の解析結果に生ずるであろう解析誤差が推定される。要素のゆがみは、有限要素解析のための後述の要素分割に起因して発生する。その推定された解析誤差が許容値を超える要素については、そのことが解析者に告知される。その結果、解析者は、解析誤差が許容値を超える要素、すなわち、解析誤差との関係においてゆがみが過大な要素が存在しないように構造物の要素を再分割することが可能となる。
【００７５】
この数値解析方法においては、さらに、解析者からの解析指令に応答し、形状が確定された複数の要素を用いることにより、構造物の変形が有限要素法によって解析される。
【００７６】
本実施形態においては、解析誤差の推定が、本発明者らにより案出されたいくつかの新規な知見を前提にして行われる。そのような知見としては次のものがある。
【００７７】
（１）要素のゆがみの定義
【００７８】
任意の要素の幾何学的形状は、互いに線形独立な複数のゆがみモード（幾何学変数）によって表現できる。換言すれば、要素の形状は、複数のゆがみモードの線形和として表現できる。
【００７９】
ここに、複数のゆがみモードは、要素の反りを無視することとすれば、要素のアスペクト比に関するモードと、２直交方向のそれぞれにおける斜傾性（スキュー量）に関するモードと、２直交方向のそれぞれにおける台形性（テーパ量）に関するモードとを含み、かつ、それら以外のモードは含まない。
【００８０】
したがって、要素のゆがみは、アスペクト比、斜傾性および台形性のみによって定義できる。
【００８１】
（２）要素の変位の定義
【００８２】
要素の変位は、互いに線形独立な複数の変形モードによって表現できる。換言すれば、要素の変位は、複数の変形モードの線形和として表現できる。
【００８３】
ここに、複数の変形モードは、２直交方向のそれぞれにおける平行移動モードと、２直交方向のそれぞれにおける引張りモードと、２直交方向のそれぞれにおけるせん断モード（単純せん断モード）と、２直交方向のそれぞれにおける曲げモード（単純曲げモード）とを含み、かつ、それら以外のモードは含まない。
【００８４】
（３）要素のゆがみとそれに起因する解析誤差との相関関係
【００８５】
要素のゆがみモードを特定すれば、そのゆがみモードとそれに起因する解析誤差との間には、要素の変形モードに依存した一定の関係が成立する。要素のゆがみに起因する解析誤差は、要素の変形モードごとに互いに分離して解析できる。
【００８６】
以下、それらの知見の理論的根拠を説明する。
【００８７】
曲線や曲面の境界を有する領域を直線や平面の境界を有する要素で理想化することはできない。そのため、曲線や曲面を有する要素が開発された。この要素は、パラメータ上においては直線要素と等価であるため、アイソパラメトリック要素と称される。
【００８８】
したがって、複雑な構造物を解析するためにアイソパラメトリック要素が使用される場合がある。この場合、構造物はしばしば、不規則な四辺形を有する面要素（パネル要素）に分割され、ゆがんだメッシュ（要素）を伴う。任意の四辺形を実用的に処理するために、そのゆがんだ要素は、自然座標（ξ，η）を用いて定義された双一次形状関数の導入により、基本の正方形要素に写像される。アイソパラメトリック要素内においては、内部変位分布も同じ双一次関数によって近似させられ、不完全な線形ひずみ変分が、物体座標（ｘ，ｙ）に関する変位の微分を採用することによって取得される。そして、対応する応力および要素剛性マトリクスが節点変位ベクトルによって計算され得る。
【００８９】
有限要素の種類を問わず、近似精度は、要素のアスペクト比が１から外れるにつれて減少するとともに、要素の形状が長方形からくずれるにつれて減少する。不規則なメッシュによる精度悪化の理由はヤコビのマトリクスにおいて見い出される。そのマトリクスは、要素剛性マトリクスの計算時、ひずみの評価における微分演算子の座標変換、または形状関数の微分における要素形状を直接に反映する。以下、特にひずみ分布において要素のゆがみが近似精度に及ぼす影響をアイソパラメトリックな四辺形要素の分野に関して検討する。
【００９０】
要素ゆがみ変数は、任意の要素の形状を、アクペクト比(縦横比、aspect ratio)に関連する斜傾性(skew)のゆがみと台形性(trapezoid)のゆがみとの組合せとして表現することによって定義される。それら斜傾性のゆがみと台形性のゆがみとはいずれも、アスペクト比(縦横比、aspect ratio)を意味する基本の長方形(rectangle)を基準に定義される。その結果、ヤコビのマトリクスは要素ゆがみ変数によって記述される。いずれの項がアイソパラメトリックな四辺形要素に付随するかを検討するために、ひずみ成分がヤコビの逆行列に関して展開される。
【００９１】
要素のゆがみは次のように定義される。すなわち、アイソパラメトリック要素の形状（ジオメトリ、幾何学的特徴）を、パラメトリックである自然座標（以下、「パラメトリック座標」という。）ξおよびηであって要素の領域内において図１の式（１．１）として一義的に定義されるものに関する双一次多項式を用いて仮定するのである。
【００９２】
この式においてＮｉ（ξ，η）は、図１の式（１．２）により定義された双一次形状関数である。この関数は、二次元要素について物体座標とパラメトリック座標との間の座標関係を定立する。前記式（１．１）における各係数は、同図の式（１．３）により表される。この式においてｈ０ないしｈ３は、直交ベクトルであり、それの値は、図２の式（１．４）により表される。前記式（１．３）における節点座標ベクトルｘ，ｙは、同図の式（１．５）により表される。
【００９３】
上述の概念に従い、微分演算子の変換において使用されるヤコビのマトリクスは、図２の式（１．６）として表現され得る。
【００９４】
ヤコビのマトリクスは、任意の四辺形要素を辺の長さが２である正方形領域に写像するものであるため、ヤコビのマトリクスの各成分は、要素の形状に直接に関連する。この相関を具体的に説明するために、次に、長方形要素、平行四辺形要素および台形要素についてのヤコビのマトリクスを検討する。
【００９５】
図３に示すように、局所座標系（ｘ，ｙ）において定義された任意の四辺形要素の形状を、基本の長方形(rectangle)のアクペクト比と共に、次の４つの変形モード、すなわち、ｘ方向およびｙ方向における斜傾性(skew)ならびにｘ方向およびｙ方向における台形性(trapezoid)に分類する。同図には、要素のゆがみのいくつかの種類が示されている。
【００９６】
図３に与えられた次元に従えば、基本の長方形の形状、平行四辺形の形状および台形の形状は、直交ベクトルｈ１，ｈ２およびｈ３を用いて図４の式（１．７）により表記される。この式において最上段の左辺は基本の長方形の節点座標を表し、２段目および３段目の両左辺は平行四辺形の節点座標を表し、４段目および５段目の両左辺は台形の節点座標を表している。
【００９７】
さらに、要素のアクペクト比は、基本長方形の縦辺と横辺との長さの比である。さらにまた、任意の四辺形要素のゆがみの程度、すなわち、ｘ方向およびｙ方向における斜傾性および台形性の程度は、その基本長方形からの差として認識することが可能である。その差は、図４の式（１．８）により記述され、これは同図の式（１．９）と等しい。
【００９８】
直交ベクトルｈ１，ｈ２およびｈ３が線形独立であるため、各ゆがみが近似精度に与える影響は互いに分離して検討し得る。ここに、直交条件ｈｉＴｈｊ＝４δｉｊ（δｉｊ：クロネッカのデルタ）により、上記節点座標を前記式（１．３）に代入すると、図５の式（１．１１）を得る。この式が示していることは、係数ａｉｘおよびａｉｙ（ｉ＝１，２，３）が実際に正確に基本長方形の形状およびゆがみ変数であり、一方、４つの座標の平均値ａ０ｘ，ａ０ｙが要素の図心の位置を示しているということである。ゆがみ変数を用いることにより、前記式（１．１）における内部座標は図５の式（１．１２）としても記述される。
【００９９】
この式における（ｘ０，ｙ０）は、要素の図心、すなわち、局所座標系の原点である。さらに、その式における幾何学変数ｈｘ，ｈｙ，α１，α２，α３，α４は図５の式（１．１３）により表される。
【０１００】
それら幾何学変数を前記式（１．１０）に代入すると、図５の式（１．１４）により表される関係が得られる。
【０１０１】
ここに、幾何学変数を実際に得るために任意の四辺形の一例を考える。それに先立ち、ある要素において一義的に対応させられた局所座標を定義することが必要である。図６には、その局所座標の設定が示されている。
【０１０２】
まず、四辺形が反っている場合、すなわち、４つの角の節点が同一平面上に存在しない場合には、図６に示すように、その四辺形は平面上に、その平面の複数の節点と四辺形の複数の節点との距離が互いに等しくなるように投影される。これにより、もとの四辺形要素（同図において破線で示す）に対応する平面要素（同図において実線で示す）が得られる。もとの四辺形の複数の節点の、上記平面要素に対して直角な方向における座標は、ｈｗｈ３である。ここに、ｈｗは投影高さを意味する。要素の反り(warp)は、構造物の形状の近似精度を低下させる主な要因であるが、今回は、上記平面の方向におけるゆがみのみを考慮する。
【０１０３】
まず、局所座標系の原点が上記平面要素の図心の位置に設定される。このことは図７の式（１．１５）として表される。この式におけるＸおよびＹは、４つの角の節点の全体座標である。
【０１０４】
次に、前記平面要素の各辺(side)が二等分され、対辺上の中点同士をつなぐ２つのベクトルＰ２４，Ｐ３１が図７の式（１．１６）によって定められる。ベクトルＰ２４は、局所的なｘ軸の方向として設定され、ｙ軸は、そのｘ軸に対して直角に定義される。このことは同図の式（１．１７）として表される。この式におけるＰ２４（ただし、特殊な記号が付されたベクトルＰ２４)は、同図の式（１．１８）により表される。
【０１０５】
局所座標{ｘｉ，ｙｉ}は、座標変換マトリクスＴ＝［ｅｘ，ｅｙ］を用いて、要素の図心に関する位置ベクトルを変換することによって得られる。このことが図８の式（１．１９）により表されている。
【０１０６】
最後に、前記平面要素に含まれるゆがみの程度が前記式（１．１３）によって定められる。
【０１０７】
任意の形状を有する四辺形要素であって、それの複数の節点(Node 1, Node 2, Node 3, and Node 4)の全体座標Ｘｉ，Ｙｉが図９の表に列挙されているものを例にとって考える。局所座標系の原点(origin)の位置は、前記式（１．１５）に示すように４つの節点座標の平均値をとることによって定められ、その単位軸ベクトルｅｘ，ｅｙは、前記式（１．１７）を計算することにより、図９の表に示すように得られる。その後、座標変換マトリクスを適用することにより、節点座標ｘｉ，ｙｉが図１０の表に示すものとして定められる。任意の四辺形およびそれの局所座標系ｘ−ｙは図１１に示されている。同図においてＸ−Ｙは全体座標系である。
【０１０８】
前記式（１．１３）を利用することにより、その四辺形の幾何学変数が次のように得られる。
【０１０９】
ｈｘ＝（１／４）ｈ１Ｔｘ＝  ２．１０３３
ｈｙ＝（１／４）ｈ２Ｔｙ＝  １．５７３８
アスペクト比＝ｈｘ／ｈｙ＝１．３３６４
α１＝（１／４）ｈ２Ｔｘ＝−０．４０４８
α２＝（１／４）ｈ１Ｔｙ＝−０．３０２９
α３＝（１／４）ｈ３Ｔｘ＝−１．１４９１
α４＝（１／４）ｈ３Ｔｙ＝−０．０８７０９
【０１１０】
さらに、前記式（１．１０）から、節点座標が、直交ベクトルｈ１，ｈ２およびｈ３の線形結合であってそれらの係数がここで得られた幾何学変数であるものとして、表現されることが証明され得る。さらに、ｘ方向における台形性ゆがみが、この特定の四辺形要素に含まれるタイプの支配的なゆがみであることも認識される。前述のように、任意の四辺形は、基本の長方形を基準として、４つのゆがみのタイプ、すなわち、図１２に示すように、ｘ方向およびｙ方向における斜傾性と、ｘ方向およびｙ方向における台形性とにより構成される。同図においては、破線で元のゆがんだ要素が示され、一方、実線で、そのゆがんだ要素に含まれる斜傾性または台形性の程度がゆがみの定義に従って示されている。
【０１１１】
以上説明した要素ゆがみのモードに従い、前記式（１．６）により表される四辺形要素のヤコビのマトリクスが、幾何学変数を用いて、図１３の式（１．２０）により表現され得る。
【０１１２】
これに対して、ヤコビのマトリクスは、図１３の式（１．２１）に示すように、３つの部分、すなわち、直交定数項、非直交定数項および線形項に分解し得る。直交成分を有する１番目の定数項は基本の四辺形を表現し、非直交成分を有する２番目の定数項は要素の斜傾性を示し、３番目の線形項は台形性のゆがみを表している。図１４には、ゆがみのモードとヤコビのマトリクスとが示されている。
【０１１３】
換言すれば、任意の要素の形状は、基本の長方形に付加される複数のゆがみ量の合計により構成されるため、ヤコビのマトリクスもまた、図１４に示すように、長方形の項と、定数かつ非直交の斜傾性の項および線形かつ台形性の項とによって構成される。ヤコビのマトリクスは、要素が長方形または平行四辺形である場合には定数であり、また、台形である場合には線形となることが観察される。
【０１１４】
ヤコビの行列式であって、座標変換が単位領域の面積を変化させる程度を表すものは、図１５の式（１．２２）に示す線形関数である。ここに、Ｊ０，Ｊ１およびＪ２は同図の式（１．２３）により表される。
【０１１５】
さらに、幾何学的な量に関するヤコビの行列式は、図１５の式（１．２４）により表現される。
【０１１６】
明らかに、ヤコビの行列式もまた、長方形および平行四辺形については定数である。さらに、ヤコビの逆行列は、図１５の式（１．２５）となる。カーテシアン座標に関する微分は、ヤコビの逆行列を用いて、パラメトリック座標に関して座標変換することによって得られる。このことが図１６の式（１．２６）により表されている。
【０１１７】
数値解析のためには、現在問題にしている有理関数よりむしろ多項式を用いて上記微分演算子を表現するほうが望ましい。そのため、ヤコビの行列式の逆数１／Ｊ（ξ，η）は、要素の図心についてのテイラー級数によって展開することができる。このことが図１６の式（１．２７）により表されている。この式において１次項まで採用すれば、ヤコビの逆行列は、同図の式（１．２８）によって近似される。
【０１１８】
これに対し、テイラーの定理はヤコビの逆行列それ自体に適用することが可能であり、その結果、図１６の式（１．２９）により表される近似式が得られる。ヤコビの逆行列をそれらの近似式によって置き換えると、形状関数の微分が多項形式で記述され得る。
【０１１９】
アイソパラメトリック要素の面内変位分布は、座標近似のために利用されたものと同じ双一次形状関数Ｎｉ（ξ，η）によって仮定される。アイソパラメトリック要素の内部変位関数は、図１７の式（１．３０）により表される。ここに、多項式の係数は同図の式（１．３１）によって与えられる。この式における節点変位ベクトルｕ，ｖは同図の式（１．３２）により表される。
【０１２０】
要素の形状が、基本の長方形に付加される斜傾性および台形性のゆがみの組合せにより表現される方法と同様に、アイソパラメトリック要素の変形は、図１８に示すように、ｘ方向およびｙ方向における４つのモード（４つの基本変形モード）、すなわち、平行移動モード(translation mode)（応力が０）、引張りモード(stretch mode)（垂直応力が一定）、せん断モード(shear mode)（せん断応力が一定）および曲げモード(bending mode)（垂直応力が線形）により構成される。剛体の運動は、平行移動と回転とを含んでいる。平行移動は、ｘ方向およびｙ方向の２つの平行移動モードにより表現される。回転は、ｘ方向におけるせん断モードとｙ方向におけるせん断モードとの差により表現される。
【０１２１】
前記式（１．２６）により表される微分演算子の座標変換を、自然座標に関する変位の微分に適用すると、図１９の式（１．３３）が得られる。ひずみ場は、同図の式（１．３４）により得られる。
【０１２２】
微分演算子についての上述の座標関係から、ｊ−１１，ｊ−１２，ｊ−２１およびｊ−２２に含まれる非直交定数項および線形関数であって、要素の斜傾性および台形性のゆがみに付随するものは、過大なひずみエネルギの発生に寄与することが分かる。その過大なひずみエネルギの発生は、要素の硬い反応を生じさせる。
【０１２３】
各ひずみ成分は実際、図１９の式（１．３５）により表される複数の関数を得るために展開される。それら複数の関数は、マトリクス表記によれば、図２０の式（１．３６）として記述され得る。ここにＧは、同図の式（１．３７）により表される。
【０１２４】
ひずみはさらに、図２１の式（１．３８）、式（１．３９）および式（１．４０）に示すように、ゆがみ変数αｉ（ｉ＝１，２，３，４）に関して記述されるように再構成される。
【０１２５】
各ひずみ成分については、１番目の部分が、基本の長方形に関連し、２番目の項は、要素が斜傾性を有するときに現れ、３番目の部分は、台形性のゆがみに対応する。例えば、要素が置かれている場が単純曲げである場合には、単純曲げモードがせん断ひずみを生じさせないことが既に知られており、したがって、この場合には、せん断成分におけるｕ３ｘおよびｕ３ｙの項は無視すべきである。その結果、要素が長方形である場合には、せん断ひずみが一定となり、不要なせん断ひずみエネルギが消去される。しかし、図１４から容易に分かるように、要素がゆがむ程度が大きいほど、ゆがみ変数（基本ゆがみ変数、幾何学変数）α１ないしα４が大きくなる。したがって、要素内における過大なひずみエネルギが発生する。換言すれば、要素に蓄えられるひずみエネルギは、その面積が変化しなくても、その要素のゆがみ程度の増加につれて増加する。
【０１２６】
各ひずみ項を詳細に検討して気が付くことは、ひずみεｘが、ｙ方向における幾何学変数α２，α４を、ｘ方向における変形モードと共に含み、一方、変形モードベクトル（節点変位ベクトル）ｕ，ｖ（図１７の式（１．３１）および式（１．３２）ならびに図１８参照）が、対応する幾何学モードベクトル（直交ベクトル）ｈ０ないしｈ３（図２の式（１．４）参照）に対して直交するということである。同様にして、ひずみεｙが、ｘ方向における幾何学変数α１，α３を、ｙ方向における変形モードと共に含む。他方、ひずみγｘｙにおける幾何学変数α１ないしα４が、同じ方向における変形モードに関連する。
【０１２７】
また、長方形の幾何学変数および平行四辺形の幾何学変数は曲げ変形モード（線形せん断）に関連し、また、台形の幾何学変数は線形変形モード（定常せん断）に対応する。換言すれば、一方向における要素の斜傾性が他の方向における曲げ変形エネルギを過大にし、一方、一方向における台形性ゆがみが他方向における引張り変形およびせん断変形の応力を過大にする。結局、現実のひずみ場が一定である場合であっても、台形性ゆがみが線形のひずみ分布を生じさせるため、台形性ゆがみを有するメッシュを避けることが望ましい。
【０１２８】
長方形要素の場合、ひずみ分布は、図２２の式（１．４１）による表される線形関数に単純化される。
【０１２９】
以上、要素のゆがみを面要素に着目して検討したが、シェル型構造物を解析する場合には、平板の曲げ成分も考慮しなければならない。この場合には、面内変形が常に面外変形を伴う。
【０１３０】
以下、ミンドリン−ライスナ型の板曲げ要素（プレート要素、ミンドリン板）の分野を説明する。
【０１３１】
まず、平板のモデル化を説明すると、全体基準フレーム｛Ｘ，Ｙ，Ｚ｝内の任意の点Ｐ（Ｘ，Ｙ，Ｚ）における平板の変位｛ｕｘ，ｕｙ，ｕｚ｝は、局所座標系｛ｘ，ｙ，ｚ｝内で表現される。この局所座標系のｘｙ平面は、各板要素の中立面上に一義的に位置決めされ、一方、ｚ軸は、その中立面に対して直角な方向に設定される。
【０１３２】
変位成分は、面内の部分と面外の部分とに分解される。２次元の面内成分は、ｘ−ｙ平面、すなわち、板の中立面上に定義され、また、横成分（面外成分）は、ｚの関数として表現される。例えば、２つの部分より成る変位のｘ成分は、図２３の式（２．１）および式（２．２）のように考える。その後、３次元場における変位が、漸近的な展開式を用いて同図の式（２．３）のように表現される。
【０１３３】
せん断効果が厚板解析に考慮されるように、高次の板要素を作る３次項まで選択することが望ましい場合がある。十分に薄い板については、変位の面内成分ｕｘおよびｕｙが、厚さ方向において線形変化するように近似させられる。垂直変位ｕｚは、板の厚さにかかわらず一定となるように仮定される。それらのことが図２３の式（２．４）により表されている。
【０１３４】
板の変位場は、図２４に示すように、５つの自由度｛ｕ，ｖ，ｗ，θｘ，θｙ｝、すなわち、３種類の平行移動および２種類の回転によって表現される。それら一般化された変位（回転も含む）は、中立面上に定義されるとともに、ｘおよびｙの関数として表現される。
【０１３５】
上述の変位近似に従い、板のひずみは、図２５の式（２．５）により表現されるものとなる。これは、マトリクス表記によれば、同図の式（２．６）により表現される。
【０１３６】
さらに、板の応力σｚ＝０である状態を仮定すると、応力は、図２６の式（２．７）の構成方程式により得られる。ここに、κｘおよびκｙは、せん断を補正する係数である。
【０１３７】
一般化された５つの変位のうち、ｕおよびｖは、面内変形であり、ｗ，θｘおよびθｙが、板の曲げ変形に対応する。このことが図２６の式（２．８）により表されている。
【０１３８】
その結果、ひずみは、面内ひずみ成分と横ひずみ成分とに分解することができる。ひずみのうちの面内部分は、板の曲げ変形｛ｗ，θｘおよびθｙ｝と共に面の変形｛ｕ，ｖ｝によって発生する。このことが図２６の式（２．９）により表されている。
【０１３９】
一方、横せん断ひずみは、横変形（面外変形）のみに関連する。このことが図２７の式（２．１０）により表されている。
【０１４０】
板のひずみエネルギは、面内の部分と面外の部分とに分解される。このことが図２７の式（２．１１）により表されている。この式における材料定数マトリクスＤｐ（面内）およびＤｓ（面外）が同図の式（２．１２）により表されている。
【０１４１】
ここで、Ａを中立面(neutral surface)の領域を表すものとし、ｈを板の厚さを表すものとし、Ｓを領域Ａの面積を表すものとし、そして、板の領域を、Ω＝Ａｘ［−ｈ／２，＋ｈ／２］とする。板の厚さｈに関して積分すると、横変形（面外変形）に起因するひずみエネルギは、図２７の式（２．１３）のように計算される。
【０１４２】
この式における２番目の項から観察できることは、板の厚さｈが薄くなるにつれて、せん断によるひずみエネルギが支配的となり、シェア・ロッキングが起こる。この現象は、薄板の実際の挙動とは矛盾する。この不要なせん断効果は、厚さ方向における横せん断ひずみの導入に起因し、そのひずみは、薄い曲げ板において消失する。換言すれば、板の厚さｈが０に近づくにつれて、横せん断ひずみγｓは０に近づかなければならず、薄い板曲げ要素における横せん断効果を考慮することは本質的に、せん断ひずみが０である状態を強制する拘束条件として理解され得る。
【０１４３】
なお付言すれば、板曲げ要素（ミンドリン板）に関しても、面要素と同様に、変位を面外の変形モードで表現し、板曲げ要素から生じるひずみも、同様にして、ゆがみモード（ゆがみ変数、幾何学変数）と変形モードとによって展開できる。
【０１４４】
次に、要素のゆがみに起因したひずみエネルギの変化をミンドリン型の要素に関して説明する。
【０１４５】
図３５に示すように、形状が変化する１つの四辺形が大きな板の中央に位置させられている。その板は、一様引張り(uniform tension)状態、単純曲げ(pure bending)状態および単純せん断(pure shear)状態で負荷を受ける。図３５は、面内変形を示すが、図２８は、面外変形を示す。図２８においては、板の角部が単純に支持された状態で、板が、単純曲げ状態(pure bending)、ねじり(torsion)状態、および、中央位置に一点で力を受ける状態(one-point center force)で負荷を受ける。
【０１４６】
ひずみエネルギ誤差Ｅｅは、図２９の式（３．１）により定義される。ここに、Ｕｅ（ｕ）は、ひずみエネルギの解析解である。これは、図２９における式（３．４）として表されるひずみエネルギ式によって定義される。ひずみエネルギの解析解Ｕｅ（ｕ）は、具体的には、応力場の厳密解をそのひずみエネルギ式に代入するとともに、それを要素の領域内にわたって積分することにより計算されるひずみエネルギである。
【０１４７】
例えば、梁状の面構造物が面内において単純曲げを受ける場合には、応力の分布が図２９の式（３．２）により表され、ひずみ場が同図の式（３．３）により表される。したがって、そのひずみ場を同図の式（３．４）により表されるひずみエネルギ式に代入し、要素を積分領域とすれば、ひずみエネルギＵｅ（ｕ）が求められる。その式（３．４）におけるＤは材料定数マトリクスを意味する。
【０１４８】
これに対し、Ｕｅ（ｕｈ）は、有限要素法により計算されるひずみエネルギの近似解である。具体的には、これは、有限要素法の剛性マトリクスと図３３の式（３．１０）におけるような節点変位の厳密解とを用いて求められるひずみエネルギである。例えば、片持ち梁の単純曲げの場合であれば、節点変位が、節点座標を図３０の式（３．５）に代入することによって求められる。この節点変位の厳密解を用いることにより、有限要素近似による要素のひずみエネルギが、図３３の式（３．１０）を適用することによって求められる。
【０１４９】
以上、変位(displacement)ｕによるひずみエネルギの計算を説明したが、変位ｖによるひずみエネルギは、図３１の式（３．８）により表される。この式（３．８）におけるＤは材料定数マトリクスを意味する。
【０１５０】
なお、表記法について説明すれば、Ｕｅ（ｕ）およびＵｅ（ｕｈ）における各添字ｅは、１要素についてのひずみエネルギであることを意味する。
【０１５１】
ひずみエネルギの近似解を計算するために必要な節点変位は、３つの負荷条件(pure bending, torsion, and 1-point loading)について図３２の表に列挙された解析的な変位関数、または３つの負荷条件(uniform tension, pure bending, and pure shear)について図４４の表に列挙された解析的な変位関数を用いて計算され得る。ここでは、図３３の式（３．９）により表される基底関数φｉｊ（ｘ，ｙ）が点負荷問題(one-point loading problem)の解析解を求めるために利用される。
【０１５２】
なお付言すれば、式（３．２）と式（３．３）は、図３５に示す面内の単純曲げの場合の応力場を表しており、そのときの厳密解が式（３．５）で表されている。すなわち、式（３．５）は、図３５に示す面内曲げ（引張り、曲げおよびせん断）の場合の厳密解を表しているのであり、また、式（３．５）に要素の節点座標値を代入し、それを式（１．３１）に入れると、式（３．６）が得られる。ただし、このことは、例えば、梁状の構造物にモーメントが与えられた場合の固定端の要素に関して成立する。
【０１５３】
ひずみエネルギの近似解は、図３３の式（３．１０）により表されるように、要素剛性マトリクス（ゆがんだ要素の剛性マトリクス）Ｋｅと共に、真の節点変位ｄ（近似解）から得られる。真の節点変位ｄは、四辺形の節点座標を前述の各負荷条件について与えられた解析的な変位関数（図３２参照）に代入することにより計算される。
【０１５４】
特に応力仮定要素である場合には、ひずみエネルギ誤差Ｅｅの計算値から、アスペクト比または斜傾性ゆがみの変化はひずみエネルギの近似解に影響しないことが分かる。台形性ゆがみだけが要素の近似精度を低下させる。図３４には、台形性ゆがみに起因するひずみエネルギ誤差Ｅｅが単純せん断、単純曲げ、点負荷およびねじりという各負荷条件ごとにプロットされている。図３５には、台形性ゆがみを有する要素が引張り、単純曲げおよびせん断のそれぞれによって変形する様子が示されている。図３４は、横変形（面外変形）の実験結果を示すが、図４５は、一様引張り(uniform tension)、単純曲げ(pure bending または in-plane bending)および単純せん断(pure shear または in-plane shear)についての面内変形の実験結果を示す。
【０１５５】
ひずみエネルギについての厳密解は、完全に長方形である要素を用いて得られる一方、台形性ゆがみが増加するにつれて、面内曲げ挙動が非常に硬くなる。曲げを受ける台形要素のその硬い反応は、その要素のアスペクト比が１より大きくなるときに強調される。この傾向は、図３５と同様に大きなアスペクト比を有する台形要素を設置することによって検証することが可能である。一様引張りおよび単純せん断の試験は、定ひずみパッチ試験であって、要素の収束を証明するために、定ひずみ状態を再現する能力を検証するものが満足させられることを示している。換言すれば、引張りおよびせん断の状態においてゆがんだ要素がひずみエネルギ誤差を発生させない場合には、その要素は、そのパッチ試験に合格する。ここで試験された応力仮定要素は、単純せん断状態においてひずみエネルギ誤差が、その要素の台形性の増加に対して一定に保たれないため、その要素は、いわゆる定ひずみパッチ試験に合格しない。
【０１５６】
台形性ゆがみが増加するにつれて、面内曲げ挙動は、せん断係数の修正なしでは非現実的となる。板の曲げほどではないが、ねじりが加わると、要素の剛性が極めて大きくなる。台形性ゆがみによって引き起こされるこの不適当な誤差の増加を改善するために、せん断係数を人為的に修正することが考えられる。横せん断ひずみは、小さな補正係数（制御変数）κｘおよびκｙを単に掛け算することによって変更される。このことが図３６の式（３．１１）により表されている。図３７には、台形性ゆがみに起因するひずみエネルギ誤差Ｅｅの、せん断係数の修正後の値が単純せん断、単純曲げ、点負荷およびねじりという各負荷条件ごとにプロットされている。
【０１５７】
ここに、制御変数κｘ，κｙは０．３１７０２と等しくなるように設定される。その結果、曲げ挙動が著しく改善される。一方、制御変数κｘ，κｙは曲げ効果とせん断効果との間のひずみエネルギバランスを制御するために設計されているから、ねじりは制御変数κｘ，κｙの影響をあまり多く受けない。
【０１５８】
以上、本実施形態に従う数値解析方法のうち、要素ゆがみに起因する解析誤差の推定の理論的根拠を詳細に説明したが、次に、この数値解析方法の具体的内容をそれを実施するコンピュータ・システムと共に詳細に説明する。
【０１５９】
図３８には、そのコンピュータ・システム１０が斜視図で示されている。このコンピュータ・システム１０は、コンピュータ１２と、入力装置１４としてのキーボード１６および図示しないマウス（ポインティング・デバイスの一例）と、出力装置２０としてのモニタ２２とを備えている。
【０１６０】
図３８に示すモニタ２２はＣＲＴ式であるが、ＬＣＤ式等、他の形式のモニタを採用することが可能である。また、同図に示すコンピュータ・システム１０は、その使用位置が固定される形式であるが、可搬式のコンピュータ・システムを採用することが可能である。
【０１６１】
図３９に示すように、コンピュータ１２は、よく知られているように、ＣＰＵ等、プロセッサ３０と、ＲＯＭ，ＲＡＭ等のメモリ３２とがバス３４により互いに接続されて構成されている。
【０１６２】
コンピュータ１２は、記憶装置４０と共に使用される。記憶装置４０は、ＣＤ−ＲＯＭ、ＭＯ等、記録媒体４２に対してデータの読み書きを行う。記録媒体４２には、この数値解析方法を実施するための複数のプログラムが予め記録されており、必要なプログラムがコンピュータ１２により記録媒体４２から読み出され、その読み出されたプログラムがコンピュータ１２により実行される。
【０１６３】
図４０には、この数値解析方法の内容がフローチャートで概念的に表されている。
【０１６４】
この数値解析方法においては、まず、ステップＳ１（以下、単に「Ｓ１」で表す。他のステップについても同じとする）において、構造物の形状を表すデータに基づき、その構造物の形状が近似的に複数の四辺形の要素（例えば、シェル要素）に自動的に分割される。この要素分割は、図示しない要素分割プログラムがコンピュータ１２により実行されることによって行われる。その分割された複数の要素をそれぞれ表すデータ（例えば、各要素の形状を定義するためのデータ）は、要素番号ｉに関連付けて、メモリ３２または記憶装置４０にストアされる。
【０１６５】
次に、Ｓ２において、その分割された各要素ごとに、その要素のゆがみに起因して、その要素を用いた有限要素法による解析結果に生ずることが予想される誤差が推定される。この解析誤差推定は、図４１にフローチャートで表される解析誤差推定プログラムがコンピュータ１２により実行されることによって行われる。この解析誤差推定プログラムは後に詳述する。
【０１６６】
続いて、Ｓ３において、その推定された解析誤差に基づき、ユーザにより、各要素ごとに、その修正が必要であるか否かが判断される。
【０１６７】
その必要がある場合には、Ｓ４において、ユーザにより、またはユーザの指示に従い、修正が必要である要素が修正される。この要素修正により、過大な解析誤差を生じさせる原因となるゆがみを有する要素の形状が、正方形に近づくように、すなわち、ゆがみが軽減されるように変更される。
【０１６８】
その後、Ｓ２に戻り、再度、解析誤差推定が自動的に行われる。それらＳ２ないしＳ４の実行が繰り返される結果、すべての要素についてゆがみの程度が許容限度以下となれば、Ｓ５において、形状が確定した複数の要素を用いることにより、有限要素法のもとに構造物の変形が解析される。この数値解析は、図示しない数値解析プログラムがコンピュータ１２により実行されることによって実行される。
【０１６９】
以上で、この数値解析法の一回の実行が終了する。
【０１７０】
なお付言すれば、図４０のフローチャートは、この数値解析方法のためにコンピュータ１２により実行される全体プログラムの内容を概念的に表すフローチャートであると考えることができる。
【０１７１】
図４１には、前記解析誤差推定プログラムの内容がフローチャートで概念的に表されている。
【０１７２】
この解析誤差推定プログラムにおいては、まず、Ｓ１０１において、要素番号ｉが１にセットされる。
【０１７３】
次に、Ｓ１０２において、今回の要素番号ｉを有する要素（以下、「今回の要素」という）の変形モードが特定される。この特定は、ユーザにより指定された変形モードを用いて行われる。ユーザの指令に従い、事前に、各要素に関連付けて変形モードがメモリ３２または記憶装置４０にストアされており、今回の要素に関連する変形モードがメモリ３２または記憶装置４０から読み出され、それにより、今回の要素の変形モードが特定されるのである。
【０１７４】
続いて、Ｓ１０３において、今回の要素を表すデータがメモリ３２または記憶装置４０から読み出される。さらに、その読み出されたデータに基づき、今回の要素のアスペクト比Ａが計算される。
【０１７５】
アスペクト比Ａは、従来からの定義に従い、例えば、次のようにして計算される。
【０１７６】
すなわち、まず、今回の要素について投影平面が、今回の要素の４つの節点から等距離にあるように想定される。次に、今回の要素に適当な固有の局所座標系が設定され、その一方、その座標系の原点が要素の図心に位置させられる。任意の形状を有する要素につき、図３における基本長方形が、式（１．１３）を適用することによって決定される。続いて、そのゆがんだ要素のアスペクト比Ａが、その基本長方形の縦辺と横辺との長さの比として決定される。斜傾性または台形性を伴うと、大きなアスペクト比Ａの影響が顕著になる。アスペクト比Ａは、図３を用いて説明すれば、ｈｘ／ｈｙとして計算される。
【０１７７】
その後、Ｓ１０４において、その計算されたアスペクト比Ａと、前記特定された変形モードとに基づき、今回の要素についての個別解析誤差ｅＡであってアスペクト比Ａにより表される要素ゆがみに起因するものが推定される。アスペクト比Ａと変形モードと個別解析誤差ｅＡとの関係であって、後述の相関推定プログラムの実行によって推定されたものがメモリ３２または記憶装置４０にストアされており、その関係に従い、今回のアスペクト比Ａと今回の変形モードとに対応する個別解析誤差ｅＡが決定されるのである。アスペクト比Ａと変形モードと個別解析誤差ｅＡとの関係は、図３７にグラフで表された関係に準じたものとして設定されている。
【０１７８】
続いて、Ｓ１０５において、今回の要素の反り度Ｗが計算される。
【０１７９】
反り度Ｗは、従前からの定義に従い、例えば、次のようにして計算される。
【０１８０】
すなわち、まず、今回の要素について投影平面が、今回の要素の４つの節点から等距離にあるように想定される。反り度Ｗは、単に、今回の要素の４つの節点のｚ座標値であり、そのｚ座標値は、平面ではない実際の要素と上記投影平面との距離に等しい。
【０１８１】
その後、Ｓ１０６において、Ｓ１０４に準じて、その計算された反り度Ｗと、前記特定された変形モードとに基づき、今回の要素についての個別解析誤差ｅＷであって反り度Ｗにより表される要素ゆがみに起因するものが、反り度Ｗと変形モードと個別解析誤差ｅＷとの関係であって後述の相関推定プログラムの実行によって推定されたものを用いて推定される。反り度Ｗと変形モードと個別解析誤差ｅＷとの関係も、図３７にグラフで表された関係に準じたものとして設定されている。
【０１８２】
続いて、Ｓ１０７において、今回の要素の斜傾度Ｓが計算される。斜傾度Ｓは、従前とは異なり、図３に示すように、今回の要素と同じアスペクト比Ａを有する長方形が基本長方形として想定され、今回の要素の一節点が基本長方形からｘ方向（横方向）にせん断変形した距離α１および、今回の要素の一節点が基本長方形からｙ方向にせん断変形した距離α２として計算される。
【０１８３】
その後、Ｓ１０８において、Ｓ１０４およびＳ１０６に準じて、その計算された斜傾度Ｓと、前記特定された変形モードとに基づき、今回の要素についての個別解析誤差ｅＳであって斜傾度Ｓにより表される要素ゆがみに起因するものが、斜傾度Ｓと変形モードと個別解析誤差ｅＳとの関係であって後述の相関推定プログラムの実行によって推定されたものを用いて推定される。斜傾度Ｓと変形モードと個別解析誤差ｅＳとの関係も、図３７にグラフで表された関係に準じたものとして設定されている。
【０１８４】
続いて、Ｓ１０９において、今回の要素の台形度Ｔが計算される。台形度Ｔは、従前とは異なり、図３に示すように、今回の要素と同じアスペクト比Ａを有する長方形が基本長方形として想定され、今回の要素の一節点が基本長方形からｘ方向（横方向）にテーパ変形した距離α３および、今回の要素の一節点が基本長方形からｙ方向にテーパ変形した距離α４として計算される。テーパ変形とは、横方向または縦方向において互いに対向する２つの要素辺間の相対的な広がりを意味する。
【０１８５】
その後、Ｓ１１０において、Ｓ１０４、Ｓ１０６およびＳ１０８に準じて、その計算された台形度Ｔと、前記特定された変形モードとに基づき、今回の要素についての個別解析誤差ｅＴであって台形度Ｔにより表される要素ゆがみに起因するものが、台形度Ｔと変形モードと個別解析誤差ｅＴとの関係であって後述の相関推定プログラムの実行によって推定されたものを用いて推定される。台形度Ｔと変形モードと個別解析誤差ｅＴとの関係は、正に、図３７にグラフで表された関係として設定されている。
【０１８６】
続いて、Ｓ１１１において、それら推定された４つの個別解析誤差ｅＡ，ｅＷ，ｅＳおよびｅＴに基づき、今回の要素の総合解析誤差Σｅが推定される。総合解析誤差Σｅは、それら４つの個別解析誤差ｅＡ，ｅＷ，ｅＳおよびｅＴの単純和として推定される。ただし、アスペクト比Ａ、反り度Ｗ、斜傾度Ｓおよび台形度Ｔ間の相乗効果を考慮し、重み付き和として総合解析誤差Σｅを推定することが可能である。
【０１８７】
その後、Ｓ１１２において、今回の要素につき、４つの個別解析誤差ｅＡ，ｅＷ，ｅＳおよびｅＴと総合解析誤差Σｅとがモニタ２２の画面上に数値により表示され、その結果、定量的に表示される。この表示は、すべての要素について一斉に行うことが可能である。
【０１８８】
なお付言すれば、解析誤差の定量表示のために、解析誤差を数値によって表示することは不可欠ではなく、例えば、グラフや図形の大きさや色を変化させることによって表示することが可能である。
【０１８９】
続いて、Ｓ１１３において、今回の要素について推定された総合解析誤差Σｅが許容値ＡＬを超えているか否かが判定される。超えている場合には、Ｓ１１４において、モニタ２２の画面上において今回の要素が他の要素から区別されるように表示される。これに対し、今回の要素について推定された総合解析誤差Σｅが許容値ＡＬを超えてはいない場合には、Ｓ１１４がスキップされる。
【０１９０】
いずれの場合にも、その後、Ｓ１１５において、今回の要素番号ｉが最大値ｉＭＡＸ以上であるか否かが判定される。Ｓ１０２ないしＳ１１４の実行が、今回の解析対象のための複数の要素すべてについて行われたか否かが判定されるのである。今回の要素番号ｉが最大値ｉＭＡＸ以上である場合には、以上で、この解析誤差推定プログラムの一回の実行が終了する。これに対し、今回の要素番号ｉが最大値以上である場合には、Ｓ１１６において、要素番号ｉが１だけインクリメントされた後、Ｓ１０２に戻る。
【０１９１】
図４２には、複数の要素について推定された解析誤差の一表示例が示されている。この表示例においては、複数の要素と複数種類の推定解析誤差との関係が表形式で表示され、さらに、総合解析誤差Σｅが許容値ＡＬを超える要素については、「適否」の欄に「Ｘ」が記入され、許容値ＡＬを超えない要素については、「Ｏ」が記入される。
【０１９２】
このように、各要素ごとに複数種類の推定解析誤差が表示されるため、ユーザは、総合解析誤差Σｅが許容値ＡＬを超えている要素につき、４つの個別解析誤差のうち最大のものを瞬時に発見でき、その要素を修正する際に参照すべき有効な要素ゆがみ情報を瞬時に発見できる。例えば、要素番号ｉが３である要素については、総合解析誤差Σｅが許容値ＡＬを超えていて、ゆがみが適当ではないことが記号によって表示されるとともに、それの４つの個別解析誤差のうちアスペクト比Ａに対応するものが最大であることが表示される。
【０１９３】
なお付言すれば、本実施形態においては、各要素のゆがみに起因する解析誤差が、そのゆがみのうち、アスペクト比Ａにより表現されるアスペクト比成分と、反り度Ｗにより表現される反り成分と、斜傾性Ｓにより表現される斜傾性成分と、台形性Ｔにより表現される台形性成分とのすべてが考慮されて推定される。
【０１９４】
ただし、それら４つの成分のうち、有限要素近似または解析のためのソフトウエア（以下、「解析ソフトウエア」という。）の種類の特徴に応じ、解析誤差に対して支配的であるもの、例えば、アスペクト比成分と台形性成分とは考慮するが、反り成分と斜傾性成分とについては考慮しないで解析誤差を推定する態様で本発明を実施することが可能である。前述の方法に従って推定された解析誤差は、要素剛性マトリクスを計算するとともに種々のゆがみ変数についてひずみエネルギ誤差を計算することにより、ユーザが関心を示す解析ソフトウエアを検証するために利用することが可能であり、その結果に基づき、ユーザは、支配的なゆがみモードを決定することが可能である。この態様によれば、本実施形態におけるより、解析誤差の計算を簡単にして高速化することが容易となる。
【０１９５】
図４３には、前記相関推定プログラムの内容がフローチャートで概念的に表されている。この相関推定プログラムは、要素のゆがみとそれに起因する解析誤差との間における関係を相関として推定するためにコンピュータ１２により実行される。
【０１９６】
ただし、この相関推定プログラムと前記解析誤差推定プログラムとを同じコンピュータによって実行することは不可欠ではない。相関は解析誤差推定プログラムの実行に不可欠であるが、解析誤差推定プログラムの実行に先立って存在すれば足り、しかも、その解析誤差推定プログラムの実行中に相関を更新することは不可欠ではない。
【０１９７】
したがって、コンピュータ１２とは別のコンピュータによる相関推定プログラムの実行によって相関を取得し、その取得された相関を表すデータをコンピュータ１２のメモリ３２または記憶装置４０に予めストアすることによっても、解析誤差を推定することが可能である。
【０１９８】
この相関推定プログラムは、この数値解析方法を実施するごとにコンピュータ１２によって実行することは不要であり、将来における数値解析方法の複数回の実施（解析対象が同じであるか否かを問わない）に先立ち、一回限り行えば足りる。なぜなら、この相関推定プログラムの実行によって推定される相関は、同じ解析対象についての複数回の数値解析に共通して使用できることはもとより、異なる解析対象についての複数回の数値解析にも共通して使用できるからである。
【０１９９】
この相関推定プログラムにおいては、概略的に説明すれば、ある四辺形要素が注目されるとともに、その要素につき、各変形モードごとに、アスペクト比Ａの設定値と反り度Ｗの設定値と斜傾度Ｓの設定値と台形度Ｔの設定値との組合せをステップ幅ずつ更新しながら解析誤差が推定される。
【０２００】
この相関推定プログラムにおいては、さらに、解析誤差の推定のため、各変形モードにつき、要素のひずみエネルギの解析解と、要素のひずみエネルギの近似解とが計算される。要素のひずみエネルギの解析解は、有限要素解析によることなく、その要素の応力およびひずみの理論解を用いて計算される。要素のひずみエネルギの近似解は、有限要素近似により、要素の節点変位の理論解と、特定の有限要素ソフトウエアのシェル要素の要素剛性マトリクスとを用いて計算される。例えば、ＭＳＣ．Ｎａｓｔｒａｎの要素剛性マトリクス、またはＡＢＡＱＵＳの要素剛性マトリクス、要素の種類で言えば変位法による要素剛性マトリクス、応力仮定法による要素剛性マトリクスなどの間において要素剛性マトリクスＫｅを交換することにより、解析ソフトウエアの種類ごとに解析ソフトウエアのゆがみによる要素の特性を調べることができる。その結果から、いずれかの解析ソフトウエアにおいては、いずれのゆがみの影響が他のゆがみより強調されるかなどを調べることが可能である。
【０２０１】
ところで、解析対象の形状を近似的に表現するメッシュモデルには多くのゆがんだ要素が含まれているが、本実施形態においては、１つの有限要素のゆがみ特性を評価するために、１つの要素に関しての要素剛性マトリクスＫｅを用いることにより、その要素内に生じるひずみエネルギの近似解が計算される。
【０２０２】
通常の有限要素解析においては、１つの解析対象を互いに共同して構成する複数の要素に関する複数の要素剛性マトリクスＫｅを合成して全体剛性マトリクスＫを作成し、それを境界条件や荷重条件と共に用いて各要素が解析されるが、本実施形態においては、上述のように、１つの要素が代表的に用いられ、それに関する要素剛性マトリクスＫｅを用いてひずみエネルギの近似解が計算される。
【０２０３】
さらに、この相関推定プログラムにおいては、それら解析解と近似解との差を用いて解析誤差が計算される。具体的には、それら解析解と近似解との差を解析解で割り算した値がひずみエネルギ誤差として計算され、それが解析誤差とされる。
【０２０４】
この相関推定プログラムにおいては、さらに具体的には、まず、Ｓ２０１において、今回注目している要素の変形モードの種類を表す変形モード値ｊが１にセットされる。変形モード値ｊは、例えば、それが１のときに引張りを意味し、２のときに曲げを意味し、３のときにせん断を意味し、４のときにねじりを意味するように設計できる。
【０２０５】
次に、Ｓ２０２ないしＳ２０５においてそれぞれ、１つの要素、すなわち、解析誤差の推定結果が将来利用されるであろう解析対象の形状を互いに共同して近似的に表現する複数の要素を代表する代表要素のアスペクト比Ａの今回値と、反り度Ｗの今回値と、斜傾度Ｓの今回値と、台形度Ｔの今回値とがそれぞれ設定される。
【０２０６】
続いて、Ｓ２０６において、代表要素につき、例えば前記式（３．４）により表されるひずみエネルギ式を用いることにより、前述のようにして、ひずみエネルギの解析解Ｕｅ（ｕ）が計算される。
【０２０７】
その後、Ｓ２０７において、代表要素につき、例えば前記式（３．１０）により表されるひずみエネルギ式を用いることにより、前述のようにして、ひずみエネルギの近似解Ｕｅ（ｕｈ）が計算される。
【０２０８】
続いて、Ｓ２０８において、代表要素につき、前記式（３．１）に従い、ひずみエネルギの解析解Ｕｅ（ｕ）の計算値からひずみエネルギの近似解Ｕｅ（ｕｈ）の計算値を差し引いた値をその解析解Ｕｅ（ｕ）の計算値で割り算することにより、ひずみエネルギ誤差Ｅｅが計算される。
【０２０９】
その後、Ｓ２０９において、その計算されたひずみエネルギ誤差Ｅｅが、解析誤差として、今回の変形モードとアスペクト比Ａの今回値と反り度Ｗの今回値と斜傾度Ｓの今回値と台形度Ｔの今回値とに関連付けてメモリ３２または記憶装置４０にストアされる。
【０２１０】
すなわち、それら変形モードｉとアスペクト比Ａの今回値と反り度Ｗの今回値と斜傾度Ｓの今回値と台形度Ｔの今回値と解析誤差の今回値との関係が、今回の変形モードについての要素ひずみ−解析誤差間関係の一例を構成しているのである。
【０２１１】
続いて、Ｓ２１０において、アスペクト比Ａと反り度Ｗと斜傾度Ｓと台形度Ｔとのすべてについて各設定値の更新が終了したか否かが判定される。すなわち、すべての幾何学変数について予定範囲での変化が与えられたか否かが判定されるのである。
【０２１２】
今回は、各設定値の更新が終了していないと仮定すれば、Ｓ２１１において、アスペクト比Ａと反り度Ｗと斜傾度Ｓと台形度Ｔとのうち該当するものの設定値が、予め定められたステップ幅で更新される。
【０２１３】
その後、Ｓ２０２ないしＳ２１１の実行が繰り返される結果、アスペクト比Ａと反り度Ｗと斜傾度Ｓと台形度Ｔとのすべてについて各設定値の更新が終了すれば、Ｓ２１２において、変形モード値ｊの今回値が最大値ｊＭＡＸであるか否かが判定される。すなわち、想定可能なすべての変形モードについて解析誤差の推定が終了したか否かが判定されるのである。
【０２１４】
今回は、変形モード値ｊの今回値が最大値ｊＭＡＸではないと仮定すれば、Ｓ２１３において、変形モード値ｊが１だけインクリメントされ、その後、Ｓ２０２に戻る。
【０２１５】
その後、Ｓ２０２ないしＳ２１３の実行が繰り返される結果、想定可能なすべての変形モードについて解析誤差の推定が終了すれば、以上で、この相関推定プログラムの一回の実行が終了する。
【０２１６】
以上の説明から明らかなように、本実施形態においては、その相関推定プログラムが前記（１）項または（９）項における「推定工程」の一例と、前記（３）項または（４）項における「計算工程」の一例とをそれぞれ構成しているのである。
【０２１７】
さらに、本実施形態においては、図４１におけるＳ１０２が前記（１１）項における「変形モード特定工程」の一例を構成し、Ｓ１１２ないしＳ１１４が同項における「表示工程」の一例を構成し、Ｓ１０３ないしＳ１１１が前記（１３）項における「推定工程」の一例を構成し、Ｓ１０３ないしＳ１１０がそれぞれ前記（１５）項における「個別解析誤差推定工程」の一例を構成し、Ｓ１１１が前記（１６）項における「総合解析誤差推定工程」の一例を構成し、Ｓ１１４が前記（１７）項における「区別表示工程」の一例を構成しているのである。
【０２１８】
さらに、本実施形態においては、図４０におけるＳ１が前記（１８）項における「要素分割工程」の一例を構成し、図４１におけるＳ１０２が同項における「変形モード特定工程」の一例を構成し、同図におけるＳ１１２ないしＳ１１４が同項における「表示工程」の一例を構成し、図４１におけるＳ５が同項における「解析工程」の一例を構成しているのである。
【０２１９】
さらに、本実施形態においては、前述の４つのモードが、前記（２）項における「複数の基本変形モード」の一例であり、前述のゆがみ変数（幾何学変数）α１ないしα４が、同項における「複数の基本ゆがみ変数」の一例であり、前述のひずみεｘおよびεｙが、同項における「垂直ひずみ」の一例であり、前述のひずみγｘｙが、同項における「せん断ひずみ」の一例である。
【０２２０】
さらに、本実施形態においては、図１９における式（１．３８）ないし式（１．４０）が、前記（２）項における「面状要素の垂直ひずみおよびせん断ひずみの各成分を定式化」した式の一例であり、同図における式（３．４）が、同項における「面状要素のひずみエネルギ」を計算する式の一例である。
【０２２１】
以上、本発明の具体的な一実施形態を図面に基づいて詳細に説明したが、これは例示であり、前記［課題を解決するための手段および発明の効果］の欄に記載の態様を始めとして、当業者の知識に基づいて種々の変形、改良を施した他の形態で本発明を実施することが可能である。
【図面の簡単な説明】
【０２２２】
【図１】本発明の一実施形態に従う数値解析方法において要素ゆがみに起因した解析誤差を推定する推定原理の理論的根拠を説明するための式（１．１）ないし（１．３）を示す図である。
【図２】上記推定原理の理論的根拠を説明するための式（１．４）ないし（１．６）を示す図である。
【図３】上記推定原理における四辺形要素のゆがみの種類を説明するための図である。
【図４】上記推定原理の理論的根拠を説明するための式（１．７）ないし（１．１０）を示す図である。
【図５】上記推定原理の理論的根拠を説明するための式（１．１１）ないし（１．１４）を示す図である。
【図６】上記推定原理における局所座標系ｘ−ｙ−ｚの設定とその局所座標系上における四辺形要素の反りを説明するための図である。
【図７】上記推定原理の理論的根拠を説明するための式（１．１５）ないし（１．１８）を示す図である。
【図８】上記推定原理の理論的根拠を説明するための式（１．１９）を示す図である。
【図９】上記推定原理の理論的根拠を説明するための任意の四角形要素の節点の全体座標、局所座標系の原点および軸ベクトルの例を表形式で示す図である。
【図１０】上記推定原理の理論的根拠を説明するための任意の四角形要素の節点の局所座標の例を表形式で示す図である。
【図１１】上記推定原理の理論的根拠を説明するための任意の四角形要素の局所座標系を示す図である。
【図１２】上記推定原理において任意の四辺形要素が複数種類のゆがみに分解される様子を説明するための図である。
【図１３】上記推定原理の理論的根拠を説明するための式（１．２０）および（１．２１）を示す図である。
【図１４】上記推定原理における複数のゆがみモードをヤコビのマトリクスと共に示す図である。
【図１５】上記推定原理の理論的根拠を説明するための式（１．２２）ないし（１．２５）を示す図である。
【図１６】上記推定原理の理論的根拠を説明するための式（１．２６）ないし（１．２９）を示す図である。
【図１７】上記推定原理の理論的根拠を説明するための式（１．３０）ないし（１．３２）を示す図である。
【図１８】上記推定原理の理論的根拠を説明するための双一次形式のアイソパラメトリックな四角形要素の複数種類の変形モードを説明するための図である。
【図１９】上記推定原理の理論的根拠を説明するための式（１．３３）ないし（１．３５）を示す図である。
【図２０】上記推定原理の理論的根拠を説明するための式（１．３６）および（１．３７）を示す図である。
【図２１】上記推定原理の理論的根拠を説明するための式（１．３８）ないし（１．４０）を示す図である。
【図２２】上記推定原理の理論的根拠を説明するための式（１．４１）を示す図である。
【図２３】上記推定原理の理論的根拠を説明するための式（２．１）ないし（２．４）を示す図である。
【図２４】上記推定原理の理論的根拠を説明するための板要素の５自由度を説明するための図である。
【図２５】上記推定原理の理論的根拠を説明するための式（２．５）および（２．６）を示す図である。
【図２６】上記推定原理の理論的根拠を説明するための式（２．７）ないし（２．９）を示す図である。
【図２７】上記推定原理の理論的根拠を説明するための式（２．１０）ないし（２．１３）を示す図である。
【図２８】上記推定原理の理論的根拠を説明するための四辺形要素が単純曲げとねじりと点負荷とによってゆがむ様子を説明するための斜視図である。
【図２９】上記推定原理の理論的根拠を説明するための式（３．１）ないし（３．４）を示す図である。
【図３０】上記推定原理の理論的根拠を説明するための式（３．５）ないし（３．７）を示す図である。
【図３１】上記推定原理の理論的根拠を説明するための式（３．８）を示す図である。
【図３２】上記推定原理の理論的根拠を説明するための複数の解析的な変位関数の例を表形式で示す図である。
【図３３】上記推定原理の理論的根拠を説明するための式（３．９）ないし（３．１０）を示す図である。
【図３４】上記推定原理に従って推定されたひずみエネルギ誤差が要素の台形性ゆがみに応じて変化する様子をせん断修正なしの状態で示すグラフである。
【図３５】上記推定原理の理論的根拠を説明するための四辺形要素が引張り、単純曲げおよびせん断によってゆがむ様子を説明するための図である。
【図３６】上記推定原理の理論的根拠を説明するための式（３．１１）を示す図である。
【図３７】上記推定原理に従って推定されたひずみエネルギ誤差が要素の台形性ゆがみに応じて変化する様子をせん断修正ありの状態で示すグラフである。
【図３８】前記数値解析方法を実施するために使用されるコンピュータ・システム１０を示す斜視図である。
【図３９】図３８のコンピュータ・システム１０のハードウエア構成を概念的に示すブロック図である。
【図４０】前記数値解析方法の内容を概念的に表すフローチャートである。
【図４１】前記数値解析方法を実施するために図３９のコンピュータ１２により実行される解析誤差推定プログラムの内容を概念的に表すフローチャートである。
【図４２】図４１におけるＳ１１２およびＳ１１４の実行により一表示例を示す図である。
【図４３】前記数値解析方法を実施するために図３９のコンピュータ１２により実行される相関推定プログラムの内容を概念的に表すフローチャートである。
【図４４】上記推定原理の理論的根拠を説明するための複数の解析的な変位関数の例を表形式で示す図である。
【図４５】上記推定原理に従って推定されたひずみエネルギ誤差が要素の台形性ゆがみに応じて変化する様子を示すグラフである。
【図１】
2004303227.tif 000003
【図２】
2004303227.tif 000004
【図３】
2004303227.tif 000005
【図４】
2004303227.tif 000006
【図５】
2004303227.tif 000007
【図６】
2004303227.tif 000008
【図７】
2004303227.tif 000009
【図８】
2004303227.tif 000010
【図９】
2004303227.tif 000011
【図１０】
2004303227.tif 000012
【図１１】
2004303227.tif 000013
【図１２】
2004303227.tif 000014
【図１３】
2004303227.tif 000015
【図１４】
2004303227.tif 000016
【図１５】
2004303227.tif 000017
【図１６】
2004303227.tif 000018
【図１７】
2004303227.tif 000019
【図１８】
2004303227.tif 000020
【図１９】
2004303227.tif 000021
【図２０】
2004303227.tif 000022
【図２１】
2004303227.tif 000023
【図２２】
2004303227.tif 000024
【図２３】
2004303227.tif 000025
【図２４】
2004303227.tif 000026
【図２５】
2004303227.tif 000027
【図２６】
2004303227.tif 000028
【図２７】
2004303227.tif 000029
【図２８】
2004303227.tif 000030
【図２９】
2004303227.tif 000031
【図３０】
2004303227.tif 000032
【図３１】
2004303227.tif 000033
【図３２】
2004303227.tif 000034
【図３３】
2004303227.tif 000035
【図３４】
2004303227.tif 000036
【図３５】
2004303227.tif 000037
【図３６】
2004303227.tif 000038
【図３７】
2004303227.tif 000039
【図３８】
2004303227.tif 000040
【図３９】
2004303227.tif 000041
【図４０】
2004303227.tif 000042
【図４１】
2004303227.tif 000043
【図４２】
2004303227.tif 000044
【図４３】
2004303227.tif 000045
【図４４】
2004303227.tif 000046
【図４５】
2004303227.tif 000047