有限要素法を用いた解析方法、及び有限要素法を用いた解析演算プログラム
【課題】メッシュ・パターンによる解析誤差を低減し、品質が低い要素（品質の低いメッシュ・パターン）からも高精度な数値解が得られる有限要素法を用いた解析方法を提供する。
【解決手段】有限要素法を用いた解析方法にあっては、解析対象の解析領域を選定する選定工程（Ｓ１）、解析領域を計算対象として複数の要素に分割する分割工程（Ｓ２）、各要素のマトリクスを作成する要素マトリクス作成工程（Ｓ３）、Galerkin重み関数と一般関数との積からなる一般関数項を積分する一般関数項積分工程（Ｓ４）、各要素のマトリクスの和と一般関数項を積分した値の和とに基づき連立方程式を作成する連立方程式作成工程（Ｓ５）、連立方程式から数値解を得る演算工程（Ｓ７，Ｓ８）を有する。該一般関数項積分工程で、Galerkin有限要素法による二階微分項の離散化結果を基に定義された節点領域の考え方を導入し、要素代表値を用いた一般関数項を積分する。
【解決手段】有限要素法を用いた解析方法にあっては、解析対象の解析領域を選定する選定工程（Ｓ１）、解析領域を計算対象として複数の要素に分割する分割工程（Ｓ２）、各要素のマトリクスを作成する要素マトリクス作成工程（Ｓ３）、Galerkin重み関数と一般関数との積からなる一般関数項を積分する一般関数項積分工程（Ｓ４）、各要素のマトリクスの和と一般関数項を積分した値の和とに基づき連立方程式を作成する連立方程式作成工程（Ｓ５）、連立方程式から数値解を得る演算工程（Ｓ７，Ｓ８）を有する。該一般関数項積分工程で、Galerkin有限要素法による二階微分項の離散化結果を基に定義された節点領域の考え方を導入し、要素代表値を用いた一般関数項を積分する。
【発明の詳細な説明】
【技術分野】
【０００１】
本発明は、有限要素法を用いた解析方法、及び有限要素法を用いた解析演算プログラムに係り、特に要素の形状（メッシュ・パターン）によって生じる誤差の低減を図った有限要素法を用いた解析方法、及び有限要素法を用いた解析演算プログラムに関する。
【背景技術】
【０００２】
一般的に有限要素法では、形状関数と同形なGalerkin重み関数を微分方程式の各項にかけ、節点周りの要素領域にわたって積分することにより数値解を求めている。有限要素法は、誤差評価が明瞭で、汎用性に富むなどの特徴があって、最も広く使用される数値解法の一つになっている。しかし、メッシュ・パターンによって解析結果が異なったり、数値誤差が大きくなる問題が依然として解決されていない。
【０００３】
［従来の有限要素法による二次元三角形の問題］
例えば、線形三角形メッシュ上でポアソン（Poisson）の方程式を離散化した場合には、もともと同じ解が出るべきところから、要素形状によって異なる解析結果が出されてしまう（非特許文献１参照）。また、非特許文献２に示す２パターンのメッシュから得た代数方程式を比較すると、二階微分項の離散化結果は同形になっているものの、ソース項の離散化結果が異なっている。ソース項に重み関数をかけて積分した結果、ソース量が要素形状と関係なくほぼ均等的に各節点の代数方程式に分配され、ほかの項との整合性が取れず、節点レベルでは保存法則が満たされない状況になっている。この不整合性の詳細は、以下のとおりである。
【０００４】
二次元ポアソンの方程式は下記のものである。
【数１】
【０００５】
この方程式は、多分野の物理現象を記述するものであるが、ここでは、熱伝導問題を例に説明を展開し、Ｔは温度で、λは熱伝導率で、Ｑはソース項と呼ばれる単位時間、単位面積あたりの発熱量である。式（１）に有限要素法を適用すると、考えている対象節点Ｐについての積分方程式として次式が得られる。
【数２】
【０００６】
ここで、Ωは図１に示す節点Ｐ周りの要素からなす領域であり、ＷはΩ上で定義された節点Ｐで１を、Ωの境界で０をとるGalerkin重み関数である。Ωにわたる積分は、各要素上の積分の和になるので、式（２）を式（３）に書き換えることができる。
【数３】
【０００７】
ここで、ｅは節点Ｐ周りの要素を表わし、
【数４】
は各要素上の演算の和を意味する。要素ｅを図２の線形三角形要素に当てはめた場合には、式（３）左辺の要素積分の結果を次式のように整理することができる。
【数５】
ここでは、局所座標系を使用しており、対象節点Ｐを局所節点１とし、数字の添え字で節点番号を表わす。
【０００８】
式（４）右辺の１項目中の
【数６】
は図３に示す中点Ｉでの節点１から節点２に向かう熱流束の２次精度の近似であり、その係数ｄＣＩは、点Ｉから外心Ｃまでの距離と等しくなっている。すなわち、式（４）右辺の１項目は点Ｉでの２次精度の熱流束を用いて、垂直な線分ＩＣを通過する熱量を評価するものである。同様に、２項目は、点Ｋでの２次精度の熱流束を用いて、垂直な線分ＣＫを通過する熱量を評価するものである。以上の考察をまとめると、式（４）の右辺は、ＩＣとＣＫという検査線を通して流出した熱量であることが分かる。この際、四辺形１ＩＣＫを節点１の節点領域と呼び、ＮＤ１で表わす。正三角形要素では外心が幾何中心と重なるため、節点領域は要素面積の３分の１になる。直角三角形要素では外心が辺中点に重なり、鈍角三角形要素では外心が要素の外に出る。これらの要素における節点１の節点領域は図４と図５に示すようになる。
【０００９】
一方、式（３）右辺のソース項の要素積分により、要素内の発熱量が要素形状と関係なくほぼ均等的に３節点に分配される。特に、Ｑが要素内で定数の場合には、ソース項の要素積分は次式になる。
【数７】
ここで、Ｓは要素面積である。すなわち、いかなる要素形状においても、節点に寄与された熱量は要素の３分の１の面積から発生したものだけである。
【００１０】
ここで、正三角形以外の要素では、検査線で決めた節点領域の面積とソース項の面積との間は不釣り合いが発生する。同じことは、節点Ｐ周りのほかの要素からも生じるので、保存法則が満たされず、ソース項と他の項の間に不整合性が表れた。その結果、新たな数値誤差が引き起こされる。図６には、有限要素法から得られた温度分布の一例を示す。０≦ｘ≦１、０≦ｙ≦１の正方形領域において、式（６）のソース分布と式（７）の境界条件を用いて、λ＝１とした解析結果を実線で示し、点線で要素形状を示す。例えば、○印を付けた中央部の４点では同値の厳密解は得られるが、局所メッシュ・パターンの影響で数値解には差が出ている。
【数８】
【数９】
【００１１】
［従来の二次元三角形問題の改善方法］
前述の不整合性問題を解決するため、以下の２つの改善策が考えられる。
（改善方法１）：式（３）右辺の要素積分
【数１０】
を
【数１１】
に置き換える。
（改善方法２）：上記改善方法１をさらに修正したものであり、鈍角三角形要素における要素外での積分を避けるため、ソース項の積分を図７に示す要素内の部分で行う。それに合わせて式（４）右辺のｄＣＩとｄＣＫを中点から外心までの長さから図７に示す長さに修正する。
【００１２】
上記改善方法１では、全ての線形三角形要素において保存法則が２次精度で満たされ、各項間の整合性が維持できるが、鈍角三角形要素の場合には、ソース項の積分範囲が要素外までに広がることが、数値振動などを引き起こす原因になっている。また、上記改善方法２では、左辺のマトリックスが強制的に変えられたため、楕円型演算子の最適な離散化方法として評価されてきた有限要素法の優位性が損なわれる。結局のところ、上記改善方法１，２では満足できる数値解が得られず、実用上では鈍角三角形要素の使用を避けており、言い換えると、鈍角三角形要素を用いないメッシュ・パターンを時間をかけて作成する必要がある。
【００１３】
［従来の有限要素法による三次元四面体の問題］
三次元のポアソンの方程式は次式になる。
【数１２】
【００１４】
有限要素法を適用すると対象節点についての積分式として次式が得られる。
【数１３】
【００１５】
式（９）は式（１０）に示す節点周りの要素積分の和で表現することができる。
【数１４】
【００１６】
式（１０）を、例えば、図８のような三つの面が節点１で垂直に交わる線形四面体要素に適用すると、左辺の要素積分は次式になる。
【数１５】
【００１７】
式（１１）の右辺から、節点１から節点２、３、４に向かう熱流束
【数１６】
の通過面積はそれぞれ
【数１７】
になることが分かる。同様な方法で、ほかの節点間に伝わる熱流束の通過面積も求められる。それらを表１にまとめて示す。
【００１８】
【表１】
【００１９】
一方、式（１０）右辺のソース項の要素積分を行うと、要素内の発熱量が要素形状と関係なく、ほぼ均等的に４節点に分配される。特に、Ｑが要素内で定数の場合には、ソース項の要素積分は次式になる。
【数１８】
【００２０】
ここでのＶは要素体積である。すなわち、いかなる要素形状においても、節点に寄与された熱量は要素の４分の１の体積から発生したものだけである。ここで、保存法則の観点から表１と式（１２）を考察すると、節点によって大きく変る熱流束の通過面積と均等的に分配されるソース項の間に不整合性が表れたことが分かる。その結果、新たな数値誤差が引き起こされる。図９には、０≦ｘ≦１、０≦ｙ≦１、０≦ｚ≦１の立方体領域において、λ＝１とし、式（１３）のソース分布と式（１４）の境界条件に有限要素法を適用して得られた温度分布を示す。例えば、図９（ｂ）の○印を付けた４点では同値の厳密解は得られるが、局所メッシュ・パターンの影響で数値解は異なっている。
【数１９】
【数２０】
【００２１】
［従来の三次元四面体問題の改善方法］
この不整合性の三次元四面体の問題の改善方法としては、図１０のように、対象節点１、四面体の外心Ｃ、要素表面三角形の外心Ｃ２、Ｃ３、Ｃ４と要素辺の中点Ｉ、Ｊ、Ｋの合計８点から成形した六面体を節点領域として定義し、それに合わせて、マトリックス成分を節点領域から算出した流束の通過面積を節点間の距離で割ったものに修正する方法が考えられる。マトリックス成分を変える理由としては、三次元有限要素法の離散化結果には二次元と同様な解釈が適用できないからである。すなわち、式（１１）左辺を積分した結果は、正四面体要素の場合を除いて、辺中点での熱流束と、外心および辺中点からなす四辺形の面積との掛け算にならない。例えば、図８の要素の場合には、節点１と節点２間の熱流束の通過面積は
【数２１】
になっているが、外心Ｃ、Ｃ３、Ｃ４と中点Ｉとによりできた四辺形の面積は
【数２２】
になる。
【００２２】
この三次元問題の改善方法は、保存法則を２次精度で満足するが、デロネー（Delaunay）分割の要素形状にしか適用できないという欠点が挙げられる。また、たとえば、図８のような四面体要素に適用すると、従来の有限要素法では対角線以外の成分がゼロになることで済むが、この改善方法では、マトリックス成分を変えたことにより、対角線以外のところで数値解析に悪い影響を及ぼす正成分が発生する。
【００２３】
［従来の有限要素法による二次元四角形の問題］
上述した二次元の線形三角形メッシュ及び三次元の線形四面体メッシュについては、上述したような改善方法も提案されているが（非特許文献３、非特許文献４参照）、二次元の線形四辺形メッシュについては議論されていないようである。有限要素法によりPoissonの方程式を離散化した場合には、全ての項に同じGalerkin重み関数をかけて積分するため、ソース項が要素形状に関係なく同要素が所有する各節点の代数方程式にほぼ均等的に分配されることがあり、ほかの項の離散化との整合性が取れず、節点レベルでは保存則が満たされない状況になっている。その詳細は、以下のとおりである。
【００２４】
上述したように、二次元Poissonの方程式は下記のものである。
【数２３】
【００２５】
この方程式は、多分野の物理現象を記述するものであるが、ここでは、熱伝導問題を例に説明を展開し、Ｔは温度で、λは熱伝導率で、Ｑはソース項と呼ばれる単位時間、単位面積あたりの発熱量である。式(1)に有限要素法を適用すると、考えている対象節点Ｐについての積分方程式は、上述した式（３）と同じ式（３）’が得られる。
【数２４】
【００２６】
ここで、ｅは図１１に示す節点Ｐ周りの要素を表わし、Σは各要素上の演算の和を意味する。また、ＷＰは節点Ｐ周りの要素領域上で定義されたGalerkin重み関数であり、節点Ｐでは１を、図１１の局所領域の境界ではゼロをとる線形関数である。部分積分とガウス・グリーン定理を利用して式（３）’左辺の積分を行い、境界積分を別途考えると、次式が得られる。
【数２５】
【００２７】
要素ｅを図１２の双一次長方形要素に当てはめた場合には、式（１５）左辺の要素積分は次式となる。ただし、対象節点Ｐの要素における局所番号を１とする。
【数２６】
【００２８】
式（１６）右辺１項目の物理的な意味は、中点Ｉでの節点１から節点２に向かう２次精度の熱流束と、中点Ｋでの節点４から節点３に向かう２次精度の熱流束との重み付き平均により、垂直線分ＩＯを通過する熱量である。図１２中の点Ｏは、四節点の座標の平均値に位置する。また、２項目の物理的な意味は、中点Ｌでの節点１から節点４に向かう２次精度の熱流束と、中点Ｊでの節点２から節点３に向かう２次精度の熱流束との重み付き平均により、垂直線分ＯＬを通過する熱量である。ここで、検査線ＩＯＬと要素辺Ｌ１Ｉで囲う領域を節点１のコントロール・ボリュームと称しＣＶ１と記す。要素形状の対称性を考えると、ほかの節点を対象節点とした検査線の検討結果が、図１２の細い実線で示すようになり、ＣＶ１＝ＣＶ２＝ＣＶ３＝ＣＶ４＝Ｓ／４となることが分かる。
【００２９】
ここで、Ｓは要素面積であり、ＣＶの添え字で所属節点を表わす。一方、Ｑが要素内で定数の場合には、式（１５）右辺の要素積分は次式になる。
【数２７】
すなわち、節点１に寄与された熱量はＳ／４から発生したものである。ほかの節点への寄与量も同じであることが重み関数の定義から分かる。各節点のコントロール・ボリュームの面積とソース項の面積とが一致するため、保存則は２次精度で満たされている。これらの結果を表２にまとめて示す。
【００３０】
【表２】
【００３１】
なお、図１３に示すような双一次平行四辺形要素の場合には、有限要素法による式（１５）左辺の要素積分は次式になる。
【数２８】
【００３２】
式（１８）右辺の１項目は、中点Ｉでの節点１から節点２に向かう２次精度の熱流束により、垂直線分ＩＪを通過する熱量である。２項目は点Ｏでの節点１から節点３に向かう２次精度の熱流束により垂直線分ＪＬを通過する熱量である。結局、検査線はＩＪＬとなり、ＣＶ１は３Ｓ／８になる。同様に、節点２を対象節点として検討すると、式（１５）左辺の要素積分は次式になる。
【数２９】
検査線はＩＪとなり、ＣＶ２はＳ／８になる。残る二節点については、要素形状の対称性を考えると、図１３の細い実線で示す検査線を得ることができる。
【００３３】
一方、Ｑが要素内で定数の場合には、ソース項の要素積分結果はＳ／４から発生した熱量になる。上記表２にこれらの結果をまとめて示す。節点に寄与された熱量の発生面積とコントロール・ボリュームの面積とが一致しないため、保存則は満たされていない。
【００３４】
図１４に示す節点３が節点４に限りなく接近する四辺形要素の場合には、有限要素法による式（１５）左辺の要素積分は次式になる。ただし、Ｔ３＝Ｔ４と扱うことができる。
【数３０】
【００３５】
式（２０）右辺の１項目は、中点Ｉでの節点１から節点２に向かう２次精度の熱流束により、垂直線分ＩＪを通過する熱量である。２項目は点Ｌでの節点１から節点４に向かう２次精度の熱流束により垂直線分ＪＬを通過する熱量である。結局、検査線はＩＪＬとなり、ＣＶ１はＳ／２になる。同様に、節点２を対象節点として検討すると、式（１５）左辺の要素積分は次式になる。
【数３１】
【００３６】
その検査線はＩＪとなり、ＣＶ２はＳ／４になる。また、Ｔ３＝Ｔ４を考慮して節点３と４を対象節点として積分した結果をそれぞれ式（２２）と式（２３）に示す。
【数３２】
【数３３】
【００３７】
両節点が重なり、Ｔ３＝Ｔ４であるため、検査線について検討するとき、代数方程式を足し合わせて扱うことができる。式（２２）と式（２３）との右辺の合計から検査線としてＪＬが得られ、ＣＶ３＋ＣＶ４＝Ｓ／４になることが分かる。図１４には細い実線で検査線を示す。
【００３８】
一方、Ｑが要素内で定数の場合には、節点１と節点２に寄与されるソース量の面積がそれぞれＳ／３になる。節点３と節点４に寄与されるソース量の面積の和がＳ／３になる。これらの結果を上記表１にまとめて示すが、保存則は満たされていないことが分かる。
【００３９】
なお、このような２次精度で保存則が満たされない問題は、ほかの一般形状の四辺形要素にも存在すると考えられるが、理論上では、コントロール・ボリュームについて調べることが困難である。
【００４０】
［従来の有限要素法による三次元六面体及び三次元五面体の問題］
上述した二次元の線形三角形メッシュ及び三次元の線形四面体メッシュについては、上述したような改善方法も提案されているが（非特許文献３、非特許文献４参照）、三次元解析によく使用される六面体要素と五面体要素については議論されていないようである。本発明では、これらの要素を検討対象とする。
【００４１】
上述したように三次元Poissonの方程式は上記式（８）のものである。
【数３４】
【００４２】
この方程式は、多分野の物理現象を記述するものであるが、ここでは、熱伝導問題を例に説明を展開し、Tは温度で、λは熱伝導率で、Ｑはソース項と呼ばれる単位時間、単位面積あたりの発熱量である。式（８）に有限要素法を適用すると、考えている対象節点Ｐについての積分方程式として次式が得られる。
【数３５】
【００４３】
ここで、ｅは図１５に示す節点Ｐ周りの要素を表わし、Σは各要素上の演算の和を意味する。また、ＷＰは節点Ｐ周りの要素領域上で定義されたGalerkin重み関数で、節点Ｐでは１を、周りの要素から構成した領域の境界ではゼロをとる線形関数である。部分積分とガウス・グリーン定理を利用して式（１０）’左辺の積分を行い、境界積分を別途考えると、次式が得られる。
【数３６】
【００４４】
［六面体要素の場合］
要素ｅを図１６（ａ）の線形長方体要素に当てはめた場合には、式（２４）左辺の要素積分結果を次式に整理することができる。ただし、対象節点Ｐの要素における局所番号を１とする。
【数３７】
【００４５】
図１６(ｂ)に示すように、点Ｉｉｊで節点ｉとｊを結ぶ線分の中点を表わしておく。式(２５)右辺１項目の物理的な意味は、中点Ｉ１２での節点１から節点２に向かう２次精度の熱流束λ（Ｔ１−Ｔ２）／２ｈｘ、中点Ｉ４３での節点４から節点３に向かう２次精度の熱流束λ（Ｔ４−Ｔ３）／２ｈｘ、中点Ｉ５６での節点５から節点６に向かう２次精度の熱流束λ（Ｔ５−Ｔ６）／２ｈｘ、と中点Ｉ８７での節点８から節点７に向かう２次精度の熱流束λ（Ｔ８−Ｔ７）／２ｈｘ、との重み付き平均により、熱流束に垂直で中点Ｉ１２を通る断面積ｈｙｈｚを通過する熱量である。また、２項目の物理的な意味は、中点Ｉ１４での節点１から節点４に向かう２次精度の熱流束、中点Ｉ２３での節点２から節点３に向かう２次精度の熱流束、中点Ｉ５８での節点５から節点８に向かう２次精度の熱流束、と中点Ｉ６７での節点６から節点７に向かう２次精度の熱流束、との重み付き平均により、熱流束に垂直で中点Ｉ１４を通る断面積ｈｘｈｚを通過する熱量である。また、３項目の物理的な意味は、中点Ｉ１５での節点１から節点５に向かう２次精度の熱流束、中点Ｉ２６での節点２から節点６に向かう２次精度の熱流束、中点Ｉ４８での節点４から節点８に向かう２次精度の熱流束、と中点Ｉ３７での節点３から節点７に向かう２次精度の熱流束、との重み付き平均により、熱流束に垂直で中点Ｉ１５を通る断面積ｈｘｈｙを通過する熱量である。ここで、断面積ｈｙｈｚ、ｈｘｈｚとｈｘｈｙを検査面として考え、これらの検査面が囲う領域を節点１のコントロール・ボリュームと称しＣＶ１と記す。このコントロール・ボリュームは、２次精度の熱流束により定義されたものである。以上の議論から、ＣＶ１は図１６(ｂ)に示す辺長ｈｘ、ｈｙ、ｈｚの長方体となり、その体積が次式となることが分かる。
【数３８】
【００４６】
ここで、Ｖは要素体積である。要素形状の対称性を考えると、ＣＶ１＝Ｖ／８(ｉ＝１，２，・・・，８)の結果が得られる。ＣＶの添え字で所属節点を表わす。一方、Ｑが要素内で定数の場合には、式(２４)右辺の要素積分は次式になる。
【数３９】
【００４７】
ここで、Ｑ(ｅ)は要素でのソース項である。式（２７）より節点１に寄与された熱量はＶ／８から発生したものであることが分かる。ほかの各節点への寄与量もＶＱ(ｅ)／８になることがGalerkin重み関数の定義から分かる。各節点のコントロール・ボリュームの体積とソース項の体積とが一致するため、保存則は２次精度で満たされている。
【００４８】
なお、図１７（ａ）に示すｙｚ平面上の平行四辺形をｘ方向に２ｈｘ分押し出して生成した六面体要素の場合には、有限要素法による式（２４）左辺の要素積分は次式になる。
【数４０】
【００４９】
式(２８)右辺の物理的な意味は、それぞれ広さがｈ２で中点Ｉ１２を通る面積、広さがｈｘｈで中点Ｉ１４を通る面積、広さが２ｈｈｘで中点Ｉ１８を通る面積、を通過する熱量である。コントロール・ボリュームは、検査面が閉じていないため厳密に求めることができないが、図１６(ｂ)に比べると、Ｖ／８以上になることが考えられる。一方、Ｑが要素内で定数の場合には、各節点のソース項の要素積分結果は依然としてＶ／８から発生した熱量になる。すなわち、節点に寄与された熱量の発生体積とコントロール・ボリュームの体積とが一致せず、保存則が満たされていないことが考えられる。
【００５０】
このような２次精度で保存則が満たされない問題は、ほかの一般形状の六面体要素にも存在すると考えられる。しかし、コントロール・ボリュームについて理論上で調べることが困難である。
【００５１】
［五面体要素の場合］
要素ｅを図１８（ａ）の線形正三角柱要素に当てはめた場合には、式（２４）左辺の要素積分は式（２９）となる。ただし、対象節点Ｐの要素における局所番号を１とする。
【数４１】
【００５２】
同様に、式（２９）右辺の物理的な意味について検討すると図１８（ｂ）に示す検査面が得られるため、節点１のコントロール・ボリュームは次式となる。
【数４２】
【００５３】
要素形状の対称性を考えると、ＣＶｉ（ｉ＝１，２，・・・，６）の結果が得られる。一方、Ｑが要素内で定数の場合には、式(２４)右辺の要素積分は次式になる。
【数４３】
【００５４】
すなわち、節点１に寄与された熱量はＶ／６から発生したものである。ほかの節点への寄与量も同じであることがGalerkin重み関数の定義から分かる。各節点のコントロール・ボリュームの体積とソース項の体積とが一致するため、保存則は２次精度で満たされている。
【００５５】
なお、図１９(ａ)に示す直角三角柱要素の場合には、有限要素法による式(２４)左辺の要素積分は次式になる。
【数４４】
【００５６】
式(３２)右辺の物理的な意味は、それぞれ広さがｈｙｈｚ／２で中点Ｉ１２を通る面積、広さがｈｘｈｚ／２で中点Ｉ１３を通る面積、広さがｈｘｈｙ／６で中点Ｉ１４を通る面積、を通過する熱量である。コントロール・ボリュームは、検査面が閉じていないため厳密に求められないが、図１８(ｂ)に比べるとＶ／６以上になることが考えられる。一方、Ｑが要素内で定数の場合には、各節点のソース項の要素積分結果は依然としてＶ／６から発生した熱量になる。節点に寄与された熱量の発生体積とコントロール・ボリュームの体積とが一致せず、保存則が満たされていないことが考えられる。
【先行技術文献】
【非特許文献】
【００５７】
【非特許文献１】O.C.Zienkiewicz，R.L.Taylor著（矢川元基ら訳）、「マトリクス有限要素法」（特に第１０章参照）、科学技術出版社、１９９６年
【非特許文献２】邵 長城著、「基本からわかる有限要素法」（特に第５章参照）、森北出版、２００８年
【非特許文献３】MING-DER DONALD HUANG, The Constant-Flow Patch Test−A Unique Guideline for the Evaluation ofDiscretization Schemes for the Current Continuity Equations. IEEE Transactions on computer-aided design, V.CAD-4, No. 4, 1985, pp.583-608.
【非特許文献４】Christian Cordes and Mario Putti, Accuracy ofGalerkin finite element for groundwater flow simulations in two andthree-dimensional triangulations. Int.J. Numer. Meth. Engng 2001, V52, P371-387.
【発明の概要】
【発明が解決しようとする課題】
【００５８】
以上説明したような問題から、本発明では以下の（ａ）〜（ｅ）を解決しようとする課題とする。
（ａ）上記従来の二次元三角形の問題の改善方法１において、鈍角三角形要素の場合には、ソース項の積分領域が要素外までに広がる問題を解決する。同時に、境界条件と線上分布ソース項の扱い方を提案する。
（ｂ）従来の有限要素法による三次元四面体の問題におけるソース項との不整合性問題を解決し、上記従来の三次元四面体の問題の改善方法における弊害を回避する。同時に、境界条件、線上分布ソース項、面上分布ソース項の扱い方を提案する。
（ｃ）上記従来の二次元四角形の問題におけるソース項との不整合性問題を解決する。
（ｄ）上記従来の三次元六面体・五面体の問題におけるソース項との不整合性問題を解決する。
（ｅ）有限要素法における荷重、体積力、質量の扱い方においても要素形状によって誤差が発生する問題も存在するため、その精度改善方法を提案する。
【００５９】
そこで本発明は、線形三角形要素、線形四面体要素、線形四角形要素、線形六面体要素、線形五面体要素を用いた有限要素法の離散化において、高精度のスキームを提案して、メッシュ・パターンによる解析誤差を低減し、例えば正三角形、正四面体、長方形、長方体などよりも品質が低い要素（品質の低いメッシュ・パターン）からも高精度な数値解が得られる有限要素法を用いた解析方法、及び有限要素法を用いた解析演算プログラムを提供することを目的とするものである。
【課題を解決するための手段】
【００６０】
本発明は、解析対象の解析領域を選定する選定工程（Ｓ１）と、
前記解析領域を計算対象として複数の要素に分割する分割工程（Ｓ２）と、
前記複数の要素のうちの、ある要素（例えばｅ）にある節点（例えばＰ，要素ｅでの節点番号を例えば１とする）についてGalerkin重み関数（Ｗ）をかけて要素積分をし、各要素のマトリクスを作成する要素マトリクス作成工程（Ｓ３）と、
Galerkin重み関数（Ｗ）と一般関数（例えばＱ，ｆ）との積からなる一般関数項を積分する一般関数項積分工程（Ｓ４）と、
前記ある節点（例えばＰ）周りの領域（例えばΩ）の要素における各要素のマトリクスの和と、前記一般関数項を積分した値の和とに基づき、連立方程式を作成する連立方程式作成工程（Ｓ５）と、
前記連立方程式に境界条件の導入をする境界条件導入工程（Ｓ６）と、
前記連立方程式を演算して数値解を得る演算工程（Ｓ７，Ｓ８）と、を備えた有限要素法を用いた解析方法において、
前記一般関数項積分工程（Ｓ４）にあって、Galerkin有限要素法による二階微分項の離散化結果を基に定義された節点領域の考え方を導入し、要素の代表値（例えばＱＧ，ＱＯ，ｆＧ，ｆＯ）を用いた一般関数項を積分することを特徴とする。
【００６１】
また、本発明は、前記一般関数項積分工程（Ｓ４）にあって、節点領域の大きさに合わせて一般関数を節点に分配することを特徴とする。
【００６２】
更に、具体的に本発明は、前記一般関数項積分工程は、前記一般関数項をソース項としたソース項積分工程（Ｓ４）であり、
前記ソース項積分工程（Ｓ４）にあって、要素の代表ソース値（例えば、幾何学中心での値や節点平均位置での値、要素平均値など）を用いたソース項を積分することを特徴とする。
【００６３】
また、具体的に本発明は、前記一般関数項積分工程は、前記一般関数項を荷重、体積力、質量（フォースと総称する）としたフォース項積分工程（Ｓ４）であり、
前記フォース項積分工程（Ｓ４）にあって、要素の代表フォース値（例えば、幾何学中心での値や節点平均位置での値、要素平均値など）を用いたフォース項を積分することを特徴とする。
【００６４】
また、詳細には本発明は、前記複数の要素は、二次元三角形の要素であり、
前記ソース項積分工程（Ｓ４）にあって、前記各要素における節点領域の面積をＮＤ、前記要素の代表ソース値をＱＧ、前記要素の面積をＳとした場合に、（ＮＤ−Ｓ／３）ＱＧを追加項としてソース項に追加することを特徴とする。
【００６５】
また、詳細には本発明は、前記複数の要素は、二次元三角形の要素であり、
前記フォース項積分工程（Ｓ４）にあって、前記各要素における節点領域の面積をＮＤ、前記要素の代表フォース値をfＧ、前記要素の面積をＳとした場合に、（ＮＤ−Ｓ／３）fＧを追加項としてフォース項に追加することを特徴とする。
【００６６】
また、詳細には本発明は、前記複数の要素は、三次元四面体の要素であり、
前記ソース項積分工程にあって、前記各要素における節点領域の体積をＮＤ、前記要素の代表ソース値をＱＧ、前記要素の体積をＶとした場合に、（ＮＤ−Ｖ／４）ＱＧを追加項としてソース項に追加し、かつ前記ＮＤとして以下の数式を用いることを特徴とする。
【数４５】
【００６７】
また、詳細には本発明は、前記複数の要素は、三次元四面体の要素であり、
前記フォース項積分工程（Ｓ４）にあって、前記各要素における節点領域の体積をＮＤ、前記要素の代表フォース値をfＧ、前記要素の体積をＶとした場合に、（ＮＤ−Ｖ／４）ＱＧを追加項としてフォース項に追加し、かつ前記ＮＤとして以下の数式を用いることを特徴とする。
【数４６】
【００６８】
また、詳細には本発明は、前記複数の要素は、二次元四角形の要素であり、
前記ソース項積分工程（Ｓ４）にあって、前記各要素における節点領域をＮＤ、前記Galerkin重み関数をＷP、前記要素の代表ソース値をＱＯとした場合に、以下の数式を追加項としてソース項に追加することを特徴とする。
【数４７】
【００６９】
また、詳細には本発明は、前記複数の要素は、二次元四角形の要素であり、
前記フォース項積分工程（Ｓ４）にあって、前記各要素における節点領域をＮＤ、前記Galerkin重み関数をＷP、前記要素の代表フォース値をfＯとした場合に、以下の数式を追加項としてフォース項に追加することを特徴とする。
【数４８】
【００７０】
また、詳細には本発明は、前記複数の要素は、三次元六面体又は三次元五面体の要素であり、
前記ソース項積分工程（Ｓ４）にあって、前記各要素における節点領域をＮＤ、前記Galerkin重み関数をＷP、前記要素の代表ソース値をＱＯとした場合に、以下の数式を追加項としてソース項に追加することを特徴とする。
【数４９】
【００７１】
また、詳細には本発明は、前記複数の要素は、三次元六面体又は三次元五面体の要素であり、
前記フォース項積分工程（Ｓ４）にあって、前記各要素における節点領域をＮＤ、前記Galerkin重み関数をＷP、前記要素の代表フォース値をfＯとした場合に、以下の数式を追加項としてフォース項に追加することを特徴とする。
【数５０】
【００７２】
また、本発明は、前記境界条件導入工程（Ｓ６）にあって、前記各要素の自然境界を通過する熱流束がゼロでない場合に、節点領域の範囲と合致した境界熱流束の節点への配分を行うことを特徴とする。
【００７３】
また、本発明は、前記境界条件導入工程（Ｓ６）にあって、前記各要素の境界での荷重、体積力、質量がゼロでない場合に、節点領域の範囲と合致した荷重、体積力、質量の節点への配分を行うことを特徴とする。
【００７４】
そして、本発明は、有限要素法を用いた解析演算をコンピュータに実行させるための有限要素法を用いた解析演算プログラムであって、
前記コンピュータに、
解析対象の解析領域を選定する選定工程（Ｓ１）と、
前記解析領域を計算対象として複数の要素に分割する分割工程（Ｓ２）と、
前記複数の要素のうちの、ある要素（例えばｅ）にある節点（例えばＰ，要素ｅでの節点番号を例えば１とする）についてGalerkin重み関数（Ｗ）をかけて要素積分をし、各要素のマトリクスを作成する要素マトリクス作成工程（Ｓ３）と、
Galerkin重み関数（Ｗ）と一般関数（例えばＱ，ｆ）との積からなり、Galerkin有限要素法による二階微分項の離散化結果を基に定義された節点領域の考え方を導入し、要素の代表値を用いた一般関数項を、積分する一般関数項積分工程（Ｓ４）と、
前記ある節点（例えばＰ）周りの領域（例えばΩ）の要素における各要素のマトリクスの和と、前記一般関数項を積分した値の和とに基づき、連立方程式を作成する連立方程式作成工程（Ｓ５）と、
前記連立方程式に境界条件の導入をする境界条件導入工程（Ｓ６）と、
前記連立方程式を演算して数値解を得る演算工程（Ｓ７，Ｓ８）と、を実行させることを特徴とする。
【００７５】
なお、上記カッコ内の符号は、図面と対照するためのものであるが、これは、発明の理解を容易にするための便宜的なものであり、特許請求の範囲の構成に何等影響を及ぼすものではない。
【発明の効果】
【００７６】
本発明によると、所定の重み関数と一般関数との積からなる一般関数項を積分する一般関数項積分工程にあって、Galerkin有限要素法による二階微分項の離散化結果を基に定義された節点領域の考え方を導入し、要素代表値（例えば幾何学中心での値や、節点座標平均位置での値、要素平均値など）を用いた一般関数項を積分し、節点領域の大きさに見合った値を取り入れることができる。これにより、メッシュ・パターンによる解析誤差を低減することができ、品質が低い要素（品質の低いメッシュ・パターン）からも高精度な数値解を得ることが可能となる。このため、メッシュ・パターンの品質を高める作業（メッシュ・パターンを正三角形、正四面体、正方形などに近づける作業）を不要とすることができるものでありながら、アルゴリズム（プログラム）の修正による計算量の増加が殆んどないので、コンピュータによる演算処理の増大も防止することができ、総じて解析時間の短縮化を図ることができる。
【図面の簡単な説明】
【００７７】
【図１】二次元三角形における節点Ｐの積分領域Ωを示す図。
【図２】線形三角形要素を示す図。
【図３】一般三角形要素での節点領域を示す図。
【図４】直角三角形要素での節点領域を示す図。
【図５】鈍角三角形要素での節点領域を示す図。
【図６】温度分布の等値線を示す図。
【図７】従来改善方法２による節点領域の修正を説明する図。
【図８】線形直角四面体要素を示す図。
【図９】四面体要素による解析結果を示す図で、（ａ）はメッシュを示す図、（ｂ）は（ａ）のある断面上の温度分布を示す図。
【図１０】従来改善方法による四面体要素での節点領域を示す図。
【図１１】二次元四角形における節点Ｐの積分領域を示す図。
【図１２】長方形要素での積分結果を示す図。
【図１３】平行四辺形要素での積分結果を示す図。
【図１４】二節点が重なる四辺形要素での積分結果を示す図。
【図１５】三次元形体における節点Ｐの積分領域を示す図。
【図１６】長方体要素を示す図で、(ａ)は要素形状を示す図、(ｂ)は積分結果を示す図。
【図１７】平行四辺形を押し出した六面体要素を示す図で、(ａ)は要素形状を示す図、(ｂ)は積分結果を示す図。
【図１８】正三角柱要素を示す図で、(ａ)は要素形状を示す図、(ｂ)は積分結果を示す図。
【図１９】直角三角柱要素を示す図で、(ａ)は要素形状を示す図、(ｂ)は積分結果を示す図。
【図２０】鋭角三角形要素での節点領域を示す図。
【図２１】直角三角形要素での節点領域を示す図。
【図２２】鈍角三角形要素での節点領域を示す図で、（ａ）は面積が正の部分を示す図、（ｂ）は面積が負の部分を示す図。
【図２３】自然境界条件及び線上に分布するソース項の処理を説明する図で、（ａ）は節点３の角度≦π／２の場合の図、（ｂ）は節点３の角度＞π／２の場合の図。
【図２４】四面体における線上分布を説明するための図。
【図２５】自重による変位問題を説明するための図。
【図２６】板の有限要素モデルを示す図。
【図２７】二次元解析メッシュのパターンを示す図で、（ａ）は直角三角形メッシュの図、（ｂ）は鈍角三角形メッシュの図、（ｃ）は細長い三角形メッシュの図。
【図２８】二次元熱伝導解析の誤差を、従来の有限要素法（ＧＦＥ）、改善方法１（ＮＤＦＥ）、本発明の方法（ＣＮＤＦＥ）について示す図で、（ａ）は直角三角形メッシュの場合の誤差を示す図、（ｂ）は鈍角三角形メッシュの場合の誤差を示す図、（ｃ）は細長い三角形メッシュの場合の誤差を示す図。
【図２９】荷重、体積力、質量が係る三角形メッシュと拘束条件を示す図。
【図３０】従来有限要素法を三角形メッシュに適用した内圧での径方向変位の解析結果を示す図。
【図３１】従来有限要素法を三角形メッシュに適用した回転での径方向変位の解析結果を示す図。
【図３２】本発明を三角形メッシュに適用した内圧での径方向変位の解析結果を示す図。
【図３３】三次元解析メッシュのパターンを示す図で、（ａ）は立方体を六個の四面体に分割したメッシュの図、（ｂ）は長方形を十四個の四面体に分割したメッシュの図、（ｃ）は（ｂ）の四面体をさらに四個の四面体に分割したメッシュの図。
【図３４】三次元熱伝導解析の誤差を、従来の有限要素法（ＧＦＥ）、本発明の方法（ＣＮＤＦＥ）について示す図で、（ａ）は立方体を六個の四面体に分割したメッシュの場合の誤差を示す図、（ｂ）は長方形を十四個の四面体に分割したメッシュの場合の誤差を示す図、（ｃ）は（ｂ）の四面体をさらに四個の四面体に分割したメッシュの場合の誤差を示す図。
【図３５】四角柱の四面体要素分割を示す図。
【図３６】従来有限要素法を四面体メッシュに適用した自重での変位の解析結果を示す図。
【図３７】本発明を四面体メッシュに適用した自重での変位の解析結果を示す図。
【図３８】Ｔ１を含む項の節点領域の定義方法を説明するための図。
【図３９】Ｔ１を含まない項の節点領域の定義方法を説明するための図で、（ａ）は正の面積を示す図、（ｂ）は負の面積を示す図。
【図４０】四辺形要素の節点位置を説明するための図。
【図４１】四辺形メッシュの種類を示す図で、（ａ）は正方形のメッシュ、（ｂ）は３，４，５，８要素からなるメッシュ、（ｃ）は菱型のメッシュ、（ｄ）は平行四辺形のメッシュ、（ｅ）台形のメッシュ、（ｆ）直角台形のメッシュ、（ｇ）直角三角形型台形のメッシュ、（ｈ）平角四辺形のメッシュ、（ｉ）低平角四辺形のメッシュ、（ｊ）高平角四辺形のメッシュ。
【図４２】基本境界条件での解誤差比較を示す図。
【図４３】基本・自然境界条件での解析誤差比較を示す図。
【図４４】基本境界条件での節点数と解析誤差の相関図。
【図４５】基本・自然境界条件での節点数と解析誤差の相関図。
【図４６】荷重、体積力、質量が係る四角形メッシュと拘束条件を示す図。
【図４７】従来有限要素法を四角形メッシュに適用した内圧での径方向変位の解析結果を示す図。
【図４８】本発明を三角形メッシュに適用した内圧での径方向変位の解析結果を示す図。
【図４９】節点１とｉの温度に関わる節点領域を示す図。
【図５０】係数Ｄｊｋの計算方法を説明する図。
【図５１】節点１以外の温度に関わる節点領域を示す図。
【図５２】六面体メッシュを示す図で、（ａ）は長方体の図、（ｂ）は６，８，１０，１６要素なるメッシュの図、（ｃ）は菱型柱の図、（ｄ）は平行四辺柱の図、（ｅ）は台形柱の図、（ｆ）は直角三角形型台形の図、（ｇ）平角柱の図、（ｈ）は低平角柱の図。
【図５３】六面体要素での解析誤差を示す図で、（ａ）は基本境界条件での解析誤差を示す図、（ｂ）は基本・自然境界条件での解析誤差を示す図。
【図５４】六面体要素での節点数と解析誤差の相関を示す図で、（ａ）は基本境界条件での解析誤差を示す図、（ｂ）は基本・自然境界条件での解析誤差を示す図。
【図５５】五面体メッシュを示す図で、（ａ）は直角三角柱の図、（ｂ）は鈍角三角柱の図、（ｃ）は細長い三角柱の図、（ｄ）は方向性を持つ直角三角柱の図。
【図５６】五面体要素での解析誤差を示す図で、（ａ）は基本境界条件での解析誤差を示す図、（ｂ）は基本・自然境界条件での解析誤差を示す図。
【図５７】本発明に係る有限要素法の解析工程を示すフローチャート。
【発明を実施するための形態】
【００７８】
以下、本発明に係る実施の形態を詳細に説明する。なお、以下の説明においては、大まかに、メッシュ・パターンが二次元三角形の場合と三次元四面体の場合と二次元四角形の場合とに分けて説明する。また、以下の説明においては、有限要素法の解析方法について説明するが、一般的にはこの解析方法は解析プログラム化され、コンピュータによって演算処理されるものであり、言い換えると、図５７に示す各工程は、コンピュータに有限要素法を用いた解析演算を実行させるための解析演算プログラムということになる。
【００７９】
まず、有限要素法を用いた解析方法の流れを、図５７に沿って大まかに説明する。有限要素法により解析対象を数値解析する際は、図５７に示すように、ステップＳ１において、解析対象（どのような形状のものであっても良い）の解析領域を選定し（選定工程）、ステップＳ２において、解析領域を計算対象として複数の要素に分割する（分割工程）。続いてステップＳ３において、要素ｅの各節点（二次元の線形三角形の場合は節点１，２，３、三次元の線形四面体の場合は節点１，２，３，４、二次元の線形四角形の場合は節点１，２，３，４）について重み関数Ｗをかけて要素積分をし、各要素のマトリクス（例えば二次元三角形の場合は上述した式（３）の左辺、三次元四面体の場合は式（１０）の左辺、二次元四角形の場合は上述した式（３）’の左辺や式（１５）の左辺、三次元六面体・三次元五面体の場合は式（２４）の左辺）を作成する（要素マトリクス作成工程）。
【００８０】
次に本発明の要部となるステップＳ４において、Galerkin重み関数Ｗと一般関数（例えば熱問題等であればＱ、荷重、体積力、質量問題であればｆ）との積からなる一般関数項を積分する（一般関数項積分工程）。このステップＳ４において例えば熱問題等を扱うときは、ポアソンの方程式のソース項を積分する（例えば二次元三角形の場合は後述する式（３３）の右辺、三次元四面体の場合は後述する式（３６）の右辺、二次元四角形の場合は後述する式（４６）の右辺、三次元六面体・五面体の場合は後述する式（７２）の右辺）（ソース項積分工程）。また、例えば荷重、体積力、質量が係る場合には、二次元三角形の場合は式（３３）’、三次元四面体の場合は式（３６）’、二次元四角形の場合は式（７１）、三次元六面体・五面体の場合は式（９６）を用いて、荷重、体積力、質量の要素積分、即ち、それぞれ式（３３）’の一項目、式（３６）’の一項目、式（７１）の一項目、式（９６）の一項目を取って代わる（フォース項積分工程）。
【００８１】
その後、ステップＳ５において、節点Ｐ周りの要素領域Ωにおける各要素のマトリクスの和と、一般関数項を積分した値の和とに基づき、連立方程式を作成する（連立方程式作成工程）。
【００８２】
この際、本発明では、詳しくは後述するように、ステップＳ６において、節点領域の範囲と合致した境界条件を、各節点周りについて作成した連立方程式のそれぞれに導入する（境界条件導入工程）。
【００８３】
そして、ステップＳ７において、連立方程式を例えば反復法等により演算して解き（演算工程）、ステップＳ８において、関連物理量の計算を行って数値解を得て（演算工程）、以上の工程によって解析を終了する。
【００８４】
［二次元三角形の解析について］
以上のような有限要素法を用いた解析方法において、本発明にあっては、ステップＳ４の一般関数項積分工程における一般関数項を、要素ｅの代表値（例えば幾何学中心での値や、節点座標平均位置での値、要素平均値など）を用いて一般関数項について計算する。
【００８５】
具体的には、ポアソンの方程式の場合には、式（１）左辺の楕円型演算子については従来の有限要素法により離散化をし、右辺のソース項については、従来の有限要素法による離散化
【数５１】
に追加項
【数５２】
を加える。すなわち、従来の式（３）の代わりに、式（３３）を提案する。
【数５３】
【００８６】
ここで、ＱＧは要素のソース代表値であり、ＮＤは対象節点の要素ｅにおける節点領域の面積である。ＮＤの定義は、鋭角三角形要素においては図２０に、直角三角形要素においては図２１に、鈍角三角形要素においては図２２に、示されるように行う。図中のＮＤの添え字で所属節点を表わし、点Ｃは外心であり、点Ｉ、Ｊ、Ｋは中点である。また、鈍角三角形要素の場合には、外心が要素の外に出ているため、図２２（ｂ）に示す節点２、節点３が所有する節点領域−ＮＤ２、−ＮＤ３が負になる。それぞれの節点領域は最終的に次式により算出される。
【数５４】
【００８７】
また、ステップＳ６の境界条件導入工程にあっては、以下のように境界条件を導入する。即ち、熱流束がゼロではない自然境界については図２３に示すように処理する。ここで、節点１と節点２を結ぶ線を境界とすると、節点３での内角がπ／２以下の場合（つまり鋭角の場合）には、図２３（ａ）のように点１と点Ｉとの間を通過し流出する熱量
【数５５】
を節点１の方程式に、点Ｉと点２との間の熱通過量
【数５６】
を節点２の方程式に加える。ここで、
【数５７】
は線積分であり、ｎは外向き法線方向である。
【００８８】
節点３での内角がπ／２を越えた場合（つまり鈍角の場合）には、図２３（ｂ）のように点１と点Ｉ１間の熱通過量
【数５８】
を節点１の方程式に、点Ｉ２と点２との間の熱通過量
【数５９】
を節点２の方程式に、点Ｉ１と点Ｉ２との間の熱通過量
【数６０】
を節点３の方程式に加える。
なお、線上に分布するソース項に対しても、図２３に示す自然境界条件の分け方と同様に節点の方程式に寄与する。
【００８９】
以上説明した本実施の形態による有限要素法を用いた二次元三角形の解析方法によれば、以下の特徴を持っている。
（１）節点領域の大きさに見合ったソース量が取り入れられたため、節点レベルで保存法則を２次精度で満足する。
（２）要素の品質が低いほど修正効果が大きい。
（３）いかなる要素形状においても、要素外のソース項を取り入れることがない。
（４）いかなる要素形状においても、要素の３つの節点が所有する節点領域の合計は
【数６１】
になるので、ソース項に追加項を入れても、要素レベルでも保存法則が満たされる。
（５）楕円型演算子に対する最適な離散化手法としての有限要素法の優位性がそのまま保持される。
（６）アルゴリズムの修正による計算量の増加は殆どない。
【００９０】
［三次元四面体の解析について］
また、三次元四面体の解析にあっても、前述のような有限要素法を用いた解析方法において、ステップＳ４の一般関数項積分工程における一般関数項を、要素ｅの代表値を用いて一般関数項について計算する。
【００９１】
具体的には、式（８）左辺の楕円型演算子については従来の有限要素法により離散化をし、右辺のソース項については追加項
【数６２】
を加えて、つまり従来の式（１０）の代わりに、式（３６）を提案する。
【数６３】
【００９２】
ここでも、ＱＧはソース代表値であり、ＮＤは対象節点の要素ｅにおける節点領域の体積である。節点領域ＮＤは以下の方法で決める。対象節点の要素ｅにおける節点番号を１とし、そこに局所座標系の原点を置くと、式（３６）左辺の積分結果を次式右辺のような係数と熱流束との掛け算の形に整理することができる。
【数６４】
【００９３】
ここで、
【数６５】
は図１０に示す要素の辺中点Ｉ、Ｊ、Ｋでの節点１から節点２、３、４に向かう熱流束を二次精度で近似するものであり、係数Ｓ１２、Ｓ１３、Ｓ１４はそれぞれの熱流束の通過面積にあたるものであり、その二つの添え字で熱流束の向きを表わす。節点１の節点領域の体積は、熱流束の通過面積を底面とし、節点１と中点間の距離を高さとする角錐の体積の合計として、次式のように求める。
【数６６】
【００９４】
また、熱流束の通過面積Ｓｉｊ（ｉ，ｊ＝１，２，３，４、ｉ≠ｊ）を表わす汎用式は、以下のように導出される。まず、式（３７）の左辺を有限要素法のアプローチにより展開すると次式が得られる。
【数６７】
【００９５】
ここで、要素体積Ｖは次式のように求められる。
【数６８】
【００９６】
式（３９）中の係数は以下のようなものである。
【数６９】
【００９７】
式（３９）よりλ（Ｔｉ−Ｔｊ）の係数は
【数７０】
になることが分かる。節点ｉとｊとの間の距離をｌｉｊで表わすと、式（４２）を得ることができ、λ（Ｔｉ−Ｔｊ）の係数を式（４３）のように書き換えることができる。
【数７１】
【数７２】
【００９８】
式（３７）の右辺を参照すると、節点ｉと節点ｊとの間に伝わる熱流束の通過面積は次式により求めることができる。
【数７３】
この式（４４）を用いて、各節点の節点領域を式（３８）と同じ考え方で算出する。
【００９９】
また、ステップＳ６の境界条件導入工程にあっては、以下のように境界条件を導入する。即ち、熱流束がゼロではない自然境界では、従来の有限要素法による処理結果
【数７４】
に、追加項
【数７５】
を加えて、
【数７６】
を境界条件として式（３６）に導入する。ここで、ｓは要素ｅに所属する境界を表わし、Ｓはｓの面積であり、添え字ＳＧはｓの幾何中心を意味し、ｎは外向き法線方向で、ＮＤはｓを上述の二次元三角形解析の方法で分割した節点領域の面積である。
【０１００】
なお、線上に分布するソース項Ｑに対しても、まず、線を共有する各要素の表面三角形に均等に分ける。例えば、図２４に示す節点１と節点２とを結ぶ線の場合には、それを有する要素表面三角形は、三角形１２４、１２３、１２３、１２５と四つがある。三角形１２３は二要素に属するため二つとして数える。各三角形にＱ／４を配分する。さらに、三角形内の三節点への配分については、上述した二次元三角形解析における線上に分布するソース項の分け方と同方法で行う。面上に分布するソース項についても、二次元三角形解析の方法で同様に配分する。
【０１０１】
以上説明した本実施の形態による有限要素法を用いた三次元四面体の解析方法によれば、以下の特徴を持っている。
（１）節点領域の大きさに見合ったソース項が分配されるため、節点レベルにおいては保存法則を従来法より高い精度で満たす。
（２）要素の品質が低いほど修正効果が大きい。
（３）任意の四面体形状においても、要素外のソース項を取り入れることがない。
（４）任意の四面体形状においても、要素の４節点に所属する節点領域の合計が要素体積に等しいため、追加項を入れても要素レベルで保存法則が満たされる。
（５）楕円型演算子に対する最適な離散化手法としての有限要素法の優位性がそのまま保持される。
（６）アルゴリズムの修正による計算量の増加は殆どない。
【０１０２】
［二次元三角形解析及び三次元四面体解析における他分野の解析］
なお、上述した二次元三角形解析と三次元四面体解析とでは、熱伝導問題を例に本解析手法の応用について説明してきたが、勿論、本手法はポアソンの方程式が記述する現象の全てにそのまま展開できる。例えば、流れ問題、電気・磁気問題、弾性変形問題、物質拡散問題などが応用できる分野として挙げられる。
【０１０３】
また、本発明の基本的な考え方は、二次精度の保存則の観点から楕円型演算子の離散化結果におけるコントロール・ボリュームを見出し、ソース項の積分領域をそれに合わせたものである。従って、この考え方は、保存則から導出された楕円型演算子を含む各種の支配方程式の有限要素解析に容易に応用できる。すなわち、物理量の発生・消失と流出・流入を算出するとき、その積分領域について、本発明で導出した節点領域の考え方を導入して積分すればよい。例えば、非定常移流拡散解析において、非定常項と移流項の離散化に対して本発明におけるソース項の積分領域の考え方が適用できる。
【０１０４】
また、荷重、体積力と質量（フォースと総称する）が係わる従来の有限要素法による解析においては、それらの扱い方はポアソンの方程式のソース項と同じである（非特許文献２参照）（すなわち、Galerkin重み関数とかけて要素積分をする。）ため、従来の技術で述べたような、要素形状と関係なく均等的に節点方程式に分配される問題も存在するが、これらの問題は、本実施の形態で説明したソース項の扱い方で改善することもできる。
【０１０５】
従来の技術での問題点を以下の例で説明する。図２５に示す縦に置かれた幅がＬ、高さがＨ、厚さが１の板の自重による変位について考察する。密度をρ、重力加速度をｇ、弾性係数をＥで表わす。この際、問題の本質をより明確にするため、平面応力問題とし、ポアソン比をν＝０とする。板全体の高さの変化量をＶで表わし、上向きの座標をｙとすると、応力がσ＝−（Ｈ−ｙ）ρｇとなり、ひずみがε＝σ／Ｅとなるため、次式の厳密解が得られる。
【数７７】
【０１０６】
一方、板を図２６のように２つの三角形要素に分割する。節点３と節点４に寄与されるｙ方向のフォース量をそれぞれｆ３とｆ４で表わすと、従来の有限要素法では、式（５）に示すように要素の重さの３分の１ずつが各節点に分配されるため、
【数７８】
となる。ここで、節点４が二つの要素を所有するため、分配された荷重は節点３の２倍になっている。解析結果として、節点３と節点４のｙ方向の変位ν３とν４が次のようになる。
【数７９】
これらの結果は、厳密解からずれていることが確認できる。
【０１０７】
［線上に分布する場合］
二次元三角形解析においては、図２３に示すように、節点３での内角がπ／２以下の場合には、
【数８０】
をそれぞれ節点１、２の方程式に加える。ここで、ｆは荷重の任意成分、体積力の任意成分、質量を意味する。節点３での内角がπ／２を越えた場合には、
【数８１】
をそれぞれ節点１、２、３の方程式に加える。
【０１０８】
三次元四面体解析においては、まず、線を共有する各要素の表面三角形に均等に分ける。例えば、図２４に示す節点１と節点２を結ぶ線の場合には、四つの表面三角形に、それぞれｆ／４を配分する。さらに、三角形内の三節点への配分については、上述した二次元三角形解析における線上に分布するソース項の分け方と同方法で行う。
【０１０９】
［面上に分布する場合］
二次元三角形、あるいは三次元四面体解析において、面ｓ上に分布する荷重、体積力、質量については、式（３３）右辺のソース項と同様に修正し、修正後の対象節点に寄与するものは
【数８２】
になる。ここで、Ｓはｓの面積であり、添え字Ｇで要素代表値（例えば幾何学中心での値や、節点座標平均位置での値、要素平均値など）を表わす。
【０１１０】
［空間に分布する場合］
空間に分布する荷重、体積力、質量については、式（３６）右辺のソース項と同様に修正し、修正後のものは
【数８３】
になる。
【０１１１】
［二次元三角形解析及び三次元四面体解析における本発明の効果］
本解析方法の導入により解析精度が改善されたことを以下の事例により示す。
【０１１２】
前記板の自重による変位問題に、節点領域の考え方を適用し、式（３３）’と図２１を参照すると、節点３と節点４に寄与されるｙ方向のフォース量は次のようになることが分かる。
【数８４】
【０１１３】
ここで、節点４に寄与されたフォース量は、二つの要素の合計になる。解析結果として、節点３と節点４のｙ方向の変位は次のようになる。
【数８５】
厳密解と一致しており、解析精度が改善されたことが確認できる。
【０１１４】
二次元三角形問題については、式（６）と（７）により定義された熱伝導問題を図２７のメッシュ・パターンを用いて解析した最大誤差を図２８に示す。この際λ＝１とし、ソース項の要素代表値として幾何学中心を採用すると、最大誤差は次式により定義される。
【数８６】
【０１１５】
ここで、ＴＡｉ：節点ｉでの厳密解、ＴＮｉ：節点ｉでの数値解、ＴＡｍａｘ：厳密解の最大値、である。図２８の横軸はメッシュの代表長さであり、同じパターンのメッシュを２倍ずつ拡大したものである。中には、図２８（ａ）の１／８、（ｂ）の１／８、（ｃ）の１／５が、図２７のメッシュそのままを使用したものである。「ＧＦＥ」で従来の有限要素法、「ＮＤＦＥ」で上記改善方法１、「ＣＮＤＦＥ」で本発明、での結果を示す。図２８より、本提案方法は、いずれのメッシュ・パターンにおいても解析誤差が従来の有限要素法の３分の１以下になったことが分かる。また、上記改善方法１は、直角三角形メッシュの場合にはよい精度が得られたが、細長いメッシュを使用すると、誤差が著しく増大することが認められた。
【０１１６】
また、荷重、体積力と質量が係わる有限要素法による解析の事例として、薄い円筒に、（１）均等な内圧を加えた変形問題、（２）円筒中心軸を回転中心とした等速回転問題、を取り上げた。図２９に解析メッシュを示しており、直角等辺三角形のシェル要素によりメッシュ分割を行った。この際、軸方向の中央に位置する節点（●で示す節点）に円周方向と軸方向の拘束条件を与えた。
【０１１７】
図３０に従来の有限要素法を用いた内圧による径方向の変位を、図３１に従来の有限要素法を用いた回転による径方向の変位を示す。点線で初期形状を示し、灰色で膨張方向の変位を、黒色で縮小方向の変位を表わす。各節点は同じ変位をする実現象と異なって、節点周辺の局所メッシュ・パターンによって、変位差が現れた。その原因は、図６の説明で述べた不整合性の問題と共通である。図３２に本発明での提案方法（上記「面上に分布する場合」）による内圧での径方向の変位結果を示しており、要素形状に応じた節点荷重を与えた結果、均一変位になっていることが確認できる。回転問題についても、解析原理が同じであるため、遠心力を同じ考え方で各節点に分配すれば、均一変位が得られることは明確的である。
【０１１８】
三次元四面体問題については、式（１３）と（１４）により定義された熱伝導問題を図３３のメッシュ・パターンを用いて解析した最大誤差を図３４に示す。この際λ＝１とし、ソース項の要素代表値として幾何学中心値を採用した。メッシュ・パターンとして、図３３（ａ）は一個の立方体を六個の四面体に、（ｂ）は一個の長方体を十四個の四面体に、（ｃ）は（ｂ）の四面体をさらに四個の四面体に分割したものである。本発明の解析方法は、いずれのメッシュ・パターンにおいても解析精度が改善され、特に品質の悪い要素では誤差が従来の有限要素法の２分の１以下になったことが分かる。
【０１１９】
荷重、体積力と質量が係わる有限要素法による解析の事例として、図３５の四角柱の自重による変形問題を取り上げ、本発明の改善効果を示す。図３６には従来の有限要素法による高さ方向の変位結果を、図３７には本発明、すなわち、式（３６）’を導入した高さ方向の変位結果を、示す。上面における４節点が同じ変位をする実現象に対して、本発明の解析結果においては、４節点の変位差は、従来の有限要素法に比べ大幅に改善されたことが確認できる。改善の理由については、上面一番左の節点を例に説明する。従来の有限要素法においては、当節点に要素自重の４分の１しか寄与されていないため、変位の絶対値が小さくなっている。本発明では、当節点に式（３６）’に基づいて要素自重の２分の１が寄与されるため、変位の改善効果が得られた。
【０１２０】
［二次元四角形の解析について］
ついで、二次元四角形の解析について説明する。二次元四角形の解析において保存則を満たさない問題を改善するため、簡単に実行できるソース項の保存型離散化手法を提案する。式（３）’の左辺については有限要素法のままで離散化をし、右辺については追加項
【数８７】
を加えて、次式により有限要素解析を行う。
【数８８】
【０１２１】
ここで、ＮＤＰは節点Ｐの要素ｅにおける節点領域(Nodal
Domain)と称するものであり、ＱＯは点Ｏでのソース値である。ＮＤＰが式（１５）左辺の積分から導かれたコントロール・ボリュームと面積上で等しくなれば、保存則を２次精度で満足することができ、メッシュ依存性を改善し、特に、低品質要素での精度向上効果が期待できる。考え方としては、ＱＯの代わりに要素ソース項のほかの代表値、例えば重心でのソース値や要素平均値などを式（４６）に入れることも考えられる。また、式（３）’右辺の要素積分を単純にＮＤＰＱＧで近似する選択もある。
【０１２２】
一般形状の四辺形要素において、節点領域を以下の普遍性のある方法により算出する。節点Ｐの要素ｅにおける局所番号を１とすると、式（１５）左辺の要素積分の結果を次の汎用式に整理することができる。
【数８９】
係数Ａｉは四節点の座標のみで決まる定数である。有限要素法の性質より、次式が成り立つ。
【数９０】
【０１２３】
［二次元四角形における節点領域の算出方法１］
ｉ，ｊ，ｋで、それぞれ節点２、３、４の中の異なる一節点を表わし、順序には制限がないこととする。係数Ａ２，Ａ３，Ａ４が持つ符号の全ての組み合わせと式（４８）を考えると、式（４７）を次の節点温度差の形に書き換えることができる。
【数９１】
【０１２４】
すなわち、Ａ１のみが正の場合には式（４９−１）を、Ａ１，Ａ２，Ａ３，Ａ４の中の二つの係数が正の場合には式（４９−２）、もしくは（４９−３）を、三つの係数が正の場合には式（４９−４）を採用する。また、式（４９−１）〜（４９−４）の右辺の各項の係数が正の値を持つことが分かる。
【０１２５】
式（４９−１）〜（４９−４）右辺の各項を、Ｔ１を含むものと含まないものに分類しておく。Ｔ１を含む項の展開については、λＢ１ｊ（Ｔ１−Ｔｊ）を例に説明する。ここで、Ｂ１ｊでλ（Ｔ１−Ｔｊ）項の係数を表わす。Ｌ１ｊで節点1と節点ｊ間の線分を表わすと、次式が得られる。
【数９２】
【０１２６】
式（５０）右辺は線分Ｌ１ｊの中点での節点1から節点ｊに向かう２次精度の熱流束により、Ｌ１ｊの中点から垂直方向に長さＢ１ｊＬ１ｊを伸びた線分を通過する熱量であることが分かる。この線分をＢ１ｊＬ１ｊで表わす。保存則の観点から、この項に相当するソース項の領域は、節点1、線分Ｌ１ｊの中点、線分Ｂ１ｊＬ１ｊの端点、といった三点から成す三角形となる。この三角形の面積をＮＤ１ｊで表わし次式により計算することができる。
【数９３】
【０１２７】
図３８には二つの例としてＮＤ１３とＮＤ１４を示す。Ｔ１を含まない項の展開については、λＢｉｊ（Ｔｉ−Ｔｊ）を例に説明する。節点ｉと節点ｊ間の線分をＬｉｊで表わすと、次式が得られる。
【数９４】
【０１２８】
式（５２）右辺は線分Ｌｉｊの中点での節点ｉから節点ｊに向かう２次精度の熱流束により、Ｌｉｊの中点から垂直方向に長さＢｉｊＬｉｊを伸びた線分を通過する熱量であることが分かる。この線分をＢｉｊＬｉｊで表わす。この項に相当するソース項の領域は、節点1、Ｌｉｊの中点（ｍと記す）、線分ＢｉｊＬｉｊのもう一つの端点（ｌと記す）といった三点から成す三角形となる。この三角形の面積をＮＤｉｊで表わし、次式により計算する。
【数９５】
式（５３）右辺の符号については、点1から中点ｍに向かうベクトルと流束ベクトルとの角度αが｜α｜＜π／２にある、すなわち、熱流束が三角形から流出する場合には正とする。αが｜α｜≧π／２にある、すなわち、熱流束が三角形に流入する場合には負とする。図３９には二つの例としてＮＤ２３＞０とＮＤ２４≦０を示す。
【０１２９】
最終的に、対象節点１の節点領域の面積は、係数Ａ２，Ａ３，Ａ４の符号に合わせて式（５４−１）〜（５４−４）中の一式により算出する。
【数９６】
【０１３０】
［二次元四角形における節点領域の算出方法２］
係数Ａ２，Ａ３，Ａ４の符号と大きさに関係なく、式（４９−１）、（５０）、（５１）と（５４−１）によりＮＤ１を求める。この際、ＮＤ１ｊ（ｊ＝２，３，４）はＢ１ｊ、すなわち−Ａｊと同じ符号になる。
【０１３１】
［本二次元四角形の解析手法の整合性、唯一性、正規性］
整合性とは、節点領域とコントロール・ボリュームと同面積であることを意味する。整合性が成立すれば、節点レベルでの保存則が２次精度で満足されることになり、精度改善効果が理論上で保証される。
【０１３２】
図１２の長方形要素に対しては、式（１６）を次式に書き換えておく。
【数９７】
【０１３３】
上記算出方法１の整合性については、ｈｙ／√２≦ｈｘ≦√２・ｈｙの場合には、Ｔ１以外の係数が負になっているため、式（４９−１）に当てはめ、次式が得られる。
【数９８】
【０１３４】
その結果、節点領域として次式が得られる。
【数９９】
【０１３５】
また、ｈｘ＞√２・ｈｙの場合には、Ｔ２の係数が正に、Ｔ３，Ｔ４の係数が負になる。Ａ１＞｜Ａ３｜であるため、式（４９−２）に当てはめると次式が得られる。
【数１００】
【０１３６】
その結果、節点領域として次式が得られる。
【数１０１】
【０１３７】
なお、ｈｙ＞√２・ｈｘの場合についてもＮＤ１＝Ｓ／４であることが証明できるため、ＮＤ１＝ＣＶ１の関係式がｈｙ／ｈｘの比例に関係なく成立していることが分かる。
【０１３８】
上記算出方法２の整合性については、式（４９−１）を採用したので、その結果は式（５７）となる。
【０１３９】
ほかの節点については、要素形状の対称性を考慮すると、ＮＤ２＝ＮＤ３＝ＮＤ４＝Ｓ／４を得ることができる。結局、長方形要素においては、式（４６）右辺の追加項がゼロとなり、保存型離散化がGalerkin有限要素法と同じものになる。
【０１４０】
さらに、図１３の平行四辺形要素、図１４の二節点が重なる四辺形要素に適用すると、上記表２に示す節点領域が得られる。以上の三種類の四辺形要素に対して、コントロール・ボリュームと完全一致していることが確認できる。しかし、一般形状の四辺形要素については、コントロール・ボリュームが明確になっていないため、代わりに後述の数値解析により精度改善効果について検証する。
【０１４１】
節点領域の唯一性とは、式（５４−１）〜（５４−４）により求められたＮＤ１の結果は、「ｉ，ｊ，ｋへの局所節点番号の与え方によって変らない」、「算出方法２を使っても変らない」、ことを意味する。例えば、ｈｘ＞√２・ｈｙの長方形要素に対して、式（５８）においては、ｉ＝２、ｊ＝３、ｋ＝４としたが、代わりにｉ＝２、ｊ＝４、ｋ＝３の選択も、次式により節点領域を求めることも考えられる。
【数１０２】
【０１４２】
その結果、次式の節点領域が得られる。
【数１０３】
これは、式（５９）の結果と同じであるため、算出方法１においては唯一性が証明された。算出方法２における唯一性は、式（５７）の結果により示されている。
【０１４３】
正規性とは、
【数１０４】
すなわち、式（４６）右辺に入れた追加項の同要素四節点についての合計がゼロになることを意味する。正規性が成立すれば、要素レベルでも、全解析領域でも保存則が満たされることが保証される。図１２、図１３、図１４の要素形状に対しては、上記表２の最後の行に示すように、
【数１０５】
になるため、正規性が保たれている。
【０１４４】
同様に図１３、図１４の要素形状についても、唯一性が保たれたことを理論的に証明できる。しかし、一般形状の四辺形要素においては、ヤコビアンが絡んでくるため、正規性と唯一性について理論上の議論は困難である。代わりに、下記の数値計算により検討する。四辺形要素の形状としては、節点１を図４０の太い●印の二通りの配置に対して、節点２を○印に置き、節点３を最下行以外の全点に置き、節点４を中心列を含めた左側の全点に置いた組み合わせで、できた各四辺形要素をテスト対象とした。これらの要素形状では実用解析で使用される要素形状をカバーできるといえる。この際、上記算出方法１においてはｉ，ｊ，ｋに許容される全ての局所節点番号の組み合わせを与えて数値解析により節点領域を算出した。その結果、面積がゼロ以下の形状と凹形状を除いて、上記算出方法１での全テストケースと上記算出方法２での全テストケースにおいて正規性と唯一性が満たされていることが確認できた。
【０１４５】
［二次元四角形における本発明の効果］
保存型離散化手法の精度改善効果を二つの熱伝導問題の解析事例により示す。この際、解析領域を正方形領域（０≦ｘ＜１、０≦ｙ≦１）に限定し、要素代表ソース値を四節点の平均位置での値とし、λ＝１とした。
【０１４６】
解析問題１としては、式（６３）のソース項に式（６４）の基本境界条件のみを加えたものである。
【数１０６】
【数１０７】
【０１４７】
その厳密解は次式になる。
【数１０８】
【０１４８】
解析問題２としては、式（６６）のソース項に式（６７）の基本境界条件と式（６８）の自然境界条件を加えたものである。
【数１０９】
【数１１０】
【数１１１】
【０１４９】
その厳密解は次式になる。
【数１１２】
【０１５０】
解析の数値誤差について次式により評価する。
【数１１３】
ここで、ＴＡｉ：節点ｉでの厳密解、ＴＮｉ：節点ｉでの数値解、ＴＡｍａｘ：厳密解の最大値、である。
【０１５１】
図４１に使用した１０パターンのメッシュを、対称性があるため０≦ｘ≦０．５、０≦ｙ≦０．５の部分のみ示す。「・」印で節点位置を、各図の表題の括弧内の数字で総節点数を表わす。図４２に基本境界条件での数値誤差、図４３に基本・自然境界条件での数値誤差を示す。これらの図より、本発明の手法（CSFEで表わす）では、正方形以外の全てのメッシュにおいて解析誤差が、従来の有限要素法（GFEで表わす）の２分の１以下に低減されたことが分かる。
【０１５２】
図４４と図４５に各メッシュにおける総節点数の逆数と解析誤差との相関を示しており、実線で二通りの正方形要素での結果を結んだ。これらの図より、従来の有限要素法にはメッシュ品質の影響が大きいが、保存型離散化手法にはメッシュ品質の影響が少なく、低平角四辺形を除けば、正方形とほぼ同程度の解析精度が得られた。
【０１５３】
結論として、線形四辺形要素によるPoissonの方程式の有限要素解析において、新たな節点領域の定義を考案し、２階微分項との整合性がとれたソース項の保存型離散化手法を提案することができた。従って、下記の（１）〜（４）の効果を得ることができた。
（１）節点方程式に寄与するソース量が、要素形状に関係する節点領域の大きさに対応させたことにより、節点レベルで保存則を２次精度で満足するようになった。しかも、要素レベルでも保存則が満たされている。
（２）鈍角要素形状においても、要素外のソース項を取り入れることがなく、ソース項が急激に変化する場合にも数値誤差を抑えることができた。
（３）要素の品質が低いほど、本手法の精度改善効果が大きい。
（４）幾つかのメッシュ・パターンを用いて数値解析を行った結果、従来の有限要素法に比べ解析精度が大幅に改善されることが確認できた。
【０１５４】
［二次元四角形解析における他分野の解析］
以上の二次元四角形解析においては、熱伝導問題を例に本発明の実施形態について説明してきたが、勿論、Poissonの方程式が記述する現象の全てにそのまま展開できる。例えば、流れ問題、電・磁気問題、弾性変形問題、物質拡散問題などが応用できる分野として挙げられる。
【０１５５】
また、本発明の基本的な考え方は、二次精度の保存則の観点から楕円型演算子の離散化結果におけるコントロール・ボリュームを見出し、ソース項の積分領域をそれに合わせたものである。従って、この考え方は、保存則から導出された楕円型演算子を含む各種の支配方程式の有限要素解析に容易に応用できる。すなわち、物理量の発生・消失と流出・流入を算出するとき、その積分領域について、本発明で導出した節点領域の考え方を導入して積分すればよい。例えば、非定常移流拡散解析において、非定常項と移流項の離散化に対して本発明におけるソース項の積分領域の考え方が適用できる。
【０１５６】
また、荷重、体積力、質量が係わる従来の有限要素法による解析においては、それらの扱い方はPoissonの方程式のソースと同じであるため、要素形状と関係なくほぼ均等的に節点方程式に分配される問題も存在する。これらの問題は本発明でのソース項の扱い方で改善することもできる。二次元（面上に分布する場合）、あるいは三次元問題（空間に分布する場合）の解析において、面ｓ上に分布する荷重、体積力、質量については、式（４６）の右辺のソース項と同様に追加項を加えて修正すると、対象節点に寄与する荷重、体積力、質量は
【数１１４】
となる。ここで、ｆは荷重の任意成分、体積力の任意成分、質量を意味する。
【０１５７】
荷重、体積力と質量が係わる有限要素法による解析の事例として、薄い円筒に、均等な内圧を加えた変形問題を取り上げた。図４６に解析メッシュを示しており、四角形のシェル要素によりメッシュ分割を行った。この際、●で示す節点に円周方向と軸方向の拘束条件を与えた。
【０１５８】
図４７に従来の有限要素法を用いた内圧による径方向の変位を示す。点線で初期形状を示し、灰色で膨張方向の変位を、黒色で縮小方向の変位を表わす。各節点は同じ変位をする実現象と異なって、節点周辺の局所メッシュ・パターンによって、変位差が現れた。図４８に本発明での提案方法（式（７１））による内圧での径方向の変位結果を示しており、要素形状に応じた節点荷重を与えた結果、変位のバラツキが大幅に改善されたことが確認できる。体積力と質量が係わる問題は、同じ原理が適用できるものである。
【０１５９】
［三次元六面体・五面体の解析について］
前述の「従来の有限要素法による三次元六面体及び三次元五面体の問題」で述べた保存則を満たさない問題を改善するため、簡単に実行できるソース項の保存型離散化手法を提案する。式（２３）の左辺については従来の有限要素法のままで離散化をし、右辺については追加項
【数１１５】
を加えて、次式により有限要素解析を行う。
【数１１６】
【０１６０】
ここで、ＮＤＰは節点Ｐの要素ｅにおける節点領域(Nodal Domain)と称するものであり、ＱＯは当要素に属する全節点の座標平均位置でのソース値である。ＮＤＰが式（２４）の左辺の要素積分から導かれたコントロール・ボリュームと体積上で等しくなれば、保存則を２次精度で満足することができ、メッシュ依存性を改善し、特に、低品質要素での精度向上効果が期待できる。また、考え方としてはＱＯの代わりに要素ソース項のほかの代表値、例えば重心でのソース値などを式（７２）に入れることも考えられる。また、式（２３）の右辺の要素積分を単純にＮＤＰＱＯで近似する選択もある。
【０１６１】
一般形状の六面体要素と五面体要素において、節点領域を以下の方法により算出する。節点Ｐの要素ｅにおける局所番号を１とすると、式（２４）左辺の要素積分の結果を次の汎用式に整理することができる。
【数１１７】
ここで、Ｎは一要素の節点数であり、六面体要素ではＮ＝８、五面体要素ではＮ＝６となる。また、係数Ａｉは節点座標のみで決まる定数である。有限要素法の性質より次式が成り立つ。
【数１１８】
【０１６２】
［三次元六面体・五面体における節点領域の算出方法１］
式（７４）を利用すると、式（７３）を次式に書き換えることができる。
【数１１９】
【０１６３】
ここで、Ｉ１ｉは節点１と節点ｉ間の距離であり、次式により算出する。
【数１２０】
【０１６４】
式(７５)の右辺の物理的な意味は、節点１と節点ｉの中点（Ｉ１ｉと記す）での節点１から節点ｉに向かう熱流束の２次近似λ（Ｔ１−Ｔｉ）／Ｉ１ｉにより、熱流束に垂直で中点Ｉ１ｉを通る面積（−Ａｉ）Ｉ１ｉを通過する熱量である。この流出熱量により節点１に割り当てられた熱源の体積を、図４９に示す節点１を頂点とし（−Ａｉ）Ｉ１ｉを底面積とした錐体とする。この錐体の体積をＮＤ１ｉで表わすと、次式が得られる。
【数１２１】
【０１６５】
各熱流束に割り当てる熱源の体積を合計したものを節点１の節点領域ＮＤ１とすると、次式が得られる。
【数１２２】
節点領域の求め方は、保存則によりある物理量の支配方程式を導出するときに使用されるコントロール・ボリュームの考え方と合致する。式（７８）の節点領域を式（７２）右辺に代入して、精度向上を図る。なお、係数Ａｉ（ｉ＝２，３，・・・，Ｎ）の中、正のものが存在するケースもあるが、この場合は、熱が錐体に流入することになっているためＮＤ１ｉが負になる。
【０１６６】
［三次元六面体・五面体における節点領域の算出方法２］
Ａｉ（ｉ＝１，２，・・・，Ｎ）中のＡｉ＞０のものをそれぞれＢｊ（ｊ＝１，２，・・・，Ｎ＋）で表わし、Ｎ＋でその数を表わす。また、Ａｉ≦０のものをそれぞれＣｋ（ｋ＝１，２，・・・，Ｎ−）で表わし、Ｎ−でその数を表わすと、式（７３）を次式で表わすことができる。
【数１２３】
【０１６７】
ここで、添え字ｉは要素に属する全節点を有限要素法で決めたルールに従って１から数えるものであるが、ｊは係数が正の節点のみを任意の順序で１から数え、ｋは係数が負の節点のみを任意の順序で１から数えるものである。従って、例えｉ，ｊ，ｋが同じ数字をとったとしても、ＴｉとＴｊ、ＴｉとＴｋとは必ずしも同一節点での温度ではなく、また、ＴｊとＴｋとはいつも異なる節点での温度である。式（７４−２）より次式が成り立つことが分かる。
【数１２４】
【０１６８】
各係数の値を棒の長さで表わすと、式（８０）の関係式を図５０で表現することができる。各係数間の繋ぎ目で区分すると、式（７３）を次式のように係数が正の節点温度と係数が負の節点温度との差の和に書き換えることができる。
【数１２５】
【０１６９】
ここで、Ｄｊｋは正の値のものでＢｊと−Ｃｋとが重なる部分である。Σはｊｋの組み合わせの合計を意味し、ｌｊｋは節点ｊとｋ間の距離である。式(８１)の右辺の物理的な意味は、節点ｊとｋの中点（Ｉｊｋと記す）での節点ｊから節点ｋに向かう熱流束の２次近似λ（Ｔｊ−Ｔｋ）／Ｉｊｋにより、熱流束に垂直で中点Ｉｊｋを通る面積Ｄｊｋｌｊｋを通過した熱量である。この熱量により節点１に割り当てられた熱源の体積は、図５１に示す節点１を頂点としＤｊｋｌｊｋを底面積とした錐体とする。この体積をＮＤｊｋで表わすと、次式が得られる。
【数１２６】
【０１７０】
式（８２）の右辺の［］内のものは、節点１から中点Ｉｊｋまでのベクトルと節点ｊから節点ｋまでのベクトルとの内積である。両ベクトルのなす角度が−π／２〜π／２の範囲内にある場合は、熱量が錐体から流出するためＮＤｊｋが正になる。その角度がこの範囲外にある場合は、熱量が錐体に流入するためＮＤｊｋが負になる。各熱流束に割り当てる熱源の体積を合計したものを節点１の節点領域ＮＤ１とすると、次式が得られる。
【数１２７】
このように導出した節点領域を式（７２）の右辺に代入して、精度向上を図る。
【０１７１】
［整合性、唯一性、正規性について］
整合性とは、節点領域とコントロール・ボリュームとが同体積であることを意味する。整合性が成立すれば、節点レベルでの保存則が２次精度で満足されることになり、精度改善効果が理論上で保証される。
【０１７２】
節点領域の考え方を図１６の長方体要素に適用した場合の整合性について検討するため、式（２５）を次式に書き換える。
【数１２８】
【０１７３】
ここで、α＝ｈｙｈｚ／１８ｈｘ、β＝ｈｘｈｚ／１８ｈｙ、γ＝ｈｘｈｙ／１８ｈｚである。上述の節点領域の算出方法１により次式を得ることができる。
【数１２９】
【０１７４】
同様にＮＤ１ｉ（ｉ＝３，４，・・・，８）を求めて、合計すると、次式が得られ、式（２６）と比較するとコントロール・ボリュームと一致することが分かる。
【数１３０】
【０１７５】
ほかの節点についても、要素形状の対称性を考慮すると整合性が満たされることが分かる。同様に図１８の正三角柱要素においても、次式が成立することが証明できる。式（３０）と比較すると整合性が満たされることが分かる。
【数１３１】
【０１７６】
また、正規性とは、
【数１３２】
すなわち、式（７２）の右辺への追加項の同要素全節点についての合計がゼロになることを意味する。正規性が成立すれば、要素レベルでも、全解析領域でも保存則が満たされることが保証される。長方体要素と正三角柱要素において、正規性が保たれていることが、式（８６）と式（８７）の結果、および要素形状の対称性から分かる。
【０１７７】
また、節点領域の唯一性とは、上記算出方法１と算出方法２で同じ結果が得られることを意味する。その理論証明は、算出方法２におけるバリエーションが多いため煩雑である。
【０１７８】
一般形状の六面体要素、五面体要素においては、ヤコビアンが絡んでくるため、整合性、正規性と唯一性について理論上の議論は困難である。代わりに、整合性については、後述の「作用及び効果」の欄で述べる数値解析において精度改善について検証する。正規性と唯一性についての検討は、以下の形状を持つ要素において行う。立方体領域（−５≦ｘ，ｙ，ｚ≦５）を幅２の立方体のメッシュに分割し、節点１をｘ＝ｙ＝ｚ＝−１に、節点２をｘ＝１、ｙ＝ｚ＝−１に置き、残る六節点（六面体要素の場合）、あるいは四節点（五面体要素の場合）はメッシュにおける全格子点に置いた組み合わせでできた全ての要素形状について検証を行った。これらの要素形状ではアスベクト比が３以下の実用解析で使用される要素形状をカバーできたといえる。この際、節点領域の算出方法２においては、正の係数の並ぶ順序については節点番号１から探し出して、負の係数の並ぶ順序については節点番号１からのと節点番号Ｎからの二通りで行った。その結果、体積がゼロ以下の要素と凹形状の要素を除いて、全テスト要素において正規性と唯一性が満たされていることが確認できた。
【０１７９】
［三次元六面体・五面体における作用及び効果］
保存型離散化手法の精度改善効果を二つの熱伝導問題の解析事例により示す。この際、解析領域を（０≦ｘ，ｙ，ｚ≦１）の立方体領域とし、λ＝１とした。基本境界条件問題としては、式（８８）のソース項分布に式（８９）の境界条件を加えた。その厳密解は式（９０）となる。また、基本・自然境界条件問題としては、式（９１）のソース項分布に式（９２）の基本境界条件と式（９３）の自然境界条件を加えた。その厳密解は式（９４）となる。
【数１３３】
【０１８０】
解析の数値誤差については次式により評価する。
【数１３４】
ここで、ＴＡｉ：節点ｉでの厳密解、ＴＮｉ：節点ｉでの数値解、ＴＡｍａｘ：厳密解の最大値、である。
【０１８１】
数値解析に使用する八パターンの六面体メッシュを図５２に示しており、ｘｙ平面上の四辺形メッシュをｚ方向に押し出して１６層を作成したものである。“●”で節点を、括弧内の数字で総節点数を表わす。メッシュが対称性を持つため、図中にはｘｙ平面上の４分の１の部分のみを示す。図５３（ａ）に基本境界条件問題、図５３（ｂ）に基本・自然境界条件問題での解析誤差を示す。長方体以外のメッシュでは、従来の有限要素法（ＧＦＥ）に比べ、ソース項の保存型離散化（ＣＳＦＥ）の解析精度が大幅に改善されたことが確認できる。また、解析誤差がメッシュ幅の自乗、すなわち、３次元問題では１／Ｎｏｄｅ２／３と比例すると指摘された（非特許文献１）ため、図５４には直線で示す長方体メッシュでの結果を基準にして、各要素形状における節点数と解析誤差との関係を示す。ソース項の保存型離散化（ＣＳＦＥ）は、長方体に近いまで解析精度が改善されたことが確認できる。
【０１８２】
図５５に使用する四パターンの五面体メッシュを示しており、ｘｙ平面上の三角形メッシュをｚ方向に押し出して１６層を作成したものである。メッシュの対称性を考慮して図中にはｘｙ平面上の４分の１のメッシュのみを示す。図５６（ａ）に基本境界条件問題、図５６（ｂ）に基本・自然境界条件問題での解析誤差を示す。特に、自然境界を含む解析条件では、従来の有限要素法（ＧＦＥ）に比べ、ソース項の保存型離散化（ＣＳＦＥ）の解析精度が大幅に改善されたことが確認できる。
【０１８３】
本提案では、解析精度の改善や、メッシュ作成工数の低減、解析時間の短縮などの効果が得られる。市販の解析ソフトに簡単な修正を加えるだけで実現でき、計算コストも増えない。
【０１８４】
［本発明の簡単に応用展開できる解析分野］
熱伝導問題を例に本発明の応用について説明してきたが、勿論、Poissonの方程式が記述する現象の全てにそのまま展開できる。例えば、流れ問題、電・磁気問題、弾性変形問題、物質拡散問題などが応用できる分野として挙げられる。
【０１８５】
本発明の基本的な考え方は、二次精度の保存則の観点から楕円型演算子の離散化結果におけるコントロール・ボリュームを見出し、ソース項の積分領域をそれに合わせたものである。したがって、この考え方は、保存則から導出された楕円型演算子を含む各種の支配方程式の有限要素解析に容易に応用できる。すなわち、物理量の発生・消失と流出・流入を算出するとき、その積分領域について、本発明で導出した節点領域の考え方を導入して積分すればよい。例えば、非定常移流拡散解析において、非定常項と移流項の離散化に対して本発明におけるソース項の積分領域の考え方が適用できる。
【０１８６】
また、荷重、体積力、質量が係わる従来の有限要素法による解析においては、それらの扱い方は、式（２３）の右辺におけるソースと同じように、
【数１３５】
により節点方程式に分配されるため、不整合性問題も存在する。これらの問題は本発明でのソース項の扱い方で改善することもできる。三次元問題の解析においては、荷重、体積力、質量について式（７２）の右辺のソース項と同様に追加項を加えて修正し、対象節点Ｐに寄与する荷重、体積力、質量は次式により算出する。
【数１３６】
ここで、ｆは荷重の任意成分、体積力の任意成分、質量を意味する。
【０１８７】
［本発明のまとめ］
以上説明したように、本発明に係る有限要素法を用いた解析方法にあっては、一般関数項（熱問題等の場合はソース項、荷重、体積力、質量が係る問題の場合はフォース項）を積分する一般関数項積分工程（ソース項積分工程、フォース項積分項手値）（ステップＳ４）にあって、Galerkin有限要素法による二階微分項の離散化結果を基に定義された節点領域の考え方を導入し、要素の代表値を用いた一般関数項（ソース項、フォース項）を積分するので、節点領域の大きさに見合った値（ソース量）を取り入れることができる。これにより、メッシュ・パターンによる解析誤差を低減することができ、品質が低い要素（品質の低いメッシュ・パターン）からも高精度な数値解を得ることが可能となる。このため、メッシュ・パターンの品質を高める作業（メッシュ・パターンを例えば正三角形、正四面体、正方形に近づける作業）を不要とすることができるものでありながら、アルゴリズム（プログラム）の修正による計算量の増加が殆んどないので、コンピュータによる演算処理の増大も防止することができ、総じて解析時間の短縮化を図ることができる。
【０１８８】
また、境界条件については、節点領域の範囲と合致したものを導入する境界条件導入工程（ステップＳ６）を設けたので、節点レベルで精度を向上させることができる。これにより、さらにメッシュ・パターンによる解析誤差を低減することができ、品質が低い要素（品質の低いメッシュ・パターン）からも高精度な数値解を得ることが可能となる。
【０１８９】
なお、以上説明した実施の形態においては、特に熱問題等において、ソース項に、二次元三角形の場合は（ＮＤ−Ｓ／３）ＱＧを追加項とし、三次元四面体の場合は（ＮＤ−Ｖ／４）ＱＧを追加項とし、二次元四角形の場合は
【数１３７】
を追加項とし、三次元六面体・五面体の場合は、
【数１３８】
を追加項としてソース項に追加したものを説明したが、これに限らず、ソース項は要素の幾何学中心でのソース値を用いたソース項に変更する手法であれば、例えばソース項（式（３３）、式（３６）、式（４６）、式（７２）の右辺全部）を単純に「ＮＤ・ＱＧ」に置き換えることも考えられる。
【０１９０】
また、本実施の形態においては、特に荷重、体積力、質量問題等において、一般関数項に、二次元三角形の場合は（ＮＤ−Ｓ／３）ｆＧを追加項とし、三次元四面体の場合は（ＮＤ−Ｖ／４）ｆＧを追加項とし、二次元四角形の場合は
【数１３９】
を追加項とし、三次元六面体・五面体の場合は
【数１４０】
を追加項として一般関数項に追加したものを説明したが、これに限らず、一般関数項は要素の幾何学中心での値を用いた一般関数項に変更する手法であれば、例えば一般関数項を単純に「ＮＤ・ｆＧ」に置き換えることも考えられる。
【０１９１】
また、本実施の形態では、有限要素法を用いた解析方法、及び有限要素法を用いた解析プログラム、として本発明を説明したが、勿論、これら解析方法や解析プログラムを実行する解析装置（解析システム）や、当該プログラムを記録した記録媒体などの形態で本発明を利用するものも、本発明の適用範囲内である。
【産業上の利用可能性】
【０１９２】
本発明に係る有限要素法を用いた解析方法、及び有限要素法を用いた解析プログラムは、解析対象の熱伝導解析、流体の流れ解析、電気流れ解析、磁気流れ解析、弾性変形解析、物質拡散解析、荷重解析、体積力解析、質量解析、等の各種解析に用いることが可能であり、特にメッシュ・パターンの品質が低いものであっても誤差の低減が求められるような有限要素法の解析に用いて好適である。
【符号の説明】
【０１９３】
１ 要素における節点
２ 要素における節点
３ 要素における節点
４ 要素のおける節点
Ｐ 解析領域における節点
ＮＤ 節点領域
Ｑ 一般関数、ソース項
ＱＧ 要素ソース代表値
ＱＯ 要素ソース代表値
Ｗ Galerkin重み関数
ｅ 要素
ｆ 一般関数、荷重、体積力、質量項
ｆＧ 要素荷重、体積力、質量の代表値
ｆＯ 要素荷重、体積力、質量の代表値
Ω １節点周りの要素領域
Ｓ１ 選定工程（ステップＳ１）
Ｓ２ 分割工程（ステップＳ２）
Ｓ３ 要素マトリクス作成工程（ステップＳ３）
Ｓ４ 一般関数項積分工程、ソース項積分工程（ステップＳ４）
Ｓ５ 連立方程式作成工程（ステップＳ５）
Ｓ６ 境界条件導入工程（ステップＳ６）
Ｓ７ 演算工程（ステップＳ７）
Ｓ８ 演算工程（ステップＳ８）
【技術分野】
【０００１】
本発明は、有限要素法を用いた解析方法、及び有限要素法を用いた解析演算プログラムに係り、特に要素の形状（メッシュ・パターン）によって生じる誤差の低減を図った有限要素法を用いた解析方法、及び有限要素法を用いた解析演算プログラムに関する。
【背景技術】
【０００２】
一般的に有限要素法では、形状関数と同形なGalerkin重み関数を微分方程式の各項にかけ、節点周りの要素領域にわたって積分することにより数値解を求めている。有限要素法は、誤差評価が明瞭で、汎用性に富むなどの特徴があって、最も広く使用される数値解法の一つになっている。しかし、メッシュ・パターンによって解析結果が異なったり、数値誤差が大きくなる問題が依然として解決されていない。
【０００３】
［従来の有限要素法による二次元三角形の問題］
例えば、線形三角形メッシュ上でポアソン（Poisson）の方程式を離散化した場合には、もともと同じ解が出るべきところから、要素形状によって異なる解析結果が出されてしまう（非特許文献１参照）。また、非特許文献２に示す２パターンのメッシュから得た代数方程式を比較すると、二階微分項の離散化結果は同形になっているものの、ソース項の離散化結果が異なっている。ソース項に重み関数をかけて積分した結果、ソース量が要素形状と関係なくほぼ均等的に各節点の代数方程式に分配され、ほかの項との整合性が取れず、節点レベルでは保存法則が満たされない状況になっている。この不整合性の詳細は、以下のとおりである。
【０００４】
二次元ポアソンの方程式は下記のものである。
【数１】
【０００５】
この方程式は、多分野の物理現象を記述するものであるが、ここでは、熱伝導問題を例に説明を展開し、Ｔは温度で、λは熱伝導率で、Ｑはソース項と呼ばれる単位時間、単位面積あたりの発熱量である。式（１）に有限要素法を適用すると、考えている対象節点Ｐについての積分方程式として次式が得られる。
【数２】
【０００６】
ここで、Ωは図１に示す節点Ｐ周りの要素からなす領域であり、ＷはΩ上で定義された節点Ｐで１を、Ωの境界で０をとるGalerkin重み関数である。Ωにわたる積分は、各要素上の積分の和になるので、式（２）を式（３）に書き換えることができる。
【数３】
【０００７】
ここで、ｅは節点Ｐ周りの要素を表わし、
【数４】
は各要素上の演算の和を意味する。要素ｅを図２の線形三角形要素に当てはめた場合には、式（３）左辺の要素積分の結果を次式のように整理することができる。
【数５】
ここでは、局所座標系を使用しており、対象節点Ｐを局所節点１とし、数字の添え字で節点番号を表わす。
【０００８】
式（４）右辺の１項目中の
【数６】
は図３に示す中点Ｉでの節点１から節点２に向かう熱流束の２次精度の近似であり、その係数ｄＣＩは、点Ｉから外心Ｃまでの距離と等しくなっている。すなわち、式（４）右辺の１項目は点Ｉでの２次精度の熱流束を用いて、垂直な線分ＩＣを通過する熱量を評価するものである。同様に、２項目は、点Ｋでの２次精度の熱流束を用いて、垂直な線分ＣＫを通過する熱量を評価するものである。以上の考察をまとめると、式（４）の右辺は、ＩＣとＣＫという検査線を通して流出した熱量であることが分かる。この際、四辺形１ＩＣＫを節点１の節点領域と呼び、ＮＤ１で表わす。正三角形要素では外心が幾何中心と重なるため、節点領域は要素面積の３分の１になる。直角三角形要素では外心が辺中点に重なり、鈍角三角形要素では外心が要素の外に出る。これらの要素における節点１の節点領域は図４と図５に示すようになる。
【０００９】
一方、式（３）右辺のソース項の要素積分により、要素内の発熱量が要素形状と関係なくほぼ均等的に３節点に分配される。特に、Ｑが要素内で定数の場合には、ソース項の要素積分は次式になる。
【数７】
ここで、Ｓは要素面積である。すなわち、いかなる要素形状においても、節点に寄与された熱量は要素の３分の１の面積から発生したものだけである。
【００１０】
ここで、正三角形以外の要素では、検査線で決めた節点領域の面積とソース項の面積との間は不釣り合いが発生する。同じことは、節点Ｐ周りのほかの要素からも生じるので、保存法則が満たされず、ソース項と他の項の間に不整合性が表れた。その結果、新たな数値誤差が引き起こされる。図６には、有限要素法から得られた温度分布の一例を示す。０≦ｘ≦１、０≦ｙ≦１の正方形領域において、式（６）のソース分布と式（７）の境界条件を用いて、λ＝１とした解析結果を実線で示し、点線で要素形状を示す。例えば、○印を付けた中央部の４点では同値の厳密解は得られるが、局所メッシュ・パターンの影響で数値解には差が出ている。
【数８】
【数９】
【００１１】
［従来の二次元三角形問題の改善方法］
前述の不整合性問題を解決するため、以下の２つの改善策が考えられる。
（改善方法１）：式（３）右辺の要素積分
【数１０】
を
【数１１】
に置き換える。
（改善方法２）：上記改善方法１をさらに修正したものであり、鈍角三角形要素における要素外での積分を避けるため、ソース項の積分を図７に示す要素内の部分で行う。それに合わせて式（４）右辺のｄＣＩとｄＣＫを中点から外心までの長さから図７に示す長さに修正する。
【００１２】
上記改善方法１では、全ての線形三角形要素において保存法則が２次精度で満たされ、各項間の整合性が維持できるが、鈍角三角形要素の場合には、ソース項の積分範囲が要素外までに広がることが、数値振動などを引き起こす原因になっている。また、上記改善方法２では、左辺のマトリックスが強制的に変えられたため、楕円型演算子の最適な離散化方法として評価されてきた有限要素法の優位性が損なわれる。結局のところ、上記改善方法１，２では満足できる数値解が得られず、実用上では鈍角三角形要素の使用を避けており、言い換えると、鈍角三角形要素を用いないメッシュ・パターンを時間をかけて作成する必要がある。
【００１３】
［従来の有限要素法による三次元四面体の問題］
三次元のポアソンの方程式は次式になる。
【数１２】
【００１４】
有限要素法を適用すると対象節点についての積分式として次式が得られる。
【数１３】
【００１５】
式（９）は式（１０）に示す節点周りの要素積分の和で表現することができる。
【数１４】
【００１６】
式（１０）を、例えば、図８のような三つの面が節点１で垂直に交わる線形四面体要素に適用すると、左辺の要素積分は次式になる。
【数１５】
【００１７】
式（１１）の右辺から、節点１から節点２、３、４に向かう熱流束
【数１６】
の通過面積はそれぞれ
【数１７】
になることが分かる。同様な方法で、ほかの節点間に伝わる熱流束の通過面積も求められる。それらを表１にまとめて示す。
【００１８】
【表１】
【００１９】
一方、式（１０）右辺のソース項の要素積分を行うと、要素内の発熱量が要素形状と関係なく、ほぼ均等的に４節点に分配される。特に、Ｑが要素内で定数の場合には、ソース項の要素積分は次式になる。
【数１８】
【００２０】
ここでのＶは要素体積である。すなわち、いかなる要素形状においても、節点に寄与された熱量は要素の４分の１の体積から発生したものだけである。ここで、保存法則の観点から表１と式（１２）を考察すると、節点によって大きく変る熱流束の通過面積と均等的に分配されるソース項の間に不整合性が表れたことが分かる。その結果、新たな数値誤差が引き起こされる。図９には、０≦ｘ≦１、０≦ｙ≦１、０≦ｚ≦１の立方体領域において、λ＝１とし、式（１３）のソース分布と式（１４）の境界条件に有限要素法を適用して得られた温度分布を示す。例えば、図９（ｂ）の○印を付けた４点では同値の厳密解は得られるが、局所メッシュ・パターンの影響で数値解は異なっている。
【数１９】
【数２０】
【００２１】
［従来の三次元四面体問題の改善方法］
この不整合性の三次元四面体の問題の改善方法としては、図１０のように、対象節点１、四面体の外心Ｃ、要素表面三角形の外心Ｃ２、Ｃ３、Ｃ４と要素辺の中点Ｉ、Ｊ、Ｋの合計８点から成形した六面体を節点領域として定義し、それに合わせて、マトリックス成分を節点領域から算出した流束の通過面積を節点間の距離で割ったものに修正する方法が考えられる。マトリックス成分を変える理由としては、三次元有限要素法の離散化結果には二次元と同様な解釈が適用できないからである。すなわち、式（１１）左辺を積分した結果は、正四面体要素の場合を除いて、辺中点での熱流束と、外心および辺中点からなす四辺形の面積との掛け算にならない。例えば、図８の要素の場合には、節点１と節点２間の熱流束の通過面積は
【数２１】
になっているが、外心Ｃ、Ｃ３、Ｃ４と中点Ｉとによりできた四辺形の面積は
【数２２】
になる。
【００２２】
この三次元問題の改善方法は、保存法則を２次精度で満足するが、デロネー（Delaunay）分割の要素形状にしか適用できないという欠点が挙げられる。また、たとえば、図８のような四面体要素に適用すると、従来の有限要素法では対角線以外の成分がゼロになることで済むが、この改善方法では、マトリックス成分を変えたことにより、対角線以外のところで数値解析に悪い影響を及ぼす正成分が発生する。
【００２３】
［従来の有限要素法による二次元四角形の問題］
上述した二次元の線形三角形メッシュ及び三次元の線形四面体メッシュについては、上述したような改善方法も提案されているが（非特許文献３、非特許文献４参照）、二次元の線形四辺形メッシュについては議論されていないようである。有限要素法によりPoissonの方程式を離散化した場合には、全ての項に同じGalerkin重み関数をかけて積分するため、ソース項が要素形状に関係なく同要素が所有する各節点の代数方程式にほぼ均等的に分配されることがあり、ほかの項の離散化との整合性が取れず、節点レベルでは保存則が満たされない状況になっている。その詳細は、以下のとおりである。
【００２４】
上述したように、二次元Poissonの方程式は下記のものである。
【数２３】
【００２５】
この方程式は、多分野の物理現象を記述するものであるが、ここでは、熱伝導問題を例に説明を展開し、Ｔは温度で、λは熱伝導率で、Ｑはソース項と呼ばれる単位時間、単位面積あたりの発熱量である。式(1)に有限要素法を適用すると、考えている対象節点Ｐについての積分方程式は、上述した式（３）と同じ式（３）’が得られる。
【数２４】
【００２６】
ここで、ｅは図１１に示す節点Ｐ周りの要素を表わし、Σは各要素上の演算の和を意味する。また、ＷＰは節点Ｐ周りの要素領域上で定義されたGalerkin重み関数であり、節点Ｐでは１を、図１１の局所領域の境界ではゼロをとる線形関数である。部分積分とガウス・グリーン定理を利用して式（３）’左辺の積分を行い、境界積分を別途考えると、次式が得られる。
【数２５】
【００２７】
要素ｅを図１２の双一次長方形要素に当てはめた場合には、式（１５）左辺の要素積分は次式となる。ただし、対象節点Ｐの要素における局所番号を１とする。
【数２６】
【００２８】
式（１６）右辺１項目の物理的な意味は、中点Ｉでの節点１から節点２に向かう２次精度の熱流束と、中点Ｋでの節点４から節点３に向かう２次精度の熱流束との重み付き平均により、垂直線分ＩＯを通過する熱量である。図１２中の点Ｏは、四節点の座標の平均値に位置する。また、２項目の物理的な意味は、中点Ｌでの節点１から節点４に向かう２次精度の熱流束と、中点Ｊでの節点２から節点３に向かう２次精度の熱流束との重み付き平均により、垂直線分ＯＬを通過する熱量である。ここで、検査線ＩＯＬと要素辺Ｌ１Ｉで囲う領域を節点１のコントロール・ボリュームと称しＣＶ１と記す。要素形状の対称性を考えると、ほかの節点を対象節点とした検査線の検討結果が、図１２の細い実線で示すようになり、ＣＶ１＝ＣＶ２＝ＣＶ３＝ＣＶ４＝Ｓ／４となることが分かる。
【００２９】
ここで、Ｓは要素面積であり、ＣＶの添え字で所属節点を表わす。一方、Ｑが要素内で定数の場合には、式（１５）右辺の要素積分は次式になる。
【数２７】
すなわち、節点１に寄与された熱量はＳ／４から発生したものである。ほかの節点への寄与量も同じであることが重み関数の定義から分かる。各節点のコントロール・ボリュームの面積とソース項の面積とが一致するため、保存則は２次精度で満たされている。これらの結果を表２にまとめて示す。
【００３０】
【表２】
【００３１】
なお、図１３に示すような双一次平行四辺形要素の場合には、有限要素法による式（１５）左辺の要素積分は次式になる。
【数２８】
【００３２】
式（１８）右辺の１項目は、中点Ｉでの節点１から節点２に向かう２次精度の熱流束により、垂直線分ＩＪを通過する熱量である。２項目は点Ｏでの節点１から節点３に向かう２次精度の熱流束により垂直線分ＪＬを通過する熱量である。結局、検査線はＩＪＬとなり、ＣＶ１は３Ｓ／８になる。同様に、節点２を対象節点として検討すると、式（１５）左辺の要素積分は次式になる。
【数２９】
検査線はＩＪとなり、ＣＶ２はＳ／８になる。残る二節点については、要素形状の対称性を考えると、図１３の細い実線で示す検査線を得ることができる。
【００３３】
一方、Ｑが要素内で定数の場合には、ソース項の要素積分結果はＳ／４から発生した熱量になる。上記表２にこれらの結果をまとめて示す。節点に寄与された熱量の発生面積とコントロール・ボリュームの面積とが一致しないため、保存則は満たされていない。
【００３４】
図１４に示す節点３が節点４に限りなく接近する四辺形要素の場合には、有限要素法による式（１５）左辺の要素積分は次式になる。ただし、Ｔ３＝Ｔ４と扱うことができる。
【数３０】
【００３５】
式（２０）右辺の１項目は、中点Ｉでの節点１から節点２に向かう２次精度の熱流束により、垂直線分ＩＪを通過する熱量である。２項目は点Ｌでの節点１から節点４に向かう２次精度の熱流束により垂直線分ＪＬを通過する熱量である。結局、検査線はＩＪＬとなり、ＣＶ１はＳ／２になる。同様に、節点２を対象節点として検討すると、式（１５）左辺の要素積分は次式になる。
【数３１】
【００３６】
その検査線はＩＪとなり、ＣＶ２はＳ／４になる。また、Ｔ３＝Ｔ４を考慮して節点３と４を対象節点として積分した結果をそれぞれ式（２２）と式（２３）に示す。
【数３２】
【数３３】
【００３７】
両節点が重なり、Ｔ３＝Ｔ４であるため、検査線について検討するとき、代数方程式を足し合わせて扱うことができる。式（２２）と式（２３）との右辺の合計から検査線としてＪＬが得られ、ＣＶ３＋ＣＶ４＝Ｓ／４になることが分かる。図１４には細い実線で検査線を示す。
【００３８】
一方、Ｑが要素内で定数の場合には、節点１と節点２に寄与されるソース量の面積がそれぞれＳ／３になる。節点３と節点４に寄与されるソース量の面積の和がＳ／３になる。これらの結果を上記表１にまとめて示すが、保存則は満たされていないことが分かる。
【００３９】
なお、このような２次精度で保存則が満たされない問題は、ほかの一般形状の四辺形要素にも存在すると考えられるが、理論上では、コントロール・ボリュームについて調べることが困難である。
【００４０】
［従来の有限要素法による三次元六面体及び三次元五面体の問題］
上述した二次元の線形三角形メッシュ及び三次元の線形四面体メッシュについては、上述したような改善方法も提案されているが（非特許文献３、非特許文献４参照）、三次元解析によく使用される六面体要素と五面体要素については議論されていないようである。本発明では、これらの要素を検討対象とする。
【００４１】
上述したように三次元Poissonの方程式は上記式（８）のものである。
【数３４】
【００４２】
この方程式は、多分野の物理現象を記述するものであるが、ここでは、熱伝導問題を例に説明を展開し、Tは温度で、λは熱伝導率で、Ｑはソース項と呼ばれる単位時間、単位面積あたりの発熱量である。式（８）に有限要素法を適用すると、考えている対象節点Ｐについての積分方程式として次式が得られる。
【数３５】
【００４３】
ここで、ｅは図１５に示す節点Ｐ周りの要素を表わし、Σは各要素上の演算の和を意味する。また、ＷＰは節点Ｐ周りの要素領域上で定義されたGalerkin重み関数で、節点Ｐでは１を、周りの要素から構成した領域の境界ではゼロをとる線形関数である。部分積分とガウス・グリーン定理を利用して式（１０）’左辺の積分を行い、境界積分を別途考えると、次式が得られる。
【数３６】
【００４４】
［六面体要素の場合］
要素ｅを図１６（ａ）の線形長方体要素に当てはめた場合には、式（２４）左辺の要素積分結果を次式に整理することができる。ただし、対象節点Ｐの要素における局所番号を１とする。
【数３７】
【００４５】
図１６(ｂ)に示すように、点Ｉｉｊで節点ｉとｊを結ぶ線分の中点を表わしておく。式(２５)右辺１項目の物理的な意味は、中点Ｉ１２での節点１から節点２に向かう２次精度の熱流束λ（Ｔ１−Ｔ２）／２ｈｘ、中点Ｉ４３での節点４から節点３に向かう２次精度の熱流束λ（Ｔ４−Ｔ３）／２ｈｘ、中点Ｉ５６での節点５から節点６に向かう２次精度の熱流束λ（Ｔ５−Ｔ６）／２ｈｘ、と中点Ｉ８７での節点８から節点７に向かう２次精度の熱流束λ（Ｔ８−Ｔ７）／２ｈｘ、との重み付き平均により、熱流束に垂直で中点Ｉ１２を通る断面積ｈｙｈｚを通過する熱量である。また、２項目の物理的な意味は、中点Ｉ１４での節点１から節点４に向かう２次精度の熱流束、中点Ｉ２３での節点２から節点３に向かう２次精度の熱流束、中点Ｉ５８での節点５から節点８に向かう２次精度の熱流束、と中点Ｉ６７での節点６から節点７に向かう２次精度の熱流束、との重み付き平均により、熱流束に垂直で中点Ｉ１４を通る断面積ｈｘｈｚを通過する熱量である。また、３項目の物理的な意味は、中点Ｉ１５での節点１から節点５に向かう２次精度の熱流束、中点Ｉ２６での節点２から節点６に向かう２次精度の熱流束、中点Ｉ４８での節点４から節点８に向かう２次精度の熱流束、と中点Ｉ３７での節点３から節点７に向かう２次精度の熱流束、との重み付き平均により、熱流束に垂直で中点Ｉ１５を通る断面積ｈｘｈｙを通過する熱量である。ここで、断面積ｈｙｈｚ、ｈｘｈｚとｈｘｈｙを検査面として考え、これらの検査面が囲う領域を節点１のコントロール・ボリュームと称しＣＶ１と記す。このコントロール・ボリュームは、２次精度の熱流束により定義されたものである。以上の議論から、ＣＶ１は図１６(ｂ)に示す辺長ｈｘ、ｈｙ、ｈｚの長方体となり、その体積が次式となることが分かる。
【数３８】
【００４６】
ここで、Ｖは要素体積である。要素形状の対称性を考えると、ＣＶ１＝Ｖ／８(ｉ＝１，２，・・・，８)の結果が得られる。ＣＶの添え字で所属節点を表わす。一方、Ｑが要素内で定数の場合には、式(２４)右辺の要素積分は次式になる。
【数３９】
【００４７】
ここで、Ｑ(ｅ)は要素でのソース項である。式（２７）より節点１に寄与された熱量はＶ／８から発生したものであることが分かる。ほかの各節点への寄与量もＶＱ(ｅ)／８になることがGalerkin重み関数の定義から分かる。各節点のコントロール・ボリュームの体積とソース項の体積とが一致するため、保存則は２次精度で満たされている。
【００４８】
なお、図１７（ａ）に示すｙｚ平面上の平行四辺形をｘ方向に２ｈｘ分押し出して生成した六面体要素の場合には、有限要素法による式（２４）左辺の要素積分は次式になる。
【数４０】
【００４９】
式(２８)右辺の物理的な意味は、それぞれ広さがｈ２で中点Ｉ１２を通る面積、広さがｈｘｈで中点Ｉ１４を通る面積、広さが２ｈｈｘで中点Ｉ１８を通る面積、を通過する熱量である。コントロール・ボリュームは、検査面が閉じていないため厳密に求めることができないが、図１６(ｂ)に比べると、Ｖ／８以上になることが考えられる。一方、Ｑが要素内で定数の場合には、各節点のソース項の要素積分結果は依然としてＶ／８から発生した熱量になる。すなわち、節点に寄与された熱量の発生体積とコントロール・ボリュームの体積とが一致せず、保存則が満たされていないことが考えられる。
【００５０】
このような２次精度で保存則が満たされない問題は、ほかの一般形状の六面体要素にも存在すると考えられる。しかし、コントロール・ボリュームについて理論上で調べることが困難である。
【００５１】
［五面体要素の場合］
要素ｅを図１８（ａ）の線形正三角柱要素に当てはめた場合には、式（２４）左辺の要素積分は式（２９）となる。ただし、対象節点Ｐの要素における局所番号を１とする。
【数４１】
【００５２】
同様に、式（２９）右辺の物理的な意味について検討すると図１８（ｂ）に示す検査面が得られるため、節点１のコントロール・ボリュームは次式となる。
【数４２】
【００５３】
要素形状の対称性を考えると、ＣＶｉ（ｉ＝１，２，・・・，６）の結果が得られる。一方、Ｑが要素内で定数の場合には、式(２４)右辺の要素積分は次式になる。
【数４３】
【００５４】
すなわち、節点１に寄与された熱量はＶ／６から発生したものである。ほかの節点への寄与量も同じであることがGalerkin重み関数の定義から分かる。各節点のコントロール・ボリュームの体積とソース項の体積とが一致するため、保存則は２次精度で満たされている。
【００５５】
なお、図１９(ａ)に示す直角三角柱要素の場合には、有限要素法による式(２４)左辺の要素積分は次式になる。
【数４４】
【００５６】
式(３２)右辺の物理的な意味は、それぞれ広さがｈｙｈｚ／２で中点Ｉ１２を通る面積、広さがｈｘｈｚ／２で中点Ｉ１３を通る面積、広さがｈｘｈｙ／６で中点Ｉ１４を通る面積、を通過する熱量である。コントロール・ボリュームは、検査面が閉じていないため厳密に求められないが、図１８(ｂ)に比べるとＶ／６以上になることが考えられる。一方、Ｑが要素内で定数の場合には、各節点のソース項の要素積分結果は依然としてＶ／６から発生した熱量になる。節点に寄与された熱量の発生体積とコントロール・ボリュームの体積とが一致せず、保存則が満たされていないことが考えられる。
【先行技術文献】
【非特許文献】
【００５７】
【非特許文献１】O.C.Zienkiewicz，R.L.Taylor著（矢川元基ら訳）、「マトリクス有限要素法」（特に第１０章参照）、科学技術出版社、１９９６年
【非特許文献２】邵 長城著、「基本からわかる有限要素法」（特に第５章参照）、森北出版、２００８年
【非特許文献３】MING-DER DONALD HUANG, The Constant-Flow Patch Test−A Unique Guideline for the Evaluation ofDiscretization Schemes for the Current Continuity Equations. IEEE Transactions on computer-aided design, V.CAD-4, No. 4, 1985, pp.583-608.
【非特許文献４】Christian Cordes and Mario Putti, Accuracy ofGalerkin finite element for groundwater flow simulations in two andthree-dimensional triangulations. Int.J. Numer. Meth. Engng 2001, V52, P371-387.
【発明の概要】
【発明が解決しようとする課題】
【００５８】
以上説明したような問題から、本発明では以下の（ａ）〜（ｅ）を解決しようとする課題とする。
（ａ）上記従来の二次元三角形の問題の改善方法１において、鈍角三角形要素の場合には、ソース項の積分領域が要素外までに広がる問題を解決する。同時に、境界条件と線上分布ソース項の扱い方を提案する。
（ｂ）従来の有限要素法による三次元四面体の問題におけるソース項との不整合性問題を解決し、上記従来の三次元四面体の問題の改善方法における弊害を回避する。同時に、境界条件、線上分布ソース項、面上分布ソース項の扱い方を提案する。
（ｃ）上記従来の二次元四角形の問題におけるソース項との不整合性問題を解決する。
（ｄ）上記従来の三次元六面体・五面体の問題におけるソース項との不整合性問題を解決する。
（ｅ）有限要素法における荷重、体積力、質量の扱い方においても要素形状によって誤差が発生する問題も存在するため、その精度改善方法を提案する。
【００５９】
そこで本発明は、線形三角形要素、線形四面体要素、線形四角形要素、線形六面体要素、線形五面体要素を用いた有限要素法の離散化において、高精度のスキームを提案して、メッシュ・パターンによる解析誤差を低減し、例えば正三角形、正四面体、長方形、長方体などよりも品質が低い要素（品質の低いメッシュ・パターン）からも高精度な数値解が得られる有限要素法を用いた解析方法、及び有限要素法を用いた解析演算プログラムを提供することを目的とするものである。
【課題を解決するための手段】
【００６０】
本発明は、解析対象の解析領域を選定する選定工程（Ｓ１）と、
前記解析領域を計算対象として複数の要素に分割する分割工程（Ｓ２）と、
前記複数の要素のうちの、ある要素（例えばｅ）にある節点（例えばＰ，要素ｅでの節点番号を例えば１とする）についてGalerkin重み関数（Ｗ）をかけて要素積分をし、各要素のマトリクスを作成する要素マトリクス作成工程（Ｓ３）と、
Galerkin重み関数（Ｗ）と一般関数（例えばＱ，ｆ）との積からなる一般関数項を積分する一般関数項積分工程（Ｓ４）と、
前記ある節点（例えばＰ）周りの領域（例えばΩ）の要素における各要素のマトリクスの和と、前記一般関数項を積分した値の和とに基づき、連立方程式を作成する連立方程式作成工程（Ｓ５）と、
前記連立方程式に境界条件の導入をする境界条件導入工程（Ｓ６）と、
前記連立方程式を演算して数値解を得る演算工程（Ｓ７，Ｓ８）と、を備えた有限要素法を用いた解析方法において、
前記一般関数項積分工程（Ｓ４）にあって、Galerkin有限要素法による二階微分項の離散化結果を基に定義された節点領域の考え方を導入し、要素の代表値（例えばＱＧ，ＱＯ，ｆＧ，ｆＯ）を用いた一般関数項を積分することを特徴とする。
【００６１】
また、本発明は、前記一般関数項積分工程（Ｓ４）にあって、節点領域の大きさに合わせて一般関数を節点に分配することを特徴とする。
【００６２】
更に、具体的に本発明は、前記一般関数項積分工程は、前記一般関数項をソース項としたソース項積分工程（Ｓ４）であり、
前記ソース項積分工程（Ｓ４）にあって、要素の代表ソース値（例えば、幾何学中心での値や節点平均位置での値、要素平均値など）を用いたソース項を積分することを特徴とする。
【００６３】
また、具体的に本発明は、前記一般関数項積分工程は、前記一般関数項を荷重、体積力、質量（フォースと総称する）としたフォース項積分工程（Ｓ４）であり、
前記フォース項積分工程（Ｓ４）にあって、要素の代表フォース値（例えば、幾何学中心での値や節点平均位置での値、要素平均値など）を用いたフォース項を積分することを特徴とする。
【００６４】
また、詳細には本発明は、前記複数の要素は、二次元三角形の要素であり、
前記ソース項積分工程（Ｓ４）にあって、前記各要素における節点領域の面積をＮＤ、前記要素の代表ソース値をＱＧ、前記要素の面積をＳとした場合に、（ＮＤ−Ｓ／３）ＱＧを追加項としてソース項に追加することを特徴とする。
【００６５】
また、詳細には本発明は、前記複数の要素は、二次元三角形の要素であり、
前記フォース項積分工程（Ｓ４）にあって、前記各要素における節点領域の面積をＮＤ、前記要素の代表フォース値をfＧ、前記要素の面積をＳとした場合に、（ＮＤ−Ｓ／３）fＧを追加項としてフォース項に追加することを特徴とする。
【００６６】
また、詳細には本発明は、前記複数の要素は、三次元四面体の要素であり、
前記ソース項積分工程にあって、前記各要素における節点領域の体積をＮＤ、前記要素の代表ソース値をＱＧ、前記要素の体積をＶとした場合に、（ＮＤ−Ｖ／４）ＱＧを追加項としてソース項に追加し、かつ前記ＮＤとして以下の数式を用いることを特徴とする。
【数４５】
【００６７】
また、詳細には本発明は、前記複数の要素は、三次元四面体の要素であり、
前記フォース項積分工程（Ｓ４）にあって、前記各要素における節点領域の体積をＮＤ、前記要素の代表フォース値をfＧ、前記要素の体積をＶとした場合に、（ＮＤ−Ｖ／４）ＱＧを追加項としてフォース項に追加し、かつ前記ＮＤとして以下の数式を用いることを特徴とする。
【数４６】
【００６８】
また、詳細には本発明は、前記複数の要素は、二次元四角形の要素であり、
前記ソース項積分工程（Ｓ４）にあって、前記各要素における節点領域をＮＤ、前記Galerkin重み関数をＷP、前記要素の代表ソース値をＱＯとした場合に、以下の数式を追加項としてソース項に追加することを特徴とする。
【数４７】
【００６９】
また、詳細には本発明は、前記複数の要素は、二次元四角形の要素であり、
前記フォース項積分工程（Ｓ４）にあって、前記各要素における節点領域をＮＤ、前記Galerkin重み関数をＷP、前記要素の代表フォース値をfＯとした場合に、以下の数式を追加項としてフォース項に追加することを特徴とする。
【数４８】
【００７０】
また、詳細には本発明は、前記複数の要素は、三次元六面体又は三次元五面体の要素であり、
前記ソース項積分工程（Ｓ４）にあって、前記各要素における節点領域をＮＤ、前記Galerkin重み関数をＷP、前記要素の代表ソース値をＱＯとした場合に、以下の数式を追加項としてソース項に追加することを特徴とする。
【数４９】
【００７１】
また、詳細には本発明は、前記複数の要素は、三次元六面体又は三次元五面体の要素であり、
前記フォース項積分工程（Ｓ４）にあって、前記各要素における節点領域をＮＤ、前記Galerkin重み関数をＷP、前記要素の代表フォース値をfＯとした場合に、以下の数式を追加項としてフォース項に追加することを特徴とする。
【数５０】
【００７２】
また、本発明は、前記境界条件導入工程（Ｓ６）にあって、前記各要素の自然境界を通過する熱流束がゼロでない場合に、節点領域の範囲と合致した境界熱流束の節点への配分を行うことを特徴とする。
【００７３】
また、本発明は、前記境界条件導入工程（Ｓ６）にあって、前記各要素の境界での荷重、体積力、質量がゼロでない場合に、節点領域の範囲と合致した荷重、体積力、質量の節点への配分を行うことを特徴とする。
【００７４】
そして、本発明は、有限要素法を用いた解析演算をコンピュータに実行させるための有限要素法を用いた解析演算プログラムであって、
前記コンピュータに、
解析対象の解析領域を選定する選定工程（Ｓ１）と、
前記解析領域を計算対象として複数の要素に分割する分割工程（Ｓ２）と、
前記複数の要素のうちの、ある要素（例えばｅ）にある節点（例えばＰ，要素ｅでの節点番号を例えば１とする）についてGalerkin重み関数（Ｗ）をかけて要素積分をし、各要素のマトリクスを作成する要素マトリクス作成工程（Ｓ３）と、
Galerkin重み関数（Ｗ）と一般関数（例えばＱ，ｆ）との積からなり、Galerkin有限要素法による二階微分項の離散化結果を基に定義された節点領域の考え方を導入し、要素の代表値を用いた一般関数項を、積分する一般関数項積分工程（Ｓ４）と、
前記ある節点（例えばＰ）周りの領域（例えばΩ）の要素における各要素のマトリクスの和と、前記一般関数項を積分した値の和とに基づき、連立方程式を作成する連立方程式作成工程（Ｓ５）と、
前記連立方程式に境界条件の導入をする境界条件導入工程（Ｓ６）と、
前記連立方程式を演算して数値解を得る演算工程（Ｓ７，Ｓ８）と、を実行させることを特徴とする。
【００７５】
なお、上記カッコ内の符号は、図面と対照するためのものであるが、これは、発明の理解を容易にするための便宜的なものであり、特許請求の範囲の構成に何等影響を及ぼすものではない。
【発明の効果】
【００７６】
本発明によると、所定の重み関数と一般関数との積からなる一般関数項を積分する一般関数項積分工程にあって、Galerkin有限要素法による二階微分項の離散化結果を基に定義された節点領域の考え方を導入し、要素代表値（例えば幾何学中心での値や、節点座標平均位置での値、要素平均値など）を用いた一般関数項を積分し、節点領域の大きさに見合った値を取り入れることができる。これにより、メッシュ・パターンによる解析誤差を低減することができ、品質が低い要素（品質の低いメッシュ・パターン）からも高精度な数値解を得ることが可能となる。このため、メッシュ・パターンの品質を高める作業（メッシュ・パターンを正三角形、正四面体、正方形などに近づける作業）を不要とすることができるものでありながら、アルゴリズム（プログラム）の修正による計算量の増加が殆んどないので、コンピュータによる演算処理の増大も防止することができ、総じて解析時間の短縮化を図ることができる。
【図面の簡単な説明】
【００７７】
【図１】二次元三角形における節点Ｐの積分領域Ωを示す図。
【図２】線形三角形要素を示す図。
【図３】一般三角形要素での節点領域を示す図。
【図４】直角三角形要素での節点領域を示す図。
【図５】鈍角三角形要素での節点領域を示す図。
【図６】温度分布の等値線を示す図。
【図７】従来改善方法２による節点領域の修正を説明する図。
【図８】線形直角四面体要素を示す図。
【図９】四面体要素による解析結果を示す図で、（ａ）はメッシュを示す図、（ｂ）は（ａ）のある断面上の温度分布を示す図。
【図１０】従来改善方法による四面体要素での節点領域を示す図。
【図１１】二次元四角形における節点Ｐの積分領域を示す図。
【図１２】長方形要素での積分結果を示す図。
【図１３】平行四辺形要素での積分結果を示す図。
【図１４】二節点が重なる四辺形要素での積分結果を示す図。
【図１５】三次元形体における節点Ｐの積分領域を示す図。
【図１６】長方体要素を示す図で、(ａ)は要素形状を示す図、(ｂ)は積分結果を示す図。
【図１７】平行四辺形を押し出した六面体要素を示す図で、(ａ)は要素形状を示す図、(ｂ)は積分結果を示す図。
【図１８】正三角柱要素を示す図で、(ａ)は要素形状を示す図、(ｂ)は積分結果を示す図。
【図１９】直角三角柱要素を示す図で、(ａ)は要素形状を示す図、(ｂ)は積分結果を示す図。
【図２０】鋭角三角形要素での節点領域を示す図。
【図２１】直角三角形要素での節点領域を示す図。
【図２２】鈍角三角形要素での節点領域を示す図で、（ａ）は面積が正の部分を示す図、（ｂ）は面積が負の部分を示す図。
【図２３】自然境界条件及び線上に分布するソース項の処理を説明する図で、（ａ）は節点３の角度≦π／２の場合の図、（ｂ）は節点３の角度＞π／２の場合の図。
【図２４】四面体における線上分布を説明するための図。
【図２５】自重による変位問題を説明するための図。
【図２６】板の有限要素モデルを示す図。
【図２７】二次元解析メッシュのパターンを示す図で、（ａ）は直角三角形メッシュの図、（ｂ）は鈍角三角形メッシュの図、（ｃ）は細長い三角形メッシュの図。
【図２８】二次元熱伝導解析の誤差を、従来の有限要素法（ＧＦＥ）、改善方法１（ＮＤＦＥ）、本発明の方法（ＣＮＤＦＥ）について示す図で、（ａ）は直角三角形メッシュの場合の誤差を示す図、（ｂ）は鈍角三角形メッシュの場合の誤差を示す図、（ｃ）は細長い三角形メッシュの場合の誤差を示す図。
【図２９】荷重、体積力、質量が係る三角形メッシュと拘束条件を示す図。
【図３０】従来有限要素法を三角形メッシュに適用した内圧での径方向変位の解析結果を示す図。
【図３１】従来有限要素法を三角形メッシュに適用した回転での径方向変位の解析結果を示す図。
【図３２】本発明を三角形メッシュに適用した内圧での径方向変位の解析結果を示す図。
【図３３】三次元解析メッシュのパターンを示す図で、（ａ）は立方体を六個の四面体に分割したメッシュの図、（ｂ）は長方形を十四個の四面体に分割したメッシュの図、（ｃ）は（ｂ）の四面体をさらに四個の四面体に分割したメッシュの図。
【図３４】三次元熱伝導解析の誤差を、従来の有限要素法（ＧＦＥ）、本発明の方法（ＣＮＤＦＥ）について示す図で、（ａ）は立方体を六個の四面体に分割したメッシュの場合の誤差を示す図、（ｂ）は長方形を十四個の四面体に分割したメッシュの場合の誤差を示す図、（ｃ）は（ｂ）の四面体をさらに四個の四面体に分割したメッシュの場合の誤差を示す図。
【図３５】四角柱の四面体要素分割を示す図。
【図３６】従来有限要素法を四面体メッシュに適用した自重での変位の解析結果を示す図。
【図３７】本発明を四面体メッシュに適用した自重での変位の解析結果を示す図。
【図３８】Ｔ１を含む項の節点領域の定義方法を説明するための図。
【図３９】Ｔ１を含まない項の節点領域の定義方法を説明するための図で、（ａ）は正の面積を示す図、（ｂ）は負の面積を示す図。
【図４０】四辺形要素の節点位置を説明するための図。
【図４１】四辺形メッシュの種類を示す図で、（ａ）は正方形のメッシュ、（ｂ）は３，４，５，８要素からなるメッシュ、（ｃ）は菱型のメッシュ、（ｄ）は平行四辺形のメッシュ、（ｅ）台形のメッシュ、（ｆ）直角台形のメッシュ、（ｇ）直角三角形型台形のメッシュ、（ｈ）平角四辺形のメッシュ、（ｉ）低平角四辺形のメッシュ、（ｊ）高平角四辺形のメッシュ。
【図４２】基本境界条件での解誤差比較を示す図。
【図４３】基本・自然境界条件での解析誤差比較を示す図。
【図４４】基本境界条件での節点数と解析誤差の相関図。
【図４５】基本・自然境界条件での節点数と解析誤差の相関図。
【図４６】荷重、体積力、質量が係る四角形メッシュと拘束条件を示す図。
【図４７】従来有限要素法を四角形メッシュに適用した内圧での径方向変位の解析結果を示す図。
【図４８】本発明を三角形メッシュに適用した内圧での径方向変位の解析結果を示す図。
【図４９】節点１とｉの温度に関わる節点領域を示す図。
【図５０】係数Ｄｊｋの計算方法を説明する図。
【図５１】節点１以外の温度に関わる節点領域を示す図。
【図５２】六面体メッシュを示す図で、（ａ）は長方体の図、（ｂ）は６，８，１０，１６要素なるメッシュの図、（ｃ）は菱型柱の図、（ｄ）は平行四辺柱の図、（ｅ）は台形柱の図、（ｆ）は直角三角形型台形の図、（ｇ）平角柱の図、（ｈ）は低平角柱の図。
【図５３】六面体要素での解析誤差を示す図で、（ａ）は基本境界条件での解析誤差を示す図、（ｂ）は基本・自然境界条件での解析誤差を示す図。
【図５４】六面体要素での節点数と解析誤差の相関を示す図で、（ａ）は基本境界条件での解析誤差を示す図、（ｂ）は基本・自然境界条件での解析誤差を示す図。
【図５５】五面体メッシュを示す図で、（ａ）は直角三角柱の図、（ｂ）は鈍角三角柱の図、（ｃ）は細長い三角柱の図、（ｄ）は方向性を持つ直角三角柱の図。
【図５６】五面体要素での解析誤差を示す図で、（ａ）は基本境界条件での解析誤差を示す図、（ｂ）は基本・自然境界条件での解析誤差を示す図。
【図５７】本発明に係る有限要素法の解析工程を示すフローチャート。
【発明を実施するための形態】
【００７８】
以下、本発明に係る実施の形態を詳細に説明する。なお、以下の説明においては、大まかに、メッシュ・パターンが二次元三角形の場合と三次元四面体の場合と二次元四角形の場合とに分けて説明する。また、以下の説明においては、有限要素法の解析方法について説明するが、一般的にはこの解析方法は解析プログラム化され、コンピュータによって演算処理されるものであり、言い換えると、図５７に示す各工程は、コンピュータに有限要素法を用いた解析演算を実行させるための解析演算プログラムということになる。
【００７９】
まず、有限要素法を用いた解析方法の流れを、図５７に沿って大まかに説明する。有限要素法により解析対象を数値解析する際は、図５７に示すように、ステップＳ１において、解析対象（どのような形状のものであっても良い）の解析領域を選定し（選定工程）、ステップＳ２において、解析領域を計算対象として複数の要素に分割する（分割工程）。続いてステップＳ３において、要素ｅの各節点（二次元の線形三角形の場合は節点１，２，３、三次元の線形四面体の場合は節点１，２，３，４、二次元の線形四角形の場合は節点１，２，３，４）について重み関数Ｗをかけて要素積分をし、各要素のマトリクス（例えば二次元三角形の場合は上述した式（３）の左辺、三次元四面体の場合は式（１０）の左辺、二次元四角形の場合は上述した式（３）’の左辺や式（１５）の左辺、三次元六面体・三次元五面体の場合は式（２４）の左辺）を作成する（要素マトリクス作成工程）。
【００８０】
次に本発明の要部となるステップＳ４において、Galerkin重み関数Ｗと一般関数（例えば熱問題等であればＱ、荷重、体積力、質量問題であればｆ）との積からなる一般関数項を積分する（一般関数項積分工程）。このステップＳ４において例えば熱問題等を扱うときは、ポアソンの方程式のソース項を積分する（例えば二次元三角形の場合は後述する式（３３）の右辺、三次元四面体の場合は後述する式（３６）の右辺、二次元四角形の場合は後述する式（４６）の右辺、三次元六面体・五面体の場合は後述する式（７２）の右辺）（ソース項積分工程）。また、例えば荷重、体積力、質量が係る場合には、二次元三角形の場合は式（３３）’、三次元四面体の場合は式（３６）’、二次元四角形の場合は式（７１）、三次元六面体・五面体の場合は式（９６）を用いて、荷重、体積力、質量の要素積分、即ち、それぞれ式（３３）’の一項目、式（３６）’の一項目、式（７１）の一項目、式（９６）の一項目を取って代わる（フォース項積分工程）。
【００８１】
その後、ステップＳ５において、節点Ｐ周りの要素領域Ωにおける各要素のマトリクスの和と、一般関数項を積分した値の和とに基づき、連立方程式を作成する（連立方程式作成工程）。
【００８２】
この際、本発明では、詳しくは後述するように、ステップＳ６において、節点領域の範囲と合致した境界条件を、各節点周りについて作成した連立方程式のそれぞれに導入する（境界条件導入工程）。
【００８３】
そして、ステップＳ７において、連立方程式を例えば反復法等により演算して解き（演算工程）、ステップＳ８において、関連物理量の計算を行って数値解を得て（演算工程）、以上の工程によって解析を終了する。
【００８４】
［二次元三角形の解析について］
以上のような有限要素法を用いた解析方法において、本発明にあっては、ステップＳ４の一般関数項積分工程における一般関数項を、要素ｅの代表値（例えば幾何学中心での値や、節点座標平均位置での値、要素平均値など）を用いて一般関数項について計算する。
【００８５】
具体的には、ポアソンの方程式の場合には、式（１）左辺の楕円型演算子については従来の有限要素法により離散化をし、右辺のソース項については、従来の有限要素法による離散化
【数５１】
に追加項
【数５２】
を加える。すなわち、従来の式（３）の代わりに、式（３３）を提案する。
【数５３】
【００８６】
ここで、ＱＧは要素のソース代表値であり、ＮＤは対象節点の要素ｅにおける節点領域の面積である。ＮＤの定義は、鋭角三角形要素においては図２０に、直角三角形要素においては図２１に、鈍角三角形要素においては図２２に、示されるように行う。図中のＮＤの添え字で所属節点を表わし、点Ｃは外心であり、点Ｉ、Ｊ、Ｋは中点である。また、鈍角三角形要素の場合には、外心が要素の外に出ているため、図２２（ｂ）に示す節点２、節点３が所有する節点領域−ＮＤ２、−ＮＤ３が負になる。それぞれの節点領域は最終的に次式により算出される。
【数５４】
【００８７】
また、ステップＳ６の境界条件導入工程にあっては、以下のように境界条件を導入する。即ち、熱流束がゼロではない自然境界については図２３に示すように処理する。ここで、節点１と節点２を結ぶ線を境界とすると、節点３での内角がπ／２以下の場合（つまり鋭角の場合）には、図２３（ａ）のように点１と点Ｉとの間を通過し流出する熱量
【数５５】
を節点１の方程式に、点Ｉと点２との間の熱通過量
【数５６】
を節点２の方程式に加える。ここで、
【数５７】
は線積分であり、ｎは外向き法線方向である。
【００８８】
節点３での内角がπ／２を越えた場合（つまり鈍角の場合）には、図２３（ｂ）のように点１と点Ｉ１間の熱通過量
【数５８】
を節点１の方程式に、点Ｉ２と点２との間の熱通過量
【数５９】
を節点２の方程式に、点Ｉ１と点Ｉ２との間の熱通過量
【数６０】
を節点３の方程式に加える。
なお、線上に分布するソース項に対しても、図２３に示す自然境界条件の分け方と同様に節点の方程式に寄与する。
【００８９】
以上説明した本実施の形態による有限要素法を用いた二次元三角形の解析方法によれば、以下の特徴を持っている。
（１）節点領域の大きさに見合ったソース量が取り入れられたため、節点レベルで保存法則を２次精度で満足する。
（２）要素の品質が低いほど修正効果が大きい。
（３）いかなる要素形状においても、要素外のソース項を取り入れることがない。
（４）いかなる要素形状においても、要素の３つの節点が所有する節点領域の合計は
【数６１】
になるので、ソース項に追加項を入れても、要素レベルでも保存法則が満たされる。
（５）楕円型演算子に対する最適な離散化手法としての有限要素法の優位性がそのまま保持される。
（６）アルゴリズムの修正による計算量の増加は殆どない。
【００９０】
［三次元四面体の解析について］
また、三次元四面体の解析にあっても、前述のような有限要素法を用いた解析方法において、ステップＳ４の一般関数項積分工程における一般関数項を、要素ｅの代表値を用いて一般関数項について計算する。
【００９１】
具体的には、式（８）左辺の楕円型演算子については従来の有限要素法により離散化をし、右辺のソース項については追加項
【数６２】
を加えて、つまり従来の式（１０）の代わりに、式（３６）を提案する。
【数６３】
【００９２】
ここでも、ＱＧはソース代表値であり、ＮＤは対象節点の要素ｅにおける節点領域の体積である。節点領域ＮＤは以下の方法で決める。対象節点の要素ｅにおける節点番号を１とし、そこに局所座標系の原点を置くと、式（３６）左辺の積分結果を次式右辺のような係数と熱流束との掛け算の形に整理することができる。
【数６４】
【００９３】
ここで、
【数６５】
は図１０に示す要素の辺中点Ｉ、Ｊ、Ｋでの節点１から節点２、３、４に向かう熱流束を二次精度で近似するものであり、係数Ｓ１２、Ｓ１３、Ｓ１４はそれぞれの熱流束の通過面積にあたるものであり、その二つの添え字で熱流束の向きを表わす。節点１の節点領域の体積は、熱流束の通過面積を底面とし、節点１と中点間の距離を高さとする角錐の体積の合計として、次式のように求める。
【数６６】
【００９４】
また、熱流束の通過面積Ｓｉｊ（ｉ，ｊ＝１，２，３，４、ｉ≠ｊ）を表わす汎用式は、以下のように導出される。まず、式（３７）の左辺を有限要素法のアプローチにより展開すると次式が得られる。
【数６７】
【００９５】
ここで、要素体積Ｖは次式のように求められる。
【数６８】
【００９６】
式（３９）中の係数は以下のようなものである。
【数６９】
【００９７】
式（３９）よりλ（Ｔｉ−Ｔｊ）の係数は
【数７０】
になることが分かる。節点ｉとｊとの間の距離をｌｉｊで表わすと、式（４２）を得ることができ、λ（Ｔｉ−Ｔｊ）の係数を式（４３）のように書き換えることができる。
【数７１】
【数７２】
【００９８】
式（３７）の右辺を参照すると、節点ｉと節点ｊとの間に伝わる熱流束の通過面積は次式により求めることができる。
【数７３】
この式（４４）を用いて、各節点の節点領域を式（３８）と同じ考え方で算出する。
【００９９】
また、ステップＳ６の境界条件導入工程にあっては、以下のように境界条件を導入する。即ち、熱流束がゼロではない自然境界では、従来の有限要素法による処理結果
【数７４】
に、追加項
【数７５】
を加えて、
【数７６】
を境界条件として式（３６）に導入する。ここで、ｓは要素ｅに所属する境界を表わし、Ｓはｓの面積であり、添え字ＳＧはｓの幾何中心を意味し、ｎは外向き法線方向で、ＮＤはｓを上述の二次元三角形解析の方法で分割した節点領域の面積である。
【０１００】
なお、線上に分布するソース項Ｑに対しても、まず、線を共有する各要素の表面三角形に均等に分ける。例えば、図２４に示す節点１と節点２とを結ぶ線の場合には、それを有する要素表面三角形は、三角形１２４、１２３、１２３、１２５と四つがある。三角形１２３は二要素に属するため二つとして数える。各三角形にＱ／４を配分する。さらに、三角形内の三節点への配分については、上述した二次元三角形解析における線上に分布するソース項の分け方と同方法で行う。面上に分布するソース項についても、二次元三角形解析の方法で同様に配分する。
【０１０１】
以上説明した本実施の形態による有限要素法を用いた三次元四面体の解析方法によれば、以下の特徴を持っている。
（１）節点領域の大きさに見合ったソース項が分配されるため、節点レベルにおいては保存法則を従来法より高い精度で満たす。
（２）要素の品質が低いほど修正効果が大きい。
（３）任意の四面体形状においても、要素外のソース項を取り入れることがない。
（４）任意の四面体形状においても、要素の４節点に所属する節点領域の合計が要素体積に等しいため、追加項を入れても要素レベルで保存法則が満たされる。
（５）楕円型演算子に対する最適な離散化手法としての有限要素法の優位性がそのまま保持される。
（６）アルゴリズムの修正による計算量の増加は殆どない。
【０１０２】
［二次元三角形解析及び三次元四面体解析における他分野の解析］
なお、上述した二次元三角形解析と三次元四面体解析とでは、熱伝導問題を例に本解析手法の応用について説明してきたが、勿論、本手法はポアソンの方程式が記述する現象の全てにそのまま展開できる。例えば、流れ問題、電気・磁気問題、弾性変形問題、物質拡散問題などが応用できる分野として挙げられる。
【０１０３】
また、本発明の基本的な考え方は、二次精度の保存則の観点から楕円型演算子の離散化結果におけるコントロール・ボリュームを見出し、ソース項の積分領域をそれに合わせたものである。従って、この考え方は、保存則から導出された楕円型演算子を含む各種の支配方程式の有限要素解析に容易に応用できる。すなわち、物理量の発生・消失と流出・流入を算出するとき、その積分領域について、本発明で導出した節点領域の考え方を導入して積分すればよい。例えば、非定常移流拡散解析において、非定常項と移流項の離散化に対して本発明におけるソース項の積分領域の考え方が適用できる。
【０１０４】
また、荷重、体積力と質量（フォースと総称する）が係わる従来の有限要素法による解析においては、それらの扱い方はポアソンの方程式のソース項と同じである（非特許文献２参照）（すなわち、Galerkin重み関数とかけて要素積分をする。）ため、従来の技術で述べたような、要素形状と関係なく均等的に節点方程式に分配される問題も存在するが、これらの問題は、本実施の形態で説明したソース項の扱い方で改善することもできる。
【０１０５】
従来の技術での問題点を以下の例で説明する。図２５に示す縦に置かれた幅がＬ、高さがＨ、厚さが１の板の自重による変位について考察する。密度をρ、重力加速度をｇ、弾性係数をＥで表わす。この際、問題の本質をより明確にするため、平面応力問題とし、ポアソン比をν＝０とする。板全体の高さの変化量をＶで表わし、上向きの座標をｙとすると、応力がσ＝−（Ｈ−ｙ）ρｇとなり、ひずみがε＝σ／Ｅとなるため、次式の厳密解が得られる。
【数７７】
【０１０６】
一方、板を図２６のように２つの三角形要素に分割する。節点３と節点４に寄与されるｙ方向のフォース量をそれぞれｆ３とｆ４で表わすと、従来の有限要素法では、式（５）に示すように要素の重さの３分の１ずつが各節点に分配されるため、
【数７８】
となる。ここで、節点４が二つの要素を所有するため、分配された荷重は節点３の２倍になっている。解析結果として、節点３と節点４のｙ方向の変位ν３とν４が次のようになる。
【数７９】
これらの結果は、厳密解からずれていることが確認できる。
【０１０７】
［線上に分布する場合］
二次元三角形解析においては、図２３に示すように、節点３での内角がπ／２以下の場合には、
【数８０】
をそれぞれ節点１、２の方程式に加える。ここで、ｆは荷重の任意成分、体積力の任意成分、質量を意味する。節点３での内角がπ／２を越えた場合には、
【数８１】
をそれぞれ節点１、２、３の方程式に加える。
【０１０８】
三次元四面体解析においては、まず、線を共有する各要素の表面三角形に均等に分ける。例えば、図２４に示す節点１と節点２を結ぶ線の場合には、四つの表面三角形に、それぞれｆ／４を配分する。さらに、三角形内の三節点への配分については、上述した二次元三角形解析における線上に分布するソース項の分け方と同方法で行う。
【０１０９】
［面上に分布する場合］
二次元三角形、あるいは三次元四面体解析において、面ｓ上に分布する荷重、体積力、質量については、式（３３）右辺のソース項と同様に修正し、修正後の対象節点に寄与するものは
【数８２】
になる。ここで、Ｓはｓの面積であり、添え字Ｇで要素代表値（例えば幾何学中心での値や、節点座標平均位置での値、要素平均値など）を表わす。
【０１１０】
［空間に分布する場合］
空間に分布する荷重、体積力、質量については、式（３６）右辺のソース項と同様に修正し、修正後のものは
【数８３】
になる。
【０１１１】
［二次元三角形解析及び三次元四面体解析における本発明の効果］
本解析方法の導入により解析精度が改善されたことを以下の事例により示す。
【０１１２】
前記板の自重による変位問題に、節点領域の考え方を適用し、式（３３）’と図２１を参照すると、節点３と節点４に寄与されるｙ方向のフォース量は次のようになることが分かる。
【数８４】
【０１１３】
ここで、節点４に寄与されたフォース量は、二つの要素の合計になる。解析結果として、節点３と節点４のｙ方向の変位は次のようになる。
【数８５】
厳密解と一致しており、解析精度が改善されたことが確認できる。
【０１１４】
二次元三角形問題については、式（６）と（７）により定義された熱伝導問題を図２７のメッシュ・パターンを用いて解析した最大誤差を図２８に示す。この際λ＝１とし、ソース項の要素代表値として幾何学中心を採用すると、最大誤差は次式により定義される。
【数８６】
【０１１５】
ここで、ＴＡｉ：節点ｉでの厳密解、ＴＮｉ：節点ｉでの数値解、ＴＡｍａｘ：厳密解の最大値、である。図２８の横軸はメッシュの代表長さであり、同じパターンのメッシュを２倍ずつ拡大したものである。中には、図２８（ａ）の１／８、（ｂ）の１／８、（ｃ）の１／５が、図２７のメッシュそのままを使用したものである。「ＧＦＥ」で従来の有限要素法、「ＮＤＦＥ」で上記改善方法１、「ＣＮＤＦＥ」で本発明、での結果を示す。図２８より、本提案方法は、いずれのメッシュ・パターンにおいても解析誤差が従来の有限要素法の３分の１以下になったことが分かる。また、上記改善方法１は、直角三角形メッシュの場合にはよい精度が得られたが、細長いメッシュを使用すると、誤差が著しく増大することが認められた。
【０１１６】
また、荷重、体積力と質量が係わる有限要素法による解析の事例として、薄い円筒に、（１）均等な内圧を加えた変形問題、（２）円筒中心軸を回転中心とした等速回転問題、を取り上げた。図２９に解析メッシュを示しており、直角等辺三角形のシェル要素によりメッシュ分割を行った。この際、軸方向の中央に位置する節点（●で示す節点）に円周方向と軸方向の拘束条件を与えた。
【０１１７】
図３０に従来の有限要素法を用いた内圧による径方向の変位を、図３１に従来の有限要素法を用いた回転による径方向の変位を示す。点線で初期形状を示し、灰色で膨張方向の変位を、黒色で縮小方向の変位を表わす。各節点は同じ変位をする実現象と異なって、節点周辺の局所メッシュ・パターンによって、変位差が現れた。その原因は、図６の説明で述べた不整合性の問題と共通である。図３２に本発明での提案方法（上記「面上に分布する場合」）による内圧での径方向の変位結果を示しており、要素形状に応じた節点荷重を与えた結果、均一変位になっていることが確認できる。回転問題についても、解析原理が同じであるため、遠心力を同じ考え方で各節点に分配すれば、均一変位が得られることは明確的である。
【０１１８】
三次元四面体問題については、式（１３）と（１４）により定義された熱伝導問題を図３３のメッシュ・パターンを用いて解析した最大誤差を図３４に示す。この際λ＝１とし、ソース項の要素代表値として幾何学中心値を採用した。メッシュ・パターンとして、図３３（ａ）は一個の立方体を六個の四面体に、（ｂ）は一個の長方体を十四個の四面体に、（ｃ）は（ｂ）の四面体をさらに四個の四面体に分割したものである。本発明の解析方法は、いずれのメッシュ・パターンにおいても解析精度が改善され、特に品質の悪い要素では誤差が従来の有限要素法の２分の１以下になったことが分かる。
【０１１９】
荷重、体積力と質量が係わる有限要素法による解析の事例として、図３５の四角柱の自重による変形問題を取り上げ、本発明の改善効果を示す。図３６には従来の有限要素法による高さ方向の変位結果を、図３７には本発明、すなわち、式（３６）’を導入した高さ方向の変位結果を、示す。上面における４節点が同じ変位をする実現象に対して、本発明の解析結果においては、４節点の変位差は、従来の有限要素法に比べ大幅に改善されたことが確認できる。改善の理由については、上面一番左の節点を例に説明する。従来の有限要素法においては、当節点に要素自重の４分の１しか寄与されていないため、変位の絶対値が小さくなっている。本発明では、当節点に式（３６）’に基づいて要素自重の２分の１が寄与されるため、変位の改善効果が得られた。
【０１２０】
［二次元四角形の解析について］
ついで、二次元四角形の解析について説明する。二次元四角形の解析において保存則を満たさない問題を改善するため、簡単に実行できるソース項の保存型離散化手法を提案する。式（３）’の左辺については有限要素法のままで離散化をし、右辺については追加項
【数８７】
を加えて、次式により有限要素解析を行う。
【数８８】
【０１２１】
ここで、ＮＤＰは節点Ｐの要素ｅにおける節点領域(Nodal
Domain)と称するものであり、ＱＯは点Ｏでのソース値である。ＮＤＰが式（１５）左辺の積分から導かれたコントロール・ボリュームと面積上で等しくなれば、保存則を２次精度で満足することができ、メッシュ依存性を改善し、特に、低品質要素での精度向上効果が期待できる。考え方としては、ＱＯの代わりに要素ソース項のほかの代表値、例えば重心でのソース値や要素平均値などを式（４６）に入れることも考えられる。また、式（３）’右辺の要素積分を単純にＮＤＰＱＧで近似する選択もある。
【０１２２】
一般形状の四辺形要素において、節点領域を以下の普遍性のある方法により算出する。節点Ｐの要素ｅにおける局所番号を１とすると、式（１５）左辺の要素積分の結果を次の汎用式に整理することができる。
【数８９】
係数Ａｉは四節点の座標のみで決まる定数である。有限要素法の性質より、次式が成り立つ。
【数９０】
【０１２３】
［二次元四角形における節点領域の算出方法１］
ｉ，ｊ，ｋで、それぞれ節点２、３、４の中の異なる一節点を表わし、順序には制限がないこととする。係数Ａ２，Ａ３，Ａ４が持つ符号の全ての組み合わせと式（４８）を考えると、式（４７）を次の節点温度差の形に書き換えることができる。
【数９１】
【０１２４】
すなわち、Ａ１のみが正の場合には式（４９−１）を、Ａ１，Ａ２，Ａ３，Ａ４の中の二つの係数が正の場合には式（４９−２）、もしくは（４９−３）を、三つの係数が正の場合には式（４９−４）を採用する。また、式（４９−１）〜（４９−４）の右辺の各項の係数が正の値を持つことが分かる。
【０１２５】
式（４９−１）〜（４９−４）右辺の各項を、Ｔ１を含むものと含まないものに分類しておく。Ｔ１を含む項の展開については、λＢ１ｊ（Ｔ１−Ｔｊ）を例に説明する。ここで、Ｂ１ｊでλ（Ｔ１−Ｔｊ）項の係数を表わす。Ｌ１ｊで節点1と節点ｊ間の線分を表わすと、次式が得られる。
【数９２】
【０１２６】
式（５０）右辺は線分Ｌ１ｊの中点での節点1から節点ｊに向かう２次精度の熱流束により、Ｌ１ｊの中点から垂直方向に長さＢ１ｊＬ１ｊを伸びた線分を通過する熱量であることが分かる。この線分をＢ１ｊＬ１ｊで表わす。保存則の観点から、この項に相当するソース項の領域は、節点1、線分Ｌ１ｊの中点、線分Ｂ１ｊＬ１ｊの端点、といった三点から成す三角形となる。この三角形の面積をＮＤ１ｊで表わし次式により計算することができる。
【数９３】
【０１２７】
図３８には二つの例としてＮＤ１３とＮＤ１４を示す。Ｔ１を含まない項の展開については、λＢｉｊ（Ｔｉ−Ｔｊ）を例に説明する。節点ｉと節点ｊ間の線分をＬｉｊで表わすと、次式が得られる。
【数９４】
【０１２８】
式（５２）右辺は線分Ｌｉｊの中点での節点ｉから節点ｊに向かう２次精度の熱流束により、Ｌｉｊの中点から垂直方向に長さＢｉｊＬｉｊを伸びた線分を通過する熱量であることが分かる。この線分をＢｉｊＬｉｊで表わす。この項に相当するソース項の領域は、節点1、Ｌｉｊの中点（ｍと記す）、線分ＢｉｊＬｉｊのもう一つの端点（ｌと記す）といった三点から成す三角形となる。この三角形の面積をＮＤｉｊで表わし、次式により計算する。
【数９５】
式（５３）右辺の符号については、点1から中点ｍに向かうベクトルと流束ベクトルとの角度αが｜α｜＜π／２にある、すなわち、熱流束が三角形から流出する場合には正とする。αが｜α｜≧π／２にある、すなわち、熱流束が三角形に流入する場合には負とする。図３９には二つの例としてＮＤ２３＞０とＮＤ２４≦０を示す。
【０１２９】
最終的に、対象節点１の節点領域の面積は、係数Ａ２，Ａ３，Ａ４の符号に合わせて式（５４−１）〜（５４−４）中の一式により算出する。
【数９６】
【０１３０】
［二次元四角形における節点領域の算出方法２］
係数Ａ２，Ａ３，Ａ４の符号と大きさに関係なく、式（４９−１）、（５０）、（５１）と（５４−１）によりＮＤ１を求める。この際、ＮＤ１ｊ（ｊ＝２，３，４）はＢ１ｊ、すなわち−Ａｊと同じ符号になる。
【０１３１】
［本二次元四角形の解析手法の整合性、唯一性、正規性］
整合性とは、節点領域とコントロール・ボリュームと同面積であることを意味する。整合性が成立すれば、節点レベルでの保存則が２次精度で満足されることになり、精度改善効果が理論上で保証される。
【０１３２】
図１２の長方形要素に対しては、式（１６）を次式に書き換えておく。
【数９７】
【０１３３】
上記算出方法１の整合性については、ｈｙ／√２≦ｈｘ≦√２・ｈｙの場合には、Ｔ１以外の係数が負になっているため、式（４９−１）に当てはめ、次式が得られる。
【数９８】
【０１３４】
その結果、節点領域として次式が得られる。
【数９９】
【０１３５】
また、ｈｘ＞√２・ｈｙの場合には、Ｔ２の係数が正に、Ｔ３，Ｔ４の係数が負になる。Ａ１＞｜Ａ３｜であるため、式（４９−２）に当てはめると次式が得られる。
【数１００】
【０１３６】
その結果、節点領域として次式が得られる。
【数１０１】
【０１３７】
なお、ｈｙ＞√２・ｈｘの場合についてもＮＤ１＝Ｓ／４であることが証明できるため、ＮＤ１＝ＣＶ１の関係式がｈｙ／ｈｘの比例に関係なく成立していることが分かる。
【０１３８】
上記算出方法２の整合性については、式（４９−１）を採用したので、その結果は式（５７）となる。
【０１３９】
ほかの節点については、要素形状の対称性を考慮すると、ＮＤ２＝ＮＤ３＝ＮＤ４＝Ｓ／４を得ることができる。結局、長方形要素においては、式（４６）右辺の追加項がゼロとなり、保存型離散化がGalerkin有限要素法と同じものになる。
【０１４０】
さらに、図１３の平行四辺形要素、図１４の二節点が重なる四辺形要素に適用すると、上記表２に示す節点領域が得られる。以上の三種類の四辺形要素に対して、コントロール・ボリュームと完全一致していることが確認できる。しかし、一般形状の四辺形要素については、コントロール・ボリュームが明確になっていないため、代わりに後述の数値解析により精度改善効果について検証する。
【０１４１】
節点領域の唯一性とは、式（５４−１）〜（５４−４）により求められたＮＤ１の結果は、「ｉ，ｊ，ｋへの局所節点番号の与え方によって変らない」、「算出方法２を使っても変らない」、ことを意味する。例えば、ｈｘ＞√２・ｈｙの長方形要素に対して、式（５８）においては、ｉ＝２、ｊ＝３、ｋ＝４としたが、代わりにｉ＝２、ｊ＝４、ｋ＝３の選択も、次式により節点領域を求めることも考えられる。
【数１０２】
【０１４２】
その結果、次式の節点領域が得られる。
【数１０３】
これは、式（５９）の結果と同じであるため、算出方法１においては唯一性が証明された。算出方法２における唯一性は、式（５７）の結果により示されている。
【０１４３】
正規性とは、
【数１０４】
すなわち、式（４６）右辺に入れた追加項の同要素四節点についての合計がゼロになることを意味する。正規性が成立すれば、要素レベルでも、全解析領域でも保存則が満たされることが保証される。図１２、図１３、図１４の要素形状に対しては、上記表２の最後の行に示すように、
【数１０５】
になるため、正規性が保たれている。
【０１４４】
同様に図１３、図１４の要素形状についても、唯一性が保たれたことを理論的に証明できる。しかし、一般形状の四辺形要素においては、ヤコビアンが絡んでくるため、正規性と唯一性について理論上の議論は困難である。代わりに、下記の数値計算により検討する。四辺形要素の形状としては、節点１を図４０の太い●印の二通りの配置に対して、節点２を○印に置き、節点３を最下行以外の全点に置き、節点４を中心列を含めた左側の全点に置いた組み合わせで、できた各四辺形要素をテスト対象とした。これらの要素形状では実用解析で使用される要素形状をカバーできるといえる。この際、上記算出方法１においてはｉ，ｊ，ｋに許容される全ての局所節点番号の組み合わせを与えて数値解析により節点領域を算出した。その結果、面積がゼロ以下の形状と凹形状を除いて、上記算出方法１での全テストケースと上記算出方法２での全テストケースにおいて正規性と唯一性が満たされていることが確認できた。
【０１４５】
［二次元四角形における本発明の効果］
保存型離散化手法の精度改善効果を二つの熱伝導問題の解析事例により示す。この際、解析領域を正方形領域（０≦ｘ＜１、０≦ｙ≦１）に限定し、要素代表ソース値を四節点の平均位置での値とし、λ＝１とした。
【０１４６】
解析問題１としては、式（６３）のソース項に式（６４）の基本境界条件のみを加えたものである。
【数１０６】
【数１０７】
【０１４７】
その厳密解は次式になる。
【数１０８】
【０１４８】
解析問題２としては、式（６６）のソース項に式（６７）の基本境界条件と式（６８）の自然境界条件を加えたものである。
【数１０９】
【数１１０】
【数１１１】
【０１４９】
その厳密解は次式になる。
【数１１２】
【０１５０】
解析の数値誤差について次式により評価する。
【数１１３】
ここで、ＴＡｉ：節点ｉでの厳密解、ＴＮｉ：節点ｉでの数値解、ＴＡｍａｘ：厳密解の最大値、である。
【０１５１】
図４１に使用した１０パターンのメッシュを、対称性があるため０≦ｘ≦０．５、０≦ｙ≦０．５の部分のみ示す。「・」印で節点位置を、各図の表題の括弧内の数字で総節点数を表わす。図４２に基本境界条件での数値誤差、図４３に基本・自然境界条件での数値誤差を示す。これらの図より、本発明の手法（CSFEで表わす）では、正方形以外の全てのメッシュにおいて解析誤差が、従来の有限要素法（GFEで表わす）の２分の１以下に低減されたことが分かる。
【０１５２】
図４４と図４５に各メッシュにおける総節点数の逆数と解析誤差との相関を示しており、実線で二通りの正方形要素での結果を結んだ。これらの図より、従来の有限要素法にはメッシュ品質の影響が大きいが、保存型離散化手法にはメッシュ品質の影響が少なく、低平角四辺形を除けば、正方形とほぼ同程度の解析精度が得られた。
【０１５３】
結論として、線形四辺形要素によるPoissonの方程式の有限要素解析において、新たな節点領域の定義を考案し、２階微分項との整合性がとれたソース項の保存型離散化手法を提案することができた。従って、下記の（１）〜（４）の効果を得ることができた。
（１）節点方程式に寄与するソース量が、要素形状に関係する節点領域の大きさに対応させたことにより、節点レベルで保存則を２次精度で満足するようになった。しかも、要素レベルでも保存則が満たされている。
（２）鈍角要素形状においても、要素外のソース項を取り入れることがなく、ソース項が急激に変化する場合にも数値誤差を抑えることができた。
（３）要素の品質が低いほど、本手法の精度改善効果が大きい。
（４）幾つかのメッシュ・パターンを用いて数値解析を行った結果、従来の有限要素法に比べ解析精度が大幅に改善されることが確認できた。
【０１５４】
［二次元四角形解析における他分野の解析］
以上の二次元四角形解析においては、熱伝導問題を例に本発明の実施形態について説明してきたが、勿論、Poissonの方程式が記述する現象の全てにそのまま展開できる。例えば、流れ問題、電・磁気問題、弾性変形問題、物質拡散問題などが応用できる分野として挙げられる。
【０１５５】
また、本発明の基本的な考え方は、二次精度の保存則の観点から楕円型演算子の離散化結果におけるコントロール・ボリュームを見出し、ソース項の積分領域をそれに合わせたものである。従って、この考え方は、保存則から導出された楕円型演算子を含む各種の支配方程式の有限要素解析に容易に応用できる。すなわち、物理量の発生・消失と流出・流入を算出するとき、その積分領域について、本発明で導出した節点領域の考え方を導入して積分すればよい。例えば、非定常移流拡散解析において、非定常項と移流項の離散化に対して本発明におけるソース項の積分領域の考え方が適用できる。
【０１５６】
また、荷重、体積力、質量が係わる従来の有限要素法による解析においては、それらの扱い方はPoissonの方程式のソースと同じであるため、要素形状と関係なくほぼ均等的に節点方程式に分配される問題も存在する。これらの問題は本発明でのソース項の扱い方で改善することもできる。二次元（面上に分布する場合）、あるいは三次元問題（空間に分布する場合）の解析において、面ｓ上に分布する荷重、体積力、質量については、式（４６）の右辺のソース項と同様に追加項を加えて修正すると、対象節点に寄与する荷重、体積力、質量は
【数１１４】
となる。ここで、ｆは荷重の任意成分、体積力の任意成分、質量を意味する。
【０１５７】
荷重、体積力と質量が係わる有限要素法による解析の事例として、薄い円筒に、均等な内圧を加えた変形問題を取り上げた。図４６に解析メッシュを示しており、四角形のシェル要素によりメッシュ分割を行った。この際、●で示す節点に円周方向と軸方向の拘束条件を与えた。
【０１５８】
図４７に従来の有限要素法を用いた内圧による径方向の変位を示す。点線で初期形状を示し、灰色で膨張方向の変位を、黒色で縮小方向の変位を表わす。各節点は同じ変位をする実現象と異なって、節点周辺の局所メッシュ・パターンによって、変位差が現れた。図４８に本発明での提案方法（式（７１））による内圧での径方向の変位結果を示しており、要素形状に応じた節点荷重を与えた結果、変位のバラツキが大幅に改善されたことが確認できる。体積力と質量が係わる問題は、同じ原理が適用できるものである。
【０１５９】
［三次元六面体・五面体の解析について］
前述の「従来の有限要素法による三次元六面体及び三次元五面体の問題」で述べた保存則を満たさない問題を改善するため、簡単に実行できるソース項の保存型離散化手法を提案する。式（２３）の左辺については従来の有限要素法のままで離散化をし、右辺については追加項
【数１１５】
を加えて、次式により有限要素解析を行う。
【数１１６】
【０１６０】
ここで、ＮＤＰは節点Ｐの要素ｅにおける節点領域(Nodal Domain)と称するものであり、ＱＯは当要素に属する全節点の座標平均位置でのソース値である。ＮＤＰが式（２４）の左辺の要素積分から導かれたコントロール・ボリュームと体積上で等しくなれば、保存則を２次精度で満足することができ、メッシュ依存性を改善し、特に、低品質要素での精度向上効果が期待できる。また、考え方としてはＱＯの代わりに要素ソース項のほかの代表値、例えば重心でのソース値などを式（７２）に入れることも考えられる。また、式（２３）の右辺の要素積分を単純にＮＤＰＱＯで近似する選択もある。
【０１６１】
一般形状の六面体要素と五面体要素において、節点領域を以下の方法により算出する。節点Ｐの要素ｅにおける局所番号を１とすると、式（２４）左辺の要素積分の結果を次の汎用式に整理することができる。
【数１１７】
ここで、Ｎは一要素の節点数であり、六面体要素ではＮ＝８、五面体要素ではＮ＝６となる。また、係数Ａｉは節点座標のみで決まる定数である。有限要素法の性質より次式が成り立つ。
【数１１８】
【０１６２】
［三次元六面体・五面体における節点領域の算出方法１］
式（７４）を利用すると、式（７３）を次式に書き換えることができる。
【数１１９】
【０１６３】
ここで、Ｉ１ｉは節点１と節点ｉ間の距離であり、次式により算出する。
【数１２０】
【０１６４】
式(７５)の右辺の物理的な意味は、節点１と節点ｉの中点（Ｉ１ｉと記す）での節点１から節点ｉに向かう熱流束の２次近似λ（Ｔ１−Ｔｉ）／Ｉ１ｉにより、熱流束に垂直で中点Ｉ１ｉを通る面積（−Ａｉ）Ｉ１ｉを通過する熱量である。この流出熱量により節点１に割り当てられた熱源の体積を、図４９に示す節点１を頂点とし（−Ａｉ）Ｉ１ｉを底面積とした錐体とする。この錐体の体積をＮＤ１ｉで表わすと、次式が得られる。
【数１２１】
【０１６５】
各熱流束に割り当てる熱源の体積を合計したものを節点１の節点領域ＮＤ１とすると、次式が得られる。
【数１２２】
節点領域の求め方は、保存則によりある物理量の支配方程式を導出するときに使用されるコントロール・ボリュームの考え方と合致する。式（７８）の節点領域を式（７２）右辺に代入して、精度向上を図る。なお、係数Ａｉ（ｉ＝２，３，・・・，Ｎ）の中、正のものが存在するケースもあるが、この場合は、熱が錐体に流入することになっているためＮＤ１ｉが負になる。
【０１６６】
［三次元六面体・五面体における節点領域の算出方法２］
Ａｉ（ｉ＝１，２，・・・，Ｎ）中のＡｉ＞０のものをそれぞれＢｊ（ｊ＝１，２，・・・，Ｎ＋）で表わし、Ｎ＋でその数を表わす。また、Ａｉ≦０のものをそれぞれＣｋ（ｋ＝１，２，・・・，Ｎ−）で表わし、Ｎ−でその数を表わすと、式（７３）を次式で表わすことができる。
【数１２３】
【０１６７】
ここで、添え字ｉは要素に属する全節点を有限要素法で決めたルールに従って１から数えるものであるが、ｊは係数が正の節点のみを任意の順序で１から数え、ｋは係数が負の節点のみを任意の順序で１から数えるものである。従って、例えｉ，ｊ，ｋが同じ数字をとったとしても、ＴｉとＴｊ、ＴｉとＴｋとは必ずしも同一節点での温度ではなく、また、ＴｊとＴｋとはいつも異なる節点での温度である。式（７４−２）より次式が成り立つことが分かる。
【数１２４】
【０１６８】
各係数の値を棒の長さで表わすと、式（８０）の関係式を図５０で表現することができる。各係数間の繋ぎ目で区分すると、式（７３）を次式のように係数が正の節点温度と係数が負の節点温度との差の和に書き換えることができる。
【数１２５】
【０１６９】
ここで、Ｄｊｋは正の値のものでＢｊと−Ｃｋとが重なる部分である。Σはｊｋの組み合わせの合計を意味し、ｌｊｋは節点ｊとｋ間の距離である。式(８１)の右辺の物理的な意味は、節点ｊとｋの中点（Ｉｊｋと記す）での節点ｊから節点ｋに向かう熱流束の２次近似λ（Ｔｊ−Ｔｋ）／Ｉｊｋにより、熱流束に垂直で中点Ｉｊｋを通る面積Ｄｊｋｌｊｋを通過した熱量である。この熱量により節点１に割り当てられた熱源の体積は、図５１に示す節点１を頂点としＤｊｋｌｊｋを底面積とした錐体とする。この体積をＮＤｊｋで表わすと、次式が得られる。
【数１２６】
【０１７０】
式（８２）の右辺の［］内のものは、節点１から中点Ｉｊｋまでのベクトルと節点ｊから節点ｋまでのベクトルとの内積である。両ベクトルのなす角度が−π／２〜π／２の範囲内にある場合は、熱量が錐体から流出するためＮＤｊｋが正になる。その角度がこの範囲外にある場合は、熱量が錐体に流入するためＮＤｊｋが負になる。各熱流束に割り当てる熱源の体積を合計したものを節点１の節点領域ＮＤ１とすると、次式が得られる。
【数１２７】
このように導出した節点領域を式（７２）の右辺に代入して、精度向上を図る。
【０１７１】
［整合性、唯一性、正規性について］
整合性とは、節点領域とコントロール・ボリュームとが同体積であることを意味する。整合性が成立すれば、節点レベルでの保存則が２次精度で満足されることになり、精度改善効果が理論上で保証される。
【０１７２】
節点領域の考え方を図１６の長方体要素に適用した場合の整合性について検討するため、式（２５）を次式に書き換える。
【数１２８】
【０１７３】
ここで、α＝ｈｙｈｚ／１８ｈｘ、β＝ｈｘｈｚ／１８ｈｙ、γ＝ｈｘｈｙ／１８ｈｚである。上述の節点領域の算出方法１により次式を得ることができる。
【数１２９】
【０１７４】
同様にＮＤ１ｉ（ｉ＝３，４，・・・，８）を求めて、合計すると、次式が得られ、式（２６）と比較するとコントロール・ボリュームと一致することが分かる。
【数１３０】
【０１７５】
ほかの節点についても、要素形状の対称性を考慮すると整合性が満たされることが分かる。同様に図１８の正三角柱要素においても、次式が成立することが証明できる。式（３０）と比較すると整合性が満たされることが分かる。
【数１３１】
【０１７６】
また、正規性とは、
【数１３２】
すなわち、式（７２）の右辺への追加項の同要素全節点についての合計がゼロになることを意味する。正規性が成立すれば、要素レベルでも、全解析領域でも保存則が満たされることが保証される。長方体要素と正三角柱要素において、正規性が保たれていることが、式（８６）と式（８７）の結果、および要素形状の対称性から分かる。
【０１７７】
また、節点領域の唯一性とは、上記算出方法１と算出方法２で同じ結果が得られることを意味する。その理論証明は、算出方法２におけるバリエーションが多いため煩雑である。
【０１７８】
一般形状の六面体要素、五面体要素においては、ヤコビアンが絡んでくるため、整合性、正規性と唯一性について理論上の議論は困難である。代わりに、整合性については、後述の「作用及び効果」の欄で述べる数値解析において精度改善について検証する。正規性と唯一性についての検討は、以下の形状を持つ要素において行う。立方体領域（−５≦ｘ，ｙ，ｚ≦５）を幅２の立方体のメッシュに分割し、節点１をｘ＝ｙ＝ｚ＝−１に、節点２をｘ＝１、ｙ＝ｚ＝−１に置き、残る六節点（六面体要素の場合）、あるいは四節点（五面体要素の場合）はメッシュにおける全格子点に置いた組み合わせでできた全ての要素形状について検証を行った。これらの要素形状ではアスベクト比が３以下の実用解析で使用される要素形状をカバーできたといえる。この際、節点領域の算出方法２においては、正の係数の並ぶ順序については節点番号１から探し出して、負の係数の並ぶ順序については節点番号１からのと節点番号Ｎからの二通りで行った。その結果、体積がゼロ以下の要素と凹形状の要素を除いて、全テスト要素において正規性と唯一性が満たされていることが確認できた。
【０１７９】
［三次元六面体・五面体における作用及び効果］
保存型離散化手法の精度改善効果を二つの熱伝導問題の解析事例により示す。この際、解析領域を（０≦ｘ，ｙ，ｚ≦１）の立方体領域とし、λ＝１とした。基本境界条件問題としては、式（８８）のソース項分布に式（８９）の境界条件を加えた。その厳密解は式（９０）となる。また、基本・自然境界条件問題としては、式（９１）のソース項分布に式（９２）の基本境界条件と式（９３）の自然境界条件を加えた。その厳密解は式（９４）となる。
【数１３３】
【０１８０】
解析の数値誤差については次式により評価する。
【数１３４】
ここで、ＴＡｉ：節点ｉでの厳密解、ＴＮｉ：節点ｉでの数値解、ＴＡｍａｘ：厳密解の最大値、である。
【０１８１】
数値解析に使用する八パターンの六面体メッシュを図５２に示しており、ｘｙ平面上の四辺形メッシュをｚ方向に押し出して１６層を作成したものである。“●”で節点を、括弧内の数字で総節点数を表わす。メッシュが対称性を持つため、図中にはｘｙ平面上の４分の１の部分のみを示す。図５３（ａ）に基本境界条件問題、図５３（ｂ）に基本・自然境界条件問題での解析誤差を示す。長方体以外のメッシュでは、従来の有限要素法（ＧＦＥ）に比べ、ソース項の保存型離散化（ＣＳＦＥ）の解析精度が大幅に改善されたことが確認できる。また、解析誤差がメッシュ幅の自乗、すなわち、３次元問題では１／Ｎｏｄｅ２／３と比例すると指摘された（非特許文献１）ため、図５４には直線で示す長方体メッシュでの結果を基準にして、各要素形状における節点数と解析誤差との関係を示す。ソース項の保存型離散化（ＣＳＦＥ）は、長方体に近いまで解析精度が改善されたことが確認できる。
【０１８２】
図５５に使用する四パターンの五面体メッシュを示しており、ｘｙ平面上の三角形メッシュをｚ方向に押し出して１６層を作成したものである。メッシュの対称性を考慮して図中にはｘｙ平面上の４分の１のメッシュのみを示す。図５６（ａ）に基本境界条件問題、図５６（ｂ）に基本・自然境界条件問題での解析誤差を示す。特に、自然境界を含む解析条件では、従来の有限要素法（ＧＦＥ）に比べ、ソース項の保存型離散化（ＣＳＦＥ）の解析精度が大幅に改善されたことが確認できる。
【０１８３】
本提案では、解析精度の改善や、メッシュ作成工数の低減、解析時間の短縮などの効果が得られる。市販の解析ソフトに簡単な修正を加えるだけで実現でき、計算コストも増えない。
【０１８４】
［本発明の簡単に応用展開できる解析分野］
熱伝導問題を例に本発明の応用について説明してきたが、勿論、Poissonの方程式が記述する現象の全てにそのまま展開できる。例えば、流れ問題、電・磁気問題、弾性変形問題、物質拡散問題などが応用できる分野として挙げられる。
【０１８５】
本発明の基本的な考え方は、二次精度の保存則の観点から楕円型演算子の離散化結果におけるコントロール・ボリュームを見出し、ソース項の積分領域をそれに合わせたものである。したがって、この考え方は、保存則から導出された楕円型演算子を含む各種の支配方程式の有限要素解析に容易に応用できる。すなわち、物理量の発生・消失と流出・流入を算出するとき、その積分領域について、本発明で導出した節点領域の考え方を導入して積分すればよい。例えば、非定常移流拡散解析において、非定常項と移流項の離散化に対して本発明におけるソース項の積分領域の考え方が適用できる。
【０１８６】
また、荷重、体積力、質量が係わる従来の有限要素法による解析においては、それらの扱い方は、式（２３）の右辺におけるソースと同じように、
【数１３５】
により節点方程式に分配されるため、不整合性問題も存在する。これらの問題は本発明でのソース項の扱い方で改善することもできる。三次元問題の解析においては、荷重、体積力、質量について式（７２）の右辺のソース項と同様に追加項を加えて修正し、対象節点Ｐに寄与する荷重、体積力、質量は次式により算出する。
【数１３６】
ここで、ｆは荷重の任意成分、体積力の任意成分、質量を意味する。
【０１８７】
［本発明のまとめ］
以上説明したように、本発明に係る有限要素法を用いた解析方法にあっては、一般関数項（熱問題等の場合はソース項、荷重、体積力、質量が係る問題の場合はフォース項）を積分する一般関数項積分工程（ソース項積分工程、フォース項積分項手値）（ステップＳ４）にあって、Galerkin有限要素法による二階微分項の離散化結果を基に定義された節点領域の考え方を導入し、要素の代表値を用いた一般関数項（ソース項、フォース項）を積分するので、節点領域の大きさに見合った値（ソース量）を取り入れることができる。これにより、メッシュ・パターンによる解析誤差を低減することができ、品質が低い要素（品質の低いメッシュ・パターン）からも高精度な数値解を得ることが可能となる。このため、メッシュ・パターンの品質を高める作業（メッシュ・パターンを例えば正三角形、正四面体、正方形に近づける作業）を不要とすることができるものでありながら、アルゴリズム（プログラム）の修正による計算量の増加が殆んどないので、コンピュータによる演算処理の増大も防止することができ、総じて解析時間の短縮化を図ることができる。
【０１８８】
また、境界条件については、節点領域の範囲と合致したものを導入する境界条件導入工程（ステップＳ６）を設けたので、節点レベルで精度を向上させることができる。これにより、さらにメッシュ・パターンによる解析誤差を低減することができ、品質が低い要素（品質の低いメッシュ・パターン）からも高精度な数値解を得ることが可能となる。
【０１８９】
なお、以上説明した実施の形態においては、特に熱問題等において、ソース項に、二次元三角形の場合は（ＮＤ−Ｓ／３）ＱＧを追加項とし、三次元四面体の場合は（ＮＤ−Ｖ／４）ＱＧを追加項とし、二次元四角形の場合は
【数１３７】
を追加項とし、三次元六面体・五面体の場合は、
【数１３８】
を追加項としてソース項に追加したものを説明したが、これに限らず、ソース項は要素の幾何学中心でのソース値を用いたソース項に変更する手法であれば、例えばソース項（式（３３）、式（３６）、式（４６）、式（７２）の右辺全部）を単純に「ＮＤ・ＱＧ」に置き換えることも考えられる。
【０１９０】
また、本実施の形態においては、特に荷重、体積力、質量問題等において、一般関数項に、二次元三角形の場合は（ＮＤ−Ｓ／３）ｆＧを追加項とし、三次元四面体の場合は（ＮＤ−Ｖ／４）ｆＧを追加項とし、二次元四角形の場合は
【数１３９】
を追加項とし、三次元六面体・五面体の場合は
【数１４０】
を追加項として一般関数項に追加したものを説明したが、これに限らず、一般関数項は要素の幾何学中心での値を用いた一般関数項に変更する手法であれば、例えば一般関数項を単純に「ＮＤ・ｆＧ」に置き換えることも考えられる。
【０１９１】
また、本実施の形態では、有限要素法を用いた解析方法、及び有限要素法を用いた解析プログラム、として本発明を説明したが、勿論、これら解析方法や解析プログラムを実行する解析装置（解析システム）や、当該プログラムを記録した記録媒体などの形態で本発明を利用するものも、本発明の適用範囲内である。
【産業上の利用可能性】
【０１９２】
本発明に係る有限要素法を用いた解析方法、及び有限要素法を用いた解析プログラムは、解析対象の熱伝導解析、流体の流れ解析、電気流れ解析、磁気流れ解析、弾性変形解析、物質拡散解析、荷重解析、体積力解析、質量解析、等の各種解析に用いることが可能であり、特にメッシュ・パターンの品質が低いものであっても誤差の低減が求められるような有限要素法の解析に用いて好適である。
【符号の説明】
【０１９３】
１ 要素における節点
２ 要素における節点
３ 要素における節点
４ 要素のおける節点
Ｐ 解析領域における節点
ＮＤ 節点領域
Ｑ 一般関数、ソース項
ＱＧ 要素ソース代表値
ＱＯ 要素ソース代表値
Ｗ Galerkin重み関数
ｅ 要素
ｆ 一般関数、荷重、体積力、質量項
ｆＧ 要素荷重、体積力、質量の代表値
ｆＯ 要素荷重、体積力、質量の代表値
Ω １節点周りの要素領域
Ｓ１ 選定工程（ステップＳ１）
Ｓ２ 分割工程（ステップＳ２）
Ｓ３ 要素マトリクス作成工程（ステップＳ３）
Ｓ４ 一般関数項積分工程、ソース項積分工程（ステップＳ４）
Ｓ５ 連立方程式作成工程（ステップＳ５）
Ｓ６ 境界条件導入工程（ステップＳ６）
Ｓ７ 演算工程（ステップＳ７）
Ｓ８ 演算工程（ステップＳ８）
【特許請求の範囲】
【請求項１】
解析対象の解析領域を選定する選定工程と、
前記解析領域を計算対象として複数の要素に分割する分割工程と、
前記複数の要素のうちの、ある要素にある節点についてGalerkin重み関数をかけて要素積分をし、各要素のマトリクスを作成する要素マトリクス作成工程と、
Galerkin重み関数と一般関数との積からなる一般関数項を積分する一般関数項積分工程と、
前記ある節点周りの領域の要素における各要素のマトリクスの和と、前記一般関数項を積分した値の和とに基づき、連立方程式を作成する連立方程式作成工程と、
前記連立方程式に境界条件の導入をする境界条件導入工程と、
前記連立方程式を演算して数値解を得る演算工程と、を備えた有限要素法を用いた解析方法において、
前記一般関数項積分工程にあって、Galerkin有限要素法による二階微分項の離散化結果を基に定義された節点領域の考え方を導入し、要素の代表値を用いた一般関数項を積分する、
ことを特徴とする有限要素法を用いた解析方法。
【請求項２】
前記一般関数項積分工程にあって、節点領域の大きさに合わせて一般関数を節点に分配する、
ことを特徴とする請求項１記載の有限要素法を用いた解析方法。
【請求項３】
前記一般関数項積分工程は、前記一般関数項をソース項としたソース項積分工程であり、
前記ソース項積分工程にあって、要素の代表ソース値を用いたソース項を積分する、
ことを特徴とする請求項１又は２記載の有限要素法を用いた解析方法。
【請求項４】
前記一般関数項積分工程は、前記一般関数項を荷重、体積力、質量としたフォース項積分工程であり、
前記フォース項積分工程にあって、要素の代表フォース値を用いたフォース項を積分する、
ことを特徴とする請求項１又は２記載の有限要素法を用いた解析方法。
【請求項５】
前記複数の要素は、二次元三角形の要素であり、
前記ソース項積分工程にあって、前記各要素における節点領域の面積をＮＤ、前記要素の代表ソース値をＱＧ、前記要素の面積をＳとした場合に、（ＮＤ−Ｓ／３）ＱＧを追加項としてソース項に追加する、
ことを特徴とする請求項３記載の有限要素法を用いた解析方法。
【請求項６】
前記複数の要素は、二次元三角形の要素であり、
前記フォース項積分工程にあって、前記各要素における節点領域の面積をＮＤ、前記要素の代表フォース値をfＧ、前記要素の面積をＳとした場合に、（ＮＤ−Ｓ／３）fＧを追加項としてフォース項に追加する、
ことを特徴とする請求項４記載の有限要素法を用いた解析方法。
【請求項７】
前記複数の要素は、三次元四面体の要素であり、
前記ソース項積分工程にあって、前記各要素における節点領域の体積をＮＤ、前記要素の代表ソース値をＱＧ、前記要素の体積をＶとした場合に、（ＮＤ−Ｖ／４）ＱＧを追加項としてソース項に追加し、かつ前記ＮＤとして以下の数式を用いる、
ことを特徴とする請求項３記載の有限要素法を用いた解析方法。
【数１４１】
【請求項８】
前記複数の要素は、三次元四面体の要素であり、
前記フォース項積分工程にあって、前記各要素における節点領域の体積をＮＤ、前記要素の代表フォース値をfＧ、前記要素の体積をＶとした場合に、（ＮＤ−Ｖ／４）ＱＧを追加項としてフォース項に追加し、かつ前記ＮＤとして以下の数式を用いる、
ことを特徴とする請求項４記載の有限要素法を用いた解析方法。
【数１４２】
【請求項９】
前記複数の要素は、二次元四角形の要素であり、
前記ソース項積分工程にあって、前記各要素における節点領域をＮＤ、前記Galerkin重み関数をＷP、前記要素の代表ソース値をＱＯとした場合に、以下の数式を追加項としてソース項に追加する、
ことを特徴とする請求項３記載の有限要素法を用いた解析方法。
【数１４３】
【請求項１０】
前記複数の要素は、二次元四角形の要素であり、
前記フォース項積分工程にあって、前記各要素における節点領域をＮＤ、前記Galerkin重み関数をＷP、前記要素の代表フォース値をfＯとした場合に、以下の数式を追加項としてフォース項に追加する、
ことを特徴とする請求項４記載の有限要素法を用いた解析方法。
【数１４４】
【請求項１１】
前記複数の要素は、三次元六面体又は三次元五面体の要素であり、
前記ソース項積分工程にあって、前記各要素における節点領域をＮＤ、前記Galerkin重み関数をＷP、前記要素の代表ソース値をＱＯとした場合に、以下の数式を追加項としてソース項に追加する、
ことを特徴とする請求項３記載の有限要素法を用いた解析方法。
【数１４５】
【請求項１２】
前記複数の要素は、三次元六面体又は三次元五面体の要素であり、
前記フォース項積分工程にあって、前記各要素における節点領域をＮＤ、前記Galerkin重み関数をＷP、前記要素の代表フォース値をfＯとした場合に、以下の数式を追加項としてフォース項に追加する、
ことを特徴とする請求項４記載の有限要素法を用いた解析方法。
【数１４６】
【請求項１３】
前記境界条件導入工程にあって、前記各要素の自然境界を通過する熱流束がゼロでない場合に、節点領域の範囲と合致した境界熱流束の節点への配分を行う、
ことを特徴とする請求項２、３、５、７、９、又は１１記載の有限要素法を用いた解析方法。
【請求項１４】
前記境界条件導入工程にあって、前記各要素の境界での荷重、体積力、質量がゼロでない場合に、節点領域の範囲と合致した荷重、体積力、質量の節点への配分を行う、
ことを特徴とする請求項２、４、６、８、１０、又は１２記載の有限要素法を用いた解析方法。
【請求項１５】
有限要素法を用いた解析演算をコンピュータに実行させるための有限要素法を用いた解析演算プログラムであって、
前記コンピュータに、
解析対象の解析領域を選定する選定工程と、
前記解析領域を計算対象として複数の要素に分割する分割工程と、
前記複数の要素のうちの、ある要素にある節点についてGalerkin重み関数をかけて要素積分をし、各要素のマトリクスを作成する要素マトリクス作成工程と、
Galerkin重み関数と一般関数との積からなり、Galerkin有限要素法による二階微分項の離散化結果を基に定義された節点領域の考え方を導入し、要素の代表値を用いた一般関数項を、積分する一般関数項積分工程と、
前記ある節点周りの領域の要素における各要素のマトリクスの和と、前記一般関数項を積分した値の和とに基づき、連立方程式を作成する連立方程式作成工程と、
前記連立方程式に境界条件の導入をする境界条件導入工程と、
前記連立方程式を演算して数値解を得る演算工程と、を実行させる、
ことを特徴とする有限要素法を用いた解析演算プログラム。
【請求項１６】
前記一般関数項積分工程にあって、節点領域の大きさに合わせて一般関数を節点に分配する
ことを特徴とする請求項１５記載の有限要素法を用いた解析演算プログラム。
【請求項１７】
前記一般関数項積分工程は、前記一般関数項をソース項とし、要素の代表ソース値を用いた前記ソース項を積分するソース項積分工程である、
ことを特徴とする請求項１５又は１６記載の有限要素法を用いた解析演算プログラム。
【請求項１８】
前記一般関数項積分工程は、前記一般関数項を荷重、体積力、質量とし、要素の代表フォース値を用いた前記フォース項を積分するフォース項積分工程である、
ことを特徴とする請求項１５又は１６記載の有限要素法を用いた解析演算プログラム。
【請求項１９】
前記複数の要素は、二次元三角形の要素であり、
前記ソース項積分工程にあって、前記各要素における節点領域の面積をＮＤ、前記要素の代表ソース値をＱＧ、前記要素の面積をＳとした場合に、（ＮＤ−Ｓ／３）ＱＧを追加項としてソース項に追加する、
ことを特徴とする請求項１７記載の有限要素法を用いた解析演算プログラム。
【請求項２０】
前記複数の要素は、二次元三角形の要素であり、
前記フォース項積分工程にあって、前記各要素における節点領域の面積をＮＤ、前記要素の代表フォース値をｆＧ、前記要素の面積をＳとした場合に、（ＮＤ−Ｓ／３）ｆＧを追加項としてフォース項に追加する、
ことを特徴とする請求項１８記載の有限要素法を用いた解析演算プログラム。
【請求項２１】
前記複数の要素は、三次元四面体の要素であり、
前記ソース項積分工程にあって、前記各要素における節点領域の体積をＮＤ、前記要素の代表ソース値をＱＧ、前記要素の体積をＶとした場合に、（ＮＤ−Ｖ／４）ＱＧを追加項としてソース項に追加し、かつ前記ＮＤとして以下の数式を用いる、
を特徴とする請求項１７記載の有限要素法を用いた解析演算プログラム。
【数１４７】
【請求項２２】
前記複数の要素は、三次元四面体の要素であり、
前記フォース項積分工程にあって、前記各要素における節点領域の体積をＮＤ、前記要素の代表フォース値をｆＧ、前記要素の体積をＶとした場合に、（ＮＤ−Ｖ／４）ＱＧを追加項としてフォース項に追加し、かつ前記ＮＤとして以下の数式を用いる、
を特徴とする請求項１８記載の有限要素法を用いた解析演算プログラム。
【数１４８】
【請求項２３】
前記複数の要素は、二次元四角形の要素であり、
前記ソース項積分工程にあって、前記各要素における節点領域をＮＤ、前記Galerkin重み関数をＷP、前記要素の代表ソース値をＱＯとした場合に、以下の数式を追加項としてソース項に追加する、
ことを特徴とする請求項１７記載の有限要素法を用いた解析演算プログラム。
【数１４９】
【請求項２４】
前記複数の要素は、二次元四角形の要素であり、
前記フォース項積分工程にあって、前記各要素における節点領域をＮＤ、前記所定の重み関数をＷP、前記要素の代表フォース値をｆＯとした場合に、以下の数式を追加項としてフォース項に追加する、
ことを特徴とする請求項１８記載の有限要素法を用いた解析演算プログラム。
【数１５０】
【請求項２５】
前記複数の要素は、三次元六面体又は三次元五面体の要素であり、
前記ソース項積分工程にあって、前記各要素における節点領域をＮＤ、前記Galerkin重み関数をＷP、前記要素の代表ソース値をＱＯとした場合に、以下の数式を追加項としてソース項に追加する、
ことを特徴とする請求項１７記載の有限要素法を用いた解析演算プログラム。
【数１５１】
【請求項２６】
前記複数の要素は、三次元六面体又は三次元五面体の要素であり、
前記フォース項積分工程にあって、前記各要素における節点領域をＮＤ、前記Galerkin重み関数をＷP、前記要素の代表フォース値をfＯとした場合に、以下の数式を追加項としてフォース項に追加する、
ことを特徴とする請求項１８記載の有限要素法を用いた解析演算プログラム。
【数１５２】
【請求項２７】
前記境界条件導入工程にあって、前記各要素の自然境界を通過する熱流束がゼロでない場合に、節点領域の範囲と合致した境界熱流束の節点への配分を実行させる、
ことを特徴とする請求項１６、１７、１９、２１、２３、又は２５記載の有限要素法を用いた解析演算プログラム。
【請求項２８】
前記境界条件導入工程にあって、前記各要素の境界で荷重、体積力、質量がゼロでない場合に、節点領域の範囲と合致した荷重、体積力、質量の節点への配分を実行させる、
ことを特徴とする請求項１６、１８、２０、２２、２４、又は２６記載の有限要素法を用いた解析演算プログラム。
【請求項１】
解析対象の解析領域を選定する選定工程と、
前記解析領域を計算対象として複数の要素に分割する分割工程と、
前記複数の要素のうちの、ある要素にある節点についてGalerkin重み関数をかけて要素積分をし、各要素のマトリクスを作成する要素マトリクス作成工程と、
Galerkin重み関数と一般関数との積からなる一般関数項を積分する一般関数項積分工程と、
前記ある節点周りの領域の要素における各要素のマトリクスの和と、前記一般関数項を積分した値の和とに基づき、連立方程式を作成する連立方程式作成工程と、
前記連立方程式に境界条件の導入をする境界条件導入工程と、
前記連立方程式を演算して数値解を得る演算工程と、を備えた有限要素法を用いた解析方法において、
前記一般関数項積分工程にあって、Galerkin有限要素法による二階微分項の離散化結果を基に定義された節点領域の考え方を導入し、要素の代表値を用いた一般関数項を積分する、
ことを特徴とする有限要素法を用いた解析方法。
【請求項２】
前記一般関数項積分工程にあって、節点領域の大きさに合わせて一般関数を節点に分配する、
ことを特徴とする請求項１記載の有限要素法を用いた解析方法。
【請求項３】
前記一般関数項積分工程は、前記一般関数項をソース項としたソース項積分工程であり、
前記ソース項積分工程にあって、要素の代表ソース値を用いたソース項を積分する、
ことを特徴とする請求項１又は２記載の有限要素法を用いた解析方法。
【請求項４】
前記一般関数項積分工程は、前記一般関数項を荷重、体積力、質量としたフォース項積分工程であり、
前記フォース項積分工程にあって、要素の代表フォース値を用いたフォース項を積分する、
ことを特徴とする請求項１又は２記載の有限要素法を用いた解析方法。
【請求項５】
前記複数の要素は、二次元三角形の要素であり、
前記ソース項積分工程にあって、前記各要素における節点領域の面積をＮＤ、前記要素の代表ソース値をＱＧ、前記要素の面積をＳとした場合に、（ＮＤ−Ｓ／３）ＱＧを追加項としてソース項に追加する、
ことを特徴とする請求項３記載の有限要素法を用いた解析方法。
【請求項６】
前記複数の要素は、二次元三角形の要素であり、
前記フォース項積分工程にあって、前記各要素における節点領域の面積をＮＤ、前記要素の代表フォース値をfＧ、前記要素の面積をＳとした場合に、（ＮＤ−Ｓ／３）fＧを追加項としてフォース項に追加する、
ことを特徴とする請求項４記載の有限要素法を用いた解析方法。
【請求項７】
前記複数の要素は、三次元四面体の要素であり、
前記ソース項積分工程にあって、前記各要素における節点領域の体積をＮＤ、前記要素の代表ソース値をＱＧ、前記要素の体積をＶとした場合に、（ＮＤ−Ｖ／４）ＱＧを追加項としてソース項に追加し、かつ前記ＮＤとして以下の数式を用いる、
ことを特徴とする請求項３記載の有限要素法を用いた解析方法。
【数１４１】
【請求項８】
前記複数の要素は、三次元四面体の要素であり、
前記フォース項積分工程にあって、前記各要素における節点領域の体積をＮＤ、前記要素の代表フォース値をfＧ、前記要素の体積をＶとした場合に、（ＮＤ−Ｖ／４）ＱＧを追加項としてフォース項に追加し、かつ前記ＮＤとして以下の数式を用いる、
ことを特徴とする請求項４記載の有限要素法を用いた解析方法。
【数１４２】
【請求項９】
前記複数の要素は、二次元四角形の要素であり、
前記ソース項積分工程にあって、前記各要素における節点領域をＮＤ、前記Galerkin重み関数をＷP、前記要素の代表ソース値をＱＯとした場合に、以下の数式を追加項としてソース項に追加する、
ことを特徴とする請求項３記載の有限要素法を用いた解析方法。
【数１４３】
【請求項１０】
前記複数の要素は、二次元四角形の要素であり、
前記フォース項積分工程にあって、前記各要素における節点領域をＮＤ、前記Galerkin重み関数をＷP、前記要素の代表フォース値をfＯとした場合に、以下の数式を追加項としてフォース項に追加する、
ことを特徴とする請求項４記載の有限要素法を用いた解析方法。
【数１４４】
【請求項１１】
前記複数の要素は、三次元六面体又は三次元五面体の要素であり、
前記ソース項積分工程にあって、前記各要素における節点領域をＮＤ、前記Galerkin重み関数をＷP、前記要素の代表ソース値をＱＯとした場合に、以下の数式を追加項としてソース項に追加する、
ことを特徴とする請求項３記載の有限要素法を用いた解析方法。
【数１４５】
【請求項１２】
前記複数の要素は、三次元六面体又は三次元五面体の要素であり、
前記フォース項積分工程にあって、前記各要素における節点領域をＮＤ、前記Galerkin重み関数をＷP、前記要素の代表フォース値をfＯとした場合に、以下の数式を追加項としてフォース項に追加する、
ことを特徴とする請求項４記載の有限要素法を用いた解析方法。
【数１４６】
【請求項１３】
前記境界条件導入工程にあって、前記各要素の自然境界を通過する熱流束がゼロでない場合に、節点領域の範囲と合致した境界熱流束の節点への配分を行う、
ことを特徴とする請求項２、３、５、７、９、又は１１記載の有限要素法を用いた解析方法。
【請求項１４】
前記境界条件導入工程にあって、前記各要素の境界での荷重、体積力、質量がゼロでない場合に、節点領域の範囲と合致した荷重、体積力、質量の節点への配分を行う、
ことを特徴とする請求項２、４、６、８、１０、又は１２記載の有限要素法を用いた解析方法。
【請求項１５】
有限要素法を用いた解析演算をコンピュータに実行させるための有限要素法を用いた解析演算プログラムであって、
前記コンピュータに、
解析対象の解析領域を選定する選定工程と、
前記解析領域を計算対象として複数の要素に分割する分割工程と、
前記複数の要素のうちの、ある要素にある節点についてGalerkin重み関数をかけて要素積分をし、各要素のマトリクスを作成する要素マトリクス作成工程と、
Galerkin重み関数と一般関数との積からなり、Galerkin有限要素法による二階微分項の離散化結果を基に定義された節点領域の考え方を導入し、要素の代表値を用いた一般関数項を、積分する一般関数項積分工程と、
前記ある節点周りの領域の要素における各要素のマトリクスの和と、前記一般関数項を積分した値の和とに基づき、連立方程式を作成する連立方程式作成工程と、
前記連立方程式に境界条件の導入をする境界条件導入工程と、
前記連立方程式を演算して数値解を得る演算工程と、を実行させる、
ことを特徴とする有限要素法を用いた解析演算プログラム。
【請求項１６】
前記一般関数項積分工程にあって、節点領域の大きさに合わせて一般関数を節点に分配する
ことを特徴とする請求項１５記載の有限要素法を用いた解析演算プログラム。
【請求項１７】
前記一般関数項積分工程は、前記一般関数項をソース項とし、要素の代表ソース値を用いた前記ソース項を積分するソース項積分工程である、
ことを特徴とする請求項１５又は１６記載の有限要素法を用いた解析演算プログラム。
【請求項１８】
前記一般関数項積分工程は、前記一般関数項を荷重、体積力、質量とし、要素の代表フォース値を用いた前記フォース項を積分するフォース項積分工程である、
ことを特徴とする請求項１５又は１６記載の有限要素法を用いた解析演算プログラム。
【請求項１９】
前記複数の要素は、二次元三角形の要素であり、
前記ソース項積分工程にあって、前記各要素における節点領域の面積をＮＤ、前記要素の代表ソース値をＱＧ、前記要素の面積をＳとした場合に、（ＮＤ−Ｓ／３）ＱＧを追加項としてソース項に追加する、
ことを特徴とする請求項１７記載の有限要素法を用いた解析演算プログラム。
【請求項２０】
前記複数の要素は、二次元三角形の要素であり、
前記フォース項積分工程にあって、前記各要素における節点領域の面積をＮＤ、前記要素の代表フォース値をｆＧ、前記要素の面積をＳとした場合に、（ＮＤ−Ｓ／３）ｆＧを追加項としてフォース項に追加する、
ことを特徴とする請求項１８記載の有限要素法を用いた解析演算プログラム。
【請求項２１】
前記複数の要素は、三次元四面体の要素であり、
前記ソース項積分工程にあって、前記各要素における節点領域の体積をＮＤ、前記要素の代表ソース値をＱＧ、前記要素の体積をＶとした場合に、（ＮＤ−Ｖ／４）ＱＧを追加項としてソース項に追加し、かつ前記ＮＤとして以下の数式を用いる、
を特徴とする請求項１７記載の有限要素法を用いた解析演算プログラム。
【数１４７】
【請求項２２】
前記複数の要素は、三次元四面体の要素であり、
前記フォース項積分工程にあって、前記各要素における節点領域の体積をＮＤ、前記要素の代表フォース値をｆＧ、前記要素の体積をＶとした場合に、（ＮＤ−Ｖ／４）ＱＧを追加項としてフォース項に追加し、かつ前記ＮＤとして以下の数式を用いる、
を特徴とする請求項１８記載の有限要素法を用いた解析演算プログラム。
【数１４８】
【請求項２３】
前記複数の要素は、二次元四角形の要素であり、
前記ソース項積分工程にあって、前記各要素における節点領域をＮＤ、前記Galerkin重み関数をＷP、前記要素の代表ソース値をＱＯとした場合に、以下の数式を追加項としてソース項に追加する、
ことを特徴とする請求項１７記載の有限要素法を用いた解析演算プログラム。
【数１４９】
【請求項２４】
前記複数の要素は、二次元四角形の要素であり、
前記フォース項積分工程にあって、前記各要素における節点領域をＮＤ、前記所定の重み関数をＷP、前記要素の代表フォース値をｆＯとした場合に、以下の数式を追加項としてフォース項に追加する、
ことを特徴とする請求項１８記載の有限要素法を用いた解析演算プログラム。
【数１５０】
【請求項２５】
前記複数の要素は、三次元六面体又は三次元五面体の要素であり、
前記ソース項積分工程にあって、前記各要素における節点領域をＮＤ、前記Galerkin重み関数をＷP、前記要素の代表ソース値をＱＯとした場合に、以下の数式を追加項としてソース項に追加する、
ことを特徴とする請求項１７記載の有限要素法を用いた解析演算プログラム。
【数１５１】
【請求項２６】
前記複数の要素は、三次元六面体又は三次元五面体の要素であり、
前記フォース項積分工程にあって、前記各要素における節点領域をＮＤ、前記Galerkin重み関数をＷP、前記要素の代表フォース値をfＯとした場合に、以下の数式を追加項としてフォース項に追加する、
ことを特徴とする請求項１８記載の有限要素法を用いた解析演算プログラム。
【数１５２】
【請求項２７】
前記境界条件導入工程にあって、前記各要素の自然境界を通過する熱流束がゼロでない場合に、節点領域の範囲と合致した境界熱流束の節点への配分を実行させる、
ことを特徴とする請求項１６、１７、１９、２１、２３、又は２５記載の有限要素法を用いた解析演算プログラム。
【請求項２８】
前記境界条件導入工程にあって、前記各要素の境界で荷重、体積力、質量がゼロでない場合に、節点領域の範囲と合致した荷重、体積力、質量の節点への配分を実行させる、
ことを特徴とする請求項１６、１８、２０、２２、２４、又は２６記載の有限要素法を用いた解析演算プログラム。
【図１】
【図２】
【図３】
【図４】
【図５】
【図６】
【図７】
【図８】
【図９】
【図１０】
【図１１】
【図１２】
【図１３】
【図１４】
【図１５】
【図１６】
【図１７】
【図１８】
【図１９】
【図２０】
【図２１】
【図２２】
【図２３】
【図２４】
【図２５】
【図２６】
【図２７】
【図２８】
【図２９】
【図３０】
【図３１】
【図３２】
【図３３】
【図３４】
【図３５】
【図３６】
【図３７】
【図３８】
【図３９】
【図４０】
【図４１】
【図４２】
【図４３】
【図４４】
【図４５】
【図４６】
【図４７】
【図４８】
【図４９】
【図５０】
【図５１】
【図５２】
【図５３】
【図５４】
【図５５】
【図５６】
【図５７】
【図２】
【図３】
【図４】
【図５】
【図６】
【図７】
【図８】
【図９】
【図１０】
【図１１】
【図１２】
【図１３】
【図１４】
【図１５】
【図１６】
【図１７】
【図１８】
【図１９】
【図２０】
【図２１】
【図２２】
【図２３】
【図２４】
【図２５】
【図２６】
【図２７】
【図２８】
【図２９】
【図３０】
【図３１】
【図３２】
【図３３】
【図３４】
【図３５】
【図３６】
【図３７】
【図３８】
【図３９】
【図４０】
【図４１】
【図４２】
【図４３】
【図４４】
【図４５】
【図４６】
【図４７】
【図４８】
【図４９】
【図５０】
【図５１】
【図５２】
【図５３】
【図５４】
【図５５】
【図５６】
【図５７】
【公開番号】特開２０１２−７４０００（Ｐ２０１２−７４０００Ａ）
【公開日】平成２４年４月１２日（２０１２．４．１２）
【国際特許分類】
【出願番号】特願２０１１−４７４４５（Ｐ２０１１−４７４４５）
【出願日】平成２３年３月４日（２０１１．３．４）
【出願人】（０００１００７６８）アイシン・エィ・ダブリュ株式会社 (3,717)
【Ｆターム（参考）】
【公開日】平成２４年４月１２日（２０１２．４．１２）
【国際特許分類】
【出願日】平成２３年３月４日（２０１１．３．４）
【出願人】（０００１００７６８）アイシン・エィ・ダブリュ株式会社 (3,717)
【Ｆターム（参考）】
[ Back to top ]