(57)【特許請求の範囲】
【請求項１】
非構造格子の有限体積法または有限要素法による数値解析における種々の係数行列を離散nabla演算子で表わして、前記係数行列に基づく連立一次方程式を解くことにより偏微分方程式の境界値問題の解を求める数値解析装置であって、前記係数行列の各要素を、１次元の線分、２次元の３角形又は３次元の４面体のいずれかである単体要素の形状に依存する係数行列基底と、前記単体要素の形状に依存しない部分係数因子の積和で表わして、前記係数行列基底を前記単体要素の形状に対応させて記憶する係数行列基底記憶手段と、前記部分係数因子を前記単体要素の種類別に計算して生成する部分係数因子生成手段と、行列計算で必要になるたびに前記係数行列基底記憶手段を参照して得た係数行列基底と前記部分係数因子生成手段で生成された部分係数因子とから前記係数行列基底と前記部分係数因子の積和を計算して係数行列要素を求める係数行列要素計算手段と、前記単体要素の種類とその座標が与えられると前記係数行列基底のうちの反変基底ベクトル（この反変基底ベクトルは、他の係数行列基底と異なり、局所座標ξiを位置ベクトルｒで偏微分したもので、離散ナブラ演算子を具体的に計算するためのものである。）を計算するベクトル計算LSIチップとを具備することを特徴とする数値解析装置。
【請求項２】
前記部分係数因子生成手段は、前記単体要素ｅの節点ａに対応する離散nabla演算子∇aのうちの前記単体要素の形状に依存しない部分係数因子を、<∂Ｎa／∂ξi>e (ｉ＝１〜ｋ、ｋは次元数)という平均計算で求める手段を備え、前記係数行列基底記憶手段は、前記単体要素ｅの節点ａに対応する離散nabla演算子∇aのうちの前記単体要素の形状に依存する係数行列基底ｇeiを、ｇei＝<ｇi>e (ｉ＝１〜ｋ、ｋは次元数)という平均計算で求めた結果を記憶する手段を備え、前記係数行列要素計算手段は、前記単体要素ｅの節点ａに対応する離散nabla演算子∇aを、∇a＝Σiｇei<∂Ｎa／∂ξi>e （和は、ｉ＝１〜ｋ、ｋは次元数）（ただし、ａは、前記単体要素の節点番号である。ｉは、座標軸の番号である。Ｎaは、形状関数である。ξiは計算座標である。ｇeiは、ｇei≡<ｇi>e＝<∇ξi>e＝(１／Ωe(ｘ))∫{∇ξidΩ(ｘ)} (積分範囲はΩe)である。ｘは物理空間の座標である。Ωは物理空間ｘの領域である。）として、係数行列基底ｇeiと部分係数因子との積和で求める手段を備えることを特徴とする請求項１記載の数値解析装置。
【請求項３】
前記部分係数因子生成手段は、前記単体要素ｅの節点ａ、ｂに対応する拡散行列Ｄabのうちの前記単体要素の形状に依存しない部分係数因子を、<(∂Ｎa／∂ξi)(∂Ｎb／∂ξj)>(ｉ、ｊ＝１〜ｋ、ｋは次元数)という平均計算で求める手段を備え、前記係数行列基底記憶手段は、前記単体要素ｅの節点ａ、ｂに対応する拡散行列Ｄabのうちの前記単体要素の形状に依存する係数行列基底を、Ωe(ｘ)ｇei・ｇej(ｉ、ｊ＝１〜ｋ、ｋは次元数)という計算で求めた結果を記憶する手段を備え、前記係数行列要素計算手段は、前記単体要素ｅの節点ａ、ｂに対応する拡散行列Ｄabを、Ｄab＝ΣiΣjΩe(ｘ)ｇei・ｇej<(∂Ｎa／∂ξi)(∂Ｎb／∂ξj)>（和は、ｉ、ｊ＝１〜ｋ、ｋは次元数）（ただし、ａ、ｂは、前記単体要素の節点番号である。ｉ、ｊは、座標軸の番号である。Ｎa、Ｎbは、形状関数である。ξi、ξjは計算座標である。ｇeiは、ｇei≡<ｇi>e＝<∇ξi>e＝(１／Ωe(ｘ))∫{∇ξidΩ(ｘ)} (積分範囲はΩe)である。ｇejも同様である。ｘは物理空間の座標である。Ωは物理空間ｘの領域である。）として、係数行列基底と部分係数因子との積和で求める手段を備えることを特徴とする請求項１記載の数値解析装置。
【請求項４】
前記部分係数因子生成手段は、前記単体要素ｅの節点ａ、ｂに対応するBTD行列Ｂabのうちの前記単体要素の形状に依存しない部分係数因子を、<(∂Ｎa／∂ξi)(∂Ｎb／∂ξj)>(ｉ、ｊ＝１〜ｋ、ｋは次元数)という平均計算で求める手段を備え、前記係数行列基底記憶手段は、前記単体要素ｅの節点ａ、ｂに対応するBTD行列Ｂabのうちの前記単体要素の形状に依存する係数行列基底を、Ωe(ｘ)ｖe・ｇeiｖe・ｇej(ｉ、ｊ＝１〜ｋ、ｋは次元数である。ｖeは要素ｅの体積である。)という計算で求めた結果を記憶する手段を備え、前記係数行列要素計算手段は、前記単体要素ｅの節点ａ、ｂに対応するBTD行列Ｂabを、Ｂab＝ΣiΣjΩe(ｘ)ｖe・ｇeiｖe・ｇej<(∂Ｎa／∂ξi)(∂Ｎb／∂ξj)>（和は、ｉ、ｊ＝１〜ｋ、ｋは次元数）（ただし、ａ、ｂは、前記単体要素の節点番号である。ｉ、ｊは、座標軸の番号である。Ｎa、Ｎbは、形状関数である。ξi、ξjは計算座標である。ｇeiは、ｇei≡<ｇi>e＝<∇ξi>e＝(１／Ωe(ｘ))∫{∇ξidΩ(ｘ)} (積分範囲はΩe)である。ｇejも同様である。ｘは物理空間の座標である。Ωは物理空間ｘの領域である。）として、係数行列基底と部分係数因子との積和で求める手段を備えることを特徴とする請求項１記載の数値解析装置。
【発明の詳細な説明】
【０００１】
【発明の属する技術分野】
本発明は、数値解析装置に関し、特に、離散nabla演算子と非記憶型係数行列を利用して、非構造格子の有限体積法または有限要素法により計算する数値解析装置に関する。
【０００２】
【従来の技術】
流体力学や構造力学などにおける問題に現れる偏微分方程式の境界値問題を、計算機による数値計算で解く方法の１つに有限要素法がある。有限要素法について、簡単に説明する。詳細については、文献１［棚橋隆彦著「流れの有限要素法解析Ｉ，II」（1997）、朝倉書店発行］などを参照されたい。
【０００３】
場の関数φの偏微分方程式
Ｌ(φ)＝ｆ  （領域Ω上で）
を、境界条件
φ＝φ0  （境界Γ1上で）
(∂φ／∂ｎ)＝ｑ  （境界Γ2上で）
で解く場合を例にして説明する。ただし、ｎは境界の法線ベクトルである。境界Γ1は、関数φの値φ0が与えられた境界（Dirichlet境界）である。境界Γ2は、関数φの微分値ｑが与えられた境界（Neumann境界）である。具体的な偏微分方程式は、例えば、
Ｌ(φ)＝∇2φ＝ｆ
である。
【０００４】
この偏微分方程式の近似解をφ*とする。近似解は、適当な基底関数の一次結合で表わす。厳密解における偏微分方程式の左辺の値との残差
ｒ＝Ｌ(φ*)−Ｌ(φ)＝Ｌ(φ*)−ｆ
に重み関数ｗj(ｊ＝１,２,・・・,ｎ；ｎは重み関数の数すなわち未定係数の数)をかけて積分したものが０となるように、
(ｗj,ｒ)≡∫ｗjｒdΩ＝０  （積分範囲は領域Ω）
として、近似解φ*を決定する。任意の重み関数に対して、この積分方程式を満足するｒは０であるから、そのとき近似解φ*は厳密解に一致する。
【０００５】
同様に、境界条件の残差についても、重み関数をかけて積分すると、
∫{Ｌ(φ*)−ｆ}ｗjdΩ  （積分範囲は領域Ω）
＋∫{φ*−φ0}ｗjdΓ1  （積分範囲は境界Γ1）
−∫{(∂φ*／∂ｎ)−ｑ}ｗjdΓ2  （積分範囲は境界Γ2）
＝０
となる。この場合も、任意の重み関数に対して、この積分方程式を満足すれば、近似解φ*は厳密解に一致する。
【０００６】
ここで、積分方程式が簡単になるように、重み関数に条件を課して、近似解φ*の制限を緩める。重み関数ｗjに、境界Γ1上で０となるという条件を課すと、
∫{Ｌ(φ*)−ｆ}ｗjdΩ  （積分範囲は領域Ω）
＝∫{(∂φ*／∂ｎ)−ｑ}ｗjdΓ  （積分範囲は境界Γ2）
となる。これを、有限要素法の基礎式とする。
【０００７】
この基礎式を部分積分すると、弱形式が得られる。例えば、
Ｌ(φ)＝∇2φ＝ｆ
の場合、
∫{∇2φ*−ｆ}ｗjdΩ  （積分範囲は領域Ω）
＝∫{(∂φ*／∂ｎ)−ｑ}ｗjdΓ  （積分範囲は境界Γ2）
であるから、Greenの第１公式を使って、
∫(∇φ*・∇ｗj)dΩ  （積分範囲は領域Ω）
＋∫∇2φ*ｗjdΩ  （積分範囲は領域Ω）
＝∫ｗj(∂φ*／∂ｎ)dΓ  （積分範囲は境界Γ2）
とすると、弱形式は、
∫{(∇φ*・∇ｗj)＋ｆｗj}dΩ  （積分範囲は領域Ω）
＝∫ｑｗjdΓ  （積分範囲は境界Γ2）
となる。変分原理により、汎関数Ｊの変分が０となる条件から、弱形式を求めることもできる。重み関数の選び方により、種々の解法がある。例えば、有限要素法のGalerkin法では、重み関数を形状関数とし、写像関数も同じ関数とする。形状関数と写像関数については、後で説明する。
【０００８】
例として、Galerkin法で、
∇2φ＝ｆ
を解く方法を説明する。領域Ωを、ｎ個の節点を頂点とする４面体の要素に分割する。節点ｉの形状関数をＮiとする。節点ｉにおける近似解φ*の値をφi（未知数）とする。節点ｊにおける重み関数ｗjの値をｗj*とすると、
φ*＝ΣiＮiφi＋ΣiＮiｈi  （和は、ｉ＝１〜ｎ）
ｗj＝Ｎjｗj*
となる。ただし、ｈiは、境界Γ1における既知の値であり、領域内部では０である。
【０００９】
これを弱形式に代入すると、
∫{(∇(ΣiＮiφi＋ΣiＮiｈi)・∇(Ｎjｗj*))  （和は、ｉ＝１〜ｎ）
＋ｆ(Ｎjｗj*)}dΩ  （積分範囲は領域Ω）
＝∫ｑ(Ｎjｗj*)dΓ  （積分範囲は境界Γ2）
となる。ｗj*は任意であるから、その係数は０である。したがって、
∫{(∇(ΣiＮiφi＋ΣiＮiｈi)・∇Ｎj)＋ｆＮj}dΩ
（積分範囲は領域Ω、和は、ｉ＝１〜ｎ）
＝∫(ｑＮj)dΓ  （積分範囲は境界Γ2）
となる。整理すると、
Σiφi∫(∇Ｎi・∇Ｎj)dΩ  （積分範囲は領域Ω、和は、ｉ＝１〜ｎ）
＝∫(ｑＮj)dΓ  （積分範囲は境界Γ2）
−∫(ｆＮj)dΩ  （積分範囲は領域Ω）
−Σiｈi∫(∇Ｎi・∇Ｎj)dΩ  （積分範囲は領域Ω、和は、ｉ＝１〜ｎ）
となる。
【００１０】
ここで、
Ａij＝∫(∇Ｎi・∇Ｎj)dΩ  （積分範囲は領域Ω）
Ｂj＝∫(ｑＮj)dΓ  （積分範囲は境界Γ2）
−∫(ｆＮj)dΩ  （積分範囲は領域Ω）
−Σiｈi∫(∇Ｎi・∇Ｎj)dΩ  （積分範囲は領域Ω、和は、ｉ＝１〜ｎ）
とすると、
ΣiＡijφi＝Ｂj  （和は、ｉ＝１〜ｎ）
となる。すなわち、
(φi)＝(Ａij)-1(Ｂj)
となる。ただし、(φi)は、φiを第ｉ成分とする縦ベクトル、(Ｂj)は、Ｂjを第ｊ成分とする縦ベクトル、(Ａij)-1は、Ａijを第ｉ,ｊ成分とする行列の逆行列である。このようにして、連立一次方程式を解くことによって、偏微分方程式の境界値問題を解くことができる。
【００１１】
この場合、係数行列要素
Ａij＝∫(∇Ｎi・∇Ｎj)dΩ  （積分範囲は領域Ω）
の計算が重要になる。これを簡単に計算するために、物理空間を計算空間に写像して、nabla演算子∇を離散化する。離散nabla演算子で係数行列要素を記述し、連立一次方程式を解くことによって、偏微分方程式の境界値問題を解く。係数行列要素のサイズは、（節点数）×（要素数）であるので、精密な計算を行なうためには、莫大な記憶容量を必要とし、計算時間も長くなる。
【００１２】
ところで、非線形力学系のシステム方程式は、時間発展方程式として記述できる。上記の例では、時間変数を含まない偏微分方程式について説明したが、有限要素法は、時間変数を含む偏微分方程式にも適用できる。時間発展方程式の中には、２種類の微分演算子、∂／∂tと∇が含まれる。∂／∂tの時間積分は、いろいろな時間進行法を産みだす。一方、nabla演算子∇の空間積分は、いろいろな離散化手法を産みだす。
【００１３】
非線形システム方程式のコンピュータ・シミュレーションにおいて、この時間進行法と空間の離散化が、複雑系大規模問題の解決において大きな課題である。特に、非構造格子を用いた有限体積法または有限要素法において、空間の離散化手続きが一般に複雑である。
【００１４】
特開平10-177591号公報に開示されている「数値解析装置及び数値解析方法並びに集積回路」において、本発明者らは、離散nabla演算子を利用して、簡易な構成でかつ短時間で数値解析できる数値解析装置を提案した。この数値解析装置では、解析対象の場を分割して得られる各要素の要素情報に基づいて、空間微分のnabla演算子を離散化した離散nabla演算子を算出する離散nabla演算子算出手段を設ける。算出した離散nabla演算子を使用して、解析対象の場の連続の式を満足するように、物理量を反復修正することにより、当該連続の式を満足する物理量を、物理量算出手段で算出する。
【００１５】
【発明が解決しようとする課題】
しかし、従来の有限要素法においては、離散nabla演算子を使用して計算しても、そのまま係数行列を計算すると長い時間がかかるという問題があった。また、一度計算した係数行列を記憶すると、大規模問題においては莫大な記憶容量を必要とするという問題があった。
【００１６】
本発明は、上記従来の問題を解決して、離散nabla演算子を使用する有限要素法の計算において、記憶容量を少なくしながらも、計算を高速化することを目的とする。
【００１７】
【課題を解決するための手段】
上記の課題を解決するために、本発明では、非構造格子の有限体積法または有限要素法により偏微分方程式の境界値問題の解を計算する数値解析装置を、離散nabla演算子で表わした係数行列要素のうちの要素形状に依存する因子を要素形状に対応させて記憶する係数行列基底記憶手段と、離散nabla演算子で表わした係数行列要素のうちの要素形状に依存しない因子を種類別に非記憶化し部分係数因子を生成する部分係数因子生成手段と、行列計算で必要になるたびに係数行列基底記憶手段を参照して得た係数行列基底と部分係数因子生成手段で生成された部分係数因子とから係数行列要素を生成する計算手段とを具備する構成とした。このように構成したことにより、記憶容量を少なくしながらも、計算を高速化できる。
【００１８】
すなわち、nabla演算子∇がベクトル量であり、座標系に依存しない量であることに着目し、離散nabla演算子を定義する。この基本アイデアに基づいて、離散nabla演算子で表現した係数行列要素を、２つの部分に分解する。その１つは、要素形状に依存する部分、他は要素形状に依存しない部分である。要素形状に依存する部分のみを記憶し、他の部分は解析的に計算しておく。計算機内部で必要になるたびに、２つの部分から係数行列要素を生成し、非記憶化する。非構造格子の有限体積法または有限要素法において、この離散nabla演算子を用いると、以下の点においてすぐれている。
【００１９】
（１）離散nabla演算子を用いると、任意の要素に対する空間の離散化が、極めて単純にできる。
（２）離散nabla演算子は、非構造格子の有限要素法および有限体積法に適している。
（３）離散nabla演算子は、object化することが極めて容易である。
（４）離散nabla演算子を用いて、非記憶型の係数行列を作ることができる。よって、低容量化と高速化に寄与する。
（５）離散nabla演算子を用いると、いままでにない次世代向けの新しいシミュレーションコードの開発が容易になる。
（６）離散nabla演算子をLSI化した専用計算機により、複雑系の大規模計算が低コストかつ高精度にできるようになる。
（７）離散nabla演算子を用いると、Poisson Solverを高速化できる。
【００２０】
【発明の実施の形態】
以下、本発明の実施の形態について、図１〜図７を参照しながら詳細に説明する。
【００２１】
（実施の形態）
本発明の実施の形態は、非構造格子の有限体積法または有限要素法により計算する場合に、離散nabla演算子のうちの要素形状に依存する因子を記憶し、離散nabla演算子のうちの要素形状に依存しない因子を解析的に計算しておき、必要になるたびに係数行列要素を生成し、係数行列要素をすべて記憶することなく係数行列を解く数値解析装置である。
【００２２】
図１は、本発明の実施の形態における数値解析装置で利用する１次元単体要素のこう配・発散・回転の説明図である。図２は、数値解析装置で利用する２次元単体要素のこう配・発散・回転の説明図である。図３は、数値解析装置で利用する３次元単体要素のこう配・発散・回転の説明図である。
【００２３】
図４は、本発明の実施の形態における数値解析装置で利用する１次元の対称要素のこう配・発散・回転の説明図である。図５は、数値解析装置で利用する２次元の対称要素のこう配・発散・回転の説明図である。図６は、数値解析装置で利用する３次元の対称要素のこう配・発散・回転の説明図である。
【００２４】
図７は、本発明の実施の形態における数値解析装置のデルチップの概念図である。図７において、要素の節点座標生成手段71は、計算対象の要素の節点座標を出力する手段である。デルチップ72は、離散nabla演算子により要素体積と反変基底ベクトルを計算するLSI回路である。データ利用手段73は、出力データｇeiとΩeを用いる計算手段である。
【００２５】
図８は、本発明の実施の形態における数値解析装置の処理手順を示す流れ図である。図９は、本発明の実施の形態における数値解析装置の機能ブロック図である。図９において、入力装置１は、パラメータなどのデータを入力するキーボードなどである。パラメータ境界条件記憶装置２は、偏微分方程式のパラメータと境界条件のデータを格納しておくメモリーである。節点位置ベクトル記憶装置３は、すべての節点の座標値を格納しておくメモリーである。既知ベクトル演算手段４は、パラメータと境界条件から既知ベクトルを計算する手段である。既知ベクトル記憶装置５は、既知ベクトルを格納しておくメモリーである。偏微分方程式情報記憶装置６は、計算可能な偏微分方程式に関する情報を格納しておくメモリーである。表示装置７は、解析結果を表示するＣＲＴディスプレイなどである。係数行列基底演算手段８は、節点座標から反変基底ベクトルと要素体積などを計算する手段である。係数行列基底記憶装置９は、係数行列基底を格納しておくメモリーである。形状関数記憶装置10は、形状関数に関するデータを格納しておくメモリーである。部分係数因子生成手段11は、解析的に計算した形状関数のデータから部分係数因子を生成する手段である。係数行列要素演算手段12は、係数行列基底と部分係数因子から係数行列要素を計算する手段である。行列演算手段13は、係数行列要素と既知ベクトルから関数ベクトルを求める手段である。行列要素記憶ワークエリア14は、係数行列要素を一時的に格納するメモリーである。
【００２６】
上記のように構成された本発明の実施の形態における数値解析装置の動作を説明する。最初に、離散nabla演算子の定義を説明する。連続空間で定義されるnabla演算子
∇≡ｉ（∂／∂ｘ）＋ｊ（∂／∂ｙ）＋ｋ（∂／∂ｚ）
を、離散空間でどのように定義したらよいかを考える。ここで、ｉ，ｊ，ｋは、それぞれｘ，ｙ，ｚの軸方向の単位ベクトルである。
【００２７】
φを任意のテンソル量とする。スカラーは０階のテンソル、ベクトルは１階のテンソル、応力テンソルは２階のテンソルである。例えば、φは、圧力ｐ（スカラー即ち０階のテンソル）または、速度ｖ（ベクトル即ち１階のテンソル）または、応力テンソルＴ（テンソル即ち２階のテンソル）を代表している。
【００２８】
また、star（＊）積は、こう配、発散、回転を代表している。例えば、∇ｐ（こう配）、∇・ｖ（発散（dot積））、∇×ｖ（回転（cross積））、∇ｖ（テンソル積（open積））、∇・Ｔ（テンソルの発散）、∇×Ｔ（テンソルの回転）、∇Ｔ（テンソルのこう配）である。これらの演算を、すべて∇＊φで代表する。
【００２９】
∇＊φの有限要素Ωe(ｒ)の要素平均値を、
<∇＊φ>e≡(１／Ωe(ｒ))∫{∇＊φdΩ(ｒ)}  (積分範囲はΩe)
で定義する。ここで、Ωe(ｒ)は、３次元の場合、４面体／５面体／６面体要素の体積を意味し、２次元の場合は、３角形／４角形要素の面積を意味し、１次元の場合は、線要素の長さを意味する。Ωeは任意の閉領域であり、その要素を十分小さくすれば、当然、
∇＊φ＝lim<∇＊φ>e  (極限はΩe→０)
が成立する。
【００３０】
要素Ωeにおけるφは、要素Ωeの節点番号をａで表わすと、
φ＝ΣaＮaφa  （和は、ａ＝１〜ｍ、ｍは要素の節点数）
で近似できる。節点ａは、隣接する要素の節点でもあるので、形状関数に要素番号も付して区別すべきであるが、煩雑であるし紛れるおそれもないので、単にＮaと書くことにする。ここで、Ｎaは形状関数である。紛れるおそれがない場合は、記号Σを省略して、
φ＝Ｎaφa
と書くこともある。この場合、添字ａは和の規約に従うという。このφを要素平均値の定義式に代入すると、
<∇＊φ>e＝Σa(１／Ωe(ｒ))∫{∇ＮadΩ(ｒ)}＊φa  (積分範囲はΩe)
＝Σa∇a＊φa（和は、ａ＝１〜ｍ）
を得る。
【００３１】
ここで、領域Ωe(ｒ)の節点ａに対応する離散nabla演算子∇aを、
∇a≡(１／Ωe(ｒ))∫{∇ＮadΩ(ｒ)}  (積分範囲はΩe)
で定義する。すなわち、離散nabla演算子は、形状関数Ｎaのこう配の要素平均値である。この定義式が基本アイデアである。そして、この定義式を用いると、要素平均値を用いて、
<∇ｐ>e＝Σa∇aｐa,  <∇・ｖ>e＝Σa∇a・ｖa,
<∇×ｖ>e＝Σa∇a×ｖa,  <∇ｖ>e＝∇aｖa,
<∇・Ｔ>e＝Σa∇a・Ｔa,  <∇×Ｔ>e＝Σa∇a×Ｔa,
<∇Ｔ>e＝Σa∇aＴa
等が直ちに計算できる。
【００３２】
同様にして、粘性応力の要素平均値は、
τe＝<μ[∇ｖ＋(∇ｖ)T]>e
＝μeΣa(∇aｖa＋ｖa∇a)
と表示できる。ここで、μは粘性係数である。μeは、μの要素平均値である。残された問題は、離散nabla演算子の計算法である。このように離散nabla演算子を定義することにより、こう配・発散・回転を、離散nabla演算子を用いて記述できる。
【００３３】
第２に、離散nabla演算子の計算法を説明する。離散nabla演算子の計算法には、体積分で計算する体計算法と、面積分で計算する面計算法がある。物理空間で歪んだ要素は、計算空間に写像することで正規化することができる。よって、離散nabla演算子の計算は、計算空間で実行するのがよい。例えば、ｒ＝(ｘ,ｙ,ｚ)＝(ｘ1,ｘ2,ｘ3)空間を、ξ＝(ξ,η,ζ)＝(ξ1,ξ2,ξ3)空間に写像する。写像は、形状関数Ｎaを用いて、
ｒ(ξ)＝ΣaＮa(ξ)ｒa  （和は、ａ＝１〜ｍ）
と実行される。ｒaは、節点ａの位置ベクトルである。
【００３４】
nabla演算子∇はベクトルであり、座標系に依存しないことを利用して、計算空間のnabla演算子を求める。すなわち、nabla演算子は、物理空間と計算空間で、３次元の場合、
∇＝∇ｘ1(∂／∂ｘ1)＋∇ｘ2(∂／∂ｘ2)＋∇ｘ3(∂／∂ｘ3)  物理空間
∇＝∇ξ1(∂／∂ξ1)＋∇ξ2(∂／∂ξ2)＋∇ξ3(∂／∂ξ3)  計算空間
と表示でき、両者は共に等しい。ここで、反変基底ベクトル
ｇi≡∇ξi  (ｉ＝１，２，３)
を定義すると、nabla演算子は、簡潔に、
∇＝Σiｇi(∂／∂ξi)
＝ｇ1(∂／∂ξ1)＋ｇ2(∂／∂ξ2)＋ｇ3(∂／∂ξ3)
と記述できる。
【００３５】
変数変換のJacobianをＪで表わせば、
ｄΩ(ｒ)＝ＪdΩ(ξ)
Ｊ＝{∂(ｘ，ｙ，ｚ)／∂(ξ，η，ζ)}
が成立する。ここで、要素平均値を、
Ｊe≡<Ｊ>e≡Ωe(ｒ)／Ωe(ξ)
ｇei≡<ｇi>e
で定義すれば、離散nabla演算子は、
∇a≡(１／Ωe(ｒ))∫{∇ＮadΩ(ｒ)}  (積分範囲はΩe)
＝(１／Ωe(ｒ))∫{∇ＮaＪdΩ(ξ)}
≒Σi(Ｊe／Ωe(ｒ))∫{ｇi(∂Ｎa／∂ξi)dΩ(ξ)}
≒Σiｇei(１／Ωe(ξ))∫{(∂Ｎa／∂ξi)dΩ(ξ)}
＝Σiｇei<∂Ｎa／∂ξi>e  （和は、ｉ＝１〜３）
と近似できる。したがって、
∇a＝Σiｇei<∂Ｎa／∂ξi>e  （体計算法）
と記述できる。これを、離散nabla演算子の体計算法と呼ぶことにする。すなわち、離散nabla演算子は、要素の形状に依存する反変基底ベクトルの要素平均値と、要素の形状に依存しない形状関数の一階微分の積和で表わすことができる。
【００３６】
もう一つの計算法が可能である。（∇＊φ）の要素平均は、Greenの定理を用いて、境界積分で表わすこともできる。すなわち、
<∇＊φ>e＝(１／Ωe(ｒ))∫{ｎ＊φdΓ(ｒ)}  (積分範囲はΓe)
である。ここで、Γe＝∂Ωeは、Ωeの境界で、ｎは、面の単位法線ベクトルである。これを、各面の総和で記述すると、
<∇＊φ>e≒(１／Ωe(ｒ))ΣsΔΓs(ｒ)＊φs
＝ΣsΔΓs(ｒ)＊φs／Ωe(ｒ)  （面計算法）
となる。△Γsは、微小要素の面積ベクトルである。Σsは、境界要素∂Ωeの分割要素の総和を意味する。ｓについての総和の意味で、Σを省略してもよい。φsは、φの面平均値である。
【００３７】
このように、離散nabla演算子は、体計算法と面計算法によって計算できる。体計算法は表示が単純で、計算機でシステム的に計算するのに便利である。面計算法は、離散nabla演算子の物理的解釈に適している。表１、表２、表３に、以上の結果をまとめておく。表２から明らかなように、離散nabla演算子は、要素の節点の位置ベクトルの差△ｒが求まれば計算できる。そして、離散nabla演算子は、非構造格子の有限体積法または有限要素法において、最も基本的な役割を演ずる。さらに、有限体積法と有限要素法の両者のもつ長所を混合する概念で、次世代の新しい計算手法を生み出す可能性を有している。
【００３８】
［表１］  離散nabla演算子
<∇＊φ>e＝(１／Ωe(ｒ)){∫∇ＮadΩ(ｒ)}＊φa  (積分範囲はΩe)
＝ΣaΣiｇei<∂Ｎa／∂ξi>e＊φa
<∇＊φ>e＝(１／Ωe(ｒ)){∫φdΓ}＊φ  (積分範囲はΓe)
＝ΣsΔΓs(ｒ)＊φs／Ωe(ｒ)
<∇＊φ>e＝Σa∇a＊φa
∇a≡(１／Ωe(ｒ)){∫∇ＮadΩ(ｒ)}
＝Σiｇei<∂Ｎa／∂ξi>e  （定義式）
【００３９】
［表２］  離散nabla演算子の表示式（６面体要素）
∇a≡Σiｇei<∂Ｎa／∂ξi>e  （和は、ｉ＝１〜３）（定義式）
∇a＝(１／Ωe(ｒ))(Ａ*1ξa＋Ａ*2ηa＋Ａ*3ζa)
∇a＝(１／４Ωe(ｒ))(Ｓ*1ξa＋Ｓ*2ηa＋Ｓ*3ζa)
ただし、
ξ1＝ξ, ξ2＝η, ξ3＝ζ
ｇe1≡<ｇ1>e＝<ｇ2×ｇ3>e／Ｊe  （反変基底ベクトル１）
ｇe2≡<ｇ2>e＝<ｇ3×ｇ1>e／Ｊe  （反変基底ベクトル２）
ｇe3≡<ｇ3>e＝<ｇ1×ｇ2>e／Ｊe  （反変基底ベクトル３）
ｇi≡∂ｒ／∂ξi  （ｉ＝１，２，３）（共変基底ベクトル）
Ａ*1≡<ｇ2×ｇ3>e＝<(∂ｒ／∂ξ2)×(∂ｒ／∂ξ3)>e
Ａ*2≡<ｇ3×ｇ1>e＝<(∂ｒ／∂ξ3)×(∂ｒ／∂ξ1)>e
Ａ*3≡<ｇ1×ｇ2>e＝<(∂ｒ／∂ξ1)×(∂ｒ／∂ξ2)>e
Ｓ*1≡(Δｒ)*2×(Δｒ)*3
Ｓ*2≡(Δｒ)*3×(Δｒ)*1
Ｓ*3≡(Δｒ)*1×(Δｒ)*2
（６面体を、ξ,η,ζがそれぞれ±１となる８点を頂点とする立方体に正規化した場合、ξaは、節点ａのξ座標である。ηaは、節点ａのη座標である。ζaは、節点ａのζ座標である。）
(Δｒ)*1＝(1/4){(Δｒ)21＋(Δｒ)34＋(Δｒ)65＋(Δｒ)78}
(Δｒ)*2＝(1/4){(Δｒ)32＋(Δｒ)41＋(Δｒ)76＋(Δｒ)85}
(Δｒ)*3＝(1/4){(Δｒ)51＋(Δｒ)62＋(Δｒ)73＋(Δｒ)84}
(Δｒ)ab＝ｒa−ｒb  （ａ，ｂは、節点番号）
【００４０】
［表３］  離散nabla演算子の表示式
∇a≡(１／Ωe(ｒ))∫{∇ＮadΩ(ｒ)}
＝Σiｇei<∂Ｎa／∂ξi>e
対称要素の場合
∇a＝Σi(１／Ωe(ｒ))ｇeiξai
∇a＊φa＝ΣiＳ*i＊(δφ)*i／Ωe(ｒ)
単体要素の場合
∇a≡(１／Ωe(ξ))ｇa
∇a＊φa＝(１／Ωe(ｒ))ｇa＊φa
ただし、ｉは、座標軸の番号、ａは、要素の節点番号、
ｇei＝<ｇi>e  （ｉ＝１，２，３）（反変基底ベクトルの要素平均）
【００４１】
第３に、単体要素のこう配・発散・回転について説明する。１次元の線分、２次元の３角形、３次元の４面体を総称して、単体要素と呼ぶ。単体要素に対して、こう配・発散・回転を計算する。反変基底ベクトル
gi＝∇ξi  （ｉ＝１，２，３）
は、単体要素に対して一定であるから、要素平均値ｇeiは、ｇiに等しい。この結果、離散nabla演算子法の体計算法と面計算法は一致し、離散nabla演算子は、近似式でなく解析表示と一致する。そして、各次元に応じて、
Σiξi＝１，  Σiｇi＝０
が成り立つ。
【００４２】
また、
Ｎa＝ξa
であるから、
<∂Ｎa／∂ξi>e＝δia
＝１（ｉ＝ａのとき）
＝０（ｉ≠ａのとき）
に等しい。このとき、
∇a＝ｇa
となる。よって、一般式の離散化は、
∇a＊φa＝ｇa＊φa
となる。これを具体的に示すと、１次元の場合、図１のようになり、２次元の場合、図２のようになり、３次元の場合、図３のようになる。
【００４３】
図１を参照して、１次元の場合を説明する。図１の▲１▼は、節点１であり、▲２▼は、節点２である。▲１▼から▲２▼までの長さＬeの線分を要素とする。▲１▼のｘ座標をｘ1とし、▲２▼のｘ座標をｘ2とする。したがって、
ｘ2−ｘ1＝Ｌe
Ｌe＝｜(1,1)T (ｘ1,ｘ2)T｜
となる。ｘ座標の単位ベクトルをｉとする。ｉは、右向きの単位ベクトルである。ξ1座標の原点を▲２▼とし、左向きを正の方向とし、長さＬeで正規化する。ξ2座標の原点を▲１▼とし、右向きを正の方向とし、長さＬeで正規化する。ξ1座標とξ2座標をｘ座標で表わすと、
ξ1＝(ｘ2−ｘ)／Ｌe
ξ2＝(ｘ−ｘ1)／Ｌe
となる。したがって、
ξ1＋ξ2＝１
となる。反変基底ベクトルｇ1とｇ2は、
ｇ1＝∇ξ1＝−ｉ／Ｌe
ｇ2＝∇ξ2＝ｉ／Ｌe
となる。したがって、
ｇ1＋ｇ2＝０
である。形状関数Ｎ1は、▲１▼で１になり、▲２▼で０になる線形関数であるから、
Ｎ1＝(ｘ2−ｘ)／Ｌe＝ξ1
となる。形状関数Ｎ2は、▲２▼で１になり、▲１▼で０になる線形関数であるから、
Ｎ2＝(ｘ−ｘ1)／Ｌe＝ξ2
となる。したがって、
Ｎa＝ξa  （ａ＝１，２）
と書ける。
【００４４】
図２を参照して、２次元の場合を説明する。▲１▼,▲２▼,▲３▼は、節点番号である。３角形▲１▼▲２▼▲３▼を要素とする。節点▲１▼の(ｘ,ｙ)座標を、（ｘ1,ｙ1）とする。節点▲２▼の(ｘ,ｙ)座標を、（ｘ2,ｙ2）とする。節点▲３▼の(ｘ,ｙ)座標を、（ｘ3,ｙ3）とする。３角形の面積をＳeとすると、
２Ｓe＝|(1,1,1)T (ｘ1, ｘ2, ｘ3)T (ｙ1,ｙ2,ｙ3)T|
となる。Δｒは、要素の位置ベクトルｒの変化量である。(Δｒ)1は、節点▲２▼から節点▲３▼に向かうベクトルである。(Δｒ)2は、節点▲３▼から節点▲１▼に向かうベクトルである。(Δｒ)3は、節点▲１▼から節点▲２▼に向かうベクトルである。ベクトルｋは、紙面に垂直で上に向かう単位ベクトルである。ξ1座標の原点を辺▲２▼▲３▼上にとり、辺▲２▼▲３▼に垂直で節点▲１▼に向かう方向を正の方向とする。ξ2座標の原点を辺▲１▼▲３▼上にとり、辺▲１▼▲３▼に垂直で節点▲２▼に向かう方向を正の方向とする。ξ3座標の原点を辺▲１▼▲２▼上にとり、辺▲１▼▲２▼に垂直で節点▲３▼に向かう方向を正の方向とする。ξ1を(ｘ,ｙ)座標で表わすと、
ξ1＝ａｘ＋ｂｙ＋ｃ
１＝ａｘ1＋ｂｙ1＋ｃ
０＝ａｘ2＋ｂｙ2＋ｃ
０＝ａｘ3＋ｂｙ3＋ｃ
ａ＝−(ｙ3−ｙ2)／２Ｓe
ｂ＝(ｘ3−ｘ2)／２Ｓe
ｃ＝(ｘ2ｙ3−ｘ3ｙ2)／２Ｓe
となる。したがって、反変基底ベクトルｇ1は、
ｇ1＝∇ξ1
＝−ｉ(ｙ3−ｙ2)／２Ｓe
＋ｊ(ｘ3−ｘ2)／２Ｓe
となる。
【００４５】
ところで、
(Δｒ)1＝ｉ(ｘ3−ｘ2)＋ｊ(ｙ3−ｙ2)
ｉ＝−ｋ×ｊ，ｊ＝ｋ×ｉ
ｋ×(Δｒ)1＝ｊ(ｘ3−ｘ2)−ｉ(ｙ3−ｙ2)
であるから、
ｇ1＝∇ξ1＝ｋ×(Δｒ)1／２Ｓe
となる。ξ2とξ3についても同様に、
ｇ2＝∇ξ2＝ｋ×(Δｒ)2／２Ｓe
ｇ3＝∇ξ3＝ｋ×(Δｒ)3／２Ｓe
となる。ξ1とξ2とξ3は一次独立でなく、
(ｘ3−ｘ2)＋(ｘ2−ｘ1)＋(ｘ1−ｘ3)＝０
(ｙ3−ｙ2)＋(ｙ2−ｙ1)＋(ｙ1−ｙ3)＝０
(ｘ2ｙ3−ｘ3ｙ2)＋(ｘ3ｙ1−ｘ1ｙ3)＋(ｘ1ｙ2−ｘ2ｙ1)＝２Ｓe
であるから、明らかに、
ξ1＋ξ2＋ξ3＝１
である。ｇ1とｇ2とｇ3も一次独立でなく、
(Δｒ)1＋(Δｒ)2＋(Δｒ)3＝０
であるから、明らかに、
ｇ1＋ｇ2＋ｇ3＝０
である。
【００４６】
節点▲１▼の形状関数Ｎ1は、節点▲１▼で１になり、節点▲２▼▲３▼で０になる線形関数であるから、明らかにξ1に一致する。ξ2とξ3についても同様であるから、
Ｎa＝ξa  （ａ＝１，２，３）
と書ける。
【００４７】
φの節点▲１▼における値をφ1とし、節点▲２▼における値をφ2とし、節点３における値をφ3とする。一般式は、
∇a＊φa＝ｇa＊φa
＝ｇ1＊φ1＋ｇ2＊φ2＋ｇ3＊φ3
＝−(ｇ2＋ｇ3)＊φ1−(ｇ3＋ｇ1)＊φ2−(ｇ1＋ｇ2)＊φ3
＝−ｇ1＊(φ2＋φ3)−ｇ2＊(φ3＋φ1)−ｇ3＊(φ1＋φ2)
＝−ｋ×(Δｒ)1／Ｓe＊(φ2＋φ3)／２
−ｋ×(Δｒ)2／Ｓe＊(φ3＋φ1)／２
−ｋ×(Δｒ)3／Ｓe＊(φ1＋φ2)／２
＝(１／Ｓe){Ｓ1＊φ1-＋Ｓ2＊φ2-＋Ｓ3＊φ3-}
となる。ただし、
Ｓ1＝(Δｒ)1×ｋ, φ1-＝(φ2＋φ3)／２
Ｓ2＝(Δｒ)2×ｋ, φ2-＝(φ3＋φ1)／２
Ｓ3＝(Δｒ)3×ｋ, φ3-＝(φ1＋φ2)／２
である。また、
ｆ1-＝(ｆ2＋ｆ3)／２
ｆ2-＝(ｆ3＋ｆ1)／２
ｆ3-＝(ｆ1＋ｆ2)／２
とすると、
∇aｆa＝(１／Ｓe){Ｓ1ｆ1-＋Ｓ2ｆ2-＋Ｓ3ｆ3-}
と書ける。同様に、
∇a・ｖa＝(１／Ｓe){Ｓ1・ｖ1-＋Ｓ2・ｖ2-＋Ｓ3・ｖ3-}
(∇a×ｖa)・ｋ＝(１／Ｓe){(Δｒ)1・ｖ1-＋(Δｒ)2・ｖ2-＋(Δｒ)3・ｖ3-}
と書ける。
【００４８】
図３を参照して、３次元の場合を説明する。▲１▼,▲２▼,▲３▼,▲４▼は、節点番号である。４面体▲１▼▲２▼▲３▼▲４▼を要素とする。節点▲１▼の(ｘ,ｙ,ｚ)座標を、（ｘ1,ｙ1,ｚ1）とする。節点▲２▼の(ｘ,ｙ,ｚ)座標を、（ｘ2,ｙ2,ｚ2）とする。節点▲３▼の(ｘ,ｙ,ｚ)座標を、（ｘ3,ｙ3,ｚ3）とする。節点▲４▼の(ｘ,ｙ,ｚ)座標を、（ｘ4,ｙ4,ｚ4）とする。４面体の体積をＶeとすると、
６Ｖe
＝|(1,1,1,1)T (ｘ1,ｘ2,ｘ3,ｘ4)T (ｙ1,ｙ2,ｙ3,ｙ4)T (ｚ1,ｚ2,ｚ3,ｚ4)T|
となる。節点▲１▼の位置ベクトルをｒ1とし、節点▲２▼の位置ベクトルをｒ2とし、節点▲３▼の位置ベクトルをｒ3とし、節点▲４▼の位置ベクトルをｒ4とする。ξ1を(ｘ,ｙ,ｚ)座標で表わすと、
ξ1＝ａｘ＋ｂｙ＋ｃｚ＋ｄ
１＝ａｘ1＋ｂｙ1＋ｃｚ1＋ｄ
０＝ａｘ2＋ｂｙ2＋ｃｚ2＋ｄ
０＝ａｘ3＋ｂｙ3＋ｃｚ3＋ｄ
０＝ａｘ3＋ｂｙ3＋ｃｚ3＋ｄ
ａ＝−(Ｓ1・ｉ)／３Ｖe
ｂ＝−(Ｓ1・ｊ)／３Ｖe
ｃ＝−(Ｓ1・ｋ)／３Ｖe
ｄ＝|Ｓ1|／３Ｖe
となる。ただし、
Ｓ1＝(ｒ3−ｒ2)×(ｒ4−ｒ2)／２
である。したがって、反変基底ベクトルｇ1は、
ｇ1＝∇ξ1
＝−ｉ(Ｓ1・ｉ)／３Ｖe
−ｊ(Ｓ1・ｊ)／３Ｖe
−ｋ(Ｓ1・ｋ)／３Ｖe
＝−Ｓ1／３Ｖe
となる。ξ2とξ3とξ4についても同様に、
ｇ2＝∇ξ2＝−Ｓ2／３Ｖe
ｇ3＝∇ξ3＝−Ｓ3／３Ｖe
ｇ4＝∇ξ4＝−Ｓ4／３Ｖe
となる。ξ1とξ2とξ3とξ4は一次独立でなく、明らかに、
ξ1＋ξ2＋ξ3＋ξ4＝１
である。ｇ1とｇ2とｇ3とｇ4も一次独立でなく、明らかに、
ｇ1＋ｇ2＋ｇ3＋ｇ4＝０
である。
【００４９】
節点▲１▼の形状関数Ｎ1は、節点▲１▼で１になり、節点▲２▼▲３▼▲４▼で０になる線形関数であるから、明らかにξ1に一致する。ξ2とξ3とξ4についても同様であるから、
Ｎa＝ξa  （ａ＝１，２，３，４）
と書ける。
【００５０】
φの節点▲１▼における値をφ1とし、節点▲２▼における値をφ2とし、節点３における値をφ3とし、節点４における値をφ4とする。一般式は、
Σa∇a＊φa＝Σaｇa＊φa
＝ｇ1＊φ1＋ｇ2＊φ2＋ｇ3＊φ3＋ｇ4＊φ4
＝−(ｇ2＋ｇ3＋ｇ4)＊φ1−(ｇ3＋ｇ4＋ｇ1)＊φ2
−(ｇ4＋ｇ1＋ｇ2)＊φ3−(ｇ1＋ｇ2＋ｇ3)＊φ4
＝−ｇ1＊(φ2＋φ3＋φ4)−ｇ2＊(φ3＋φ4＋φ1)
−ｇ3＊(φ4＋φ1＋φ2)−ｇ4＊(φ1＋φ2＋φ3)
＝Ｓ1／Ｖe＊(φ2＋φ3＋φ4)／３
＋Ｓ2／Ｖe＊(φ3＋φ4＋φ1)／３
＋Ｓ3／Ｖe＊(φ4＋φ1＋φ2)／３
＋Ｓ4／Ｖe＊(φ1＋φ2＋φ3)／３
＝(１／Ｖe){Ｓ1＊φ1-＋Ｓ2＊φ2-＋Ｓ3＊φ3-＋Ｓ4＊φ4-}
となる。ただし、
Ｓ1＝(ｒ3−ｒ2)×(ｒ4−ｒ2)／２, φ1-＝(φ2＋φ3＋φ4)／３
Ｓ2＝(ｒ3−ｒ1)×(ｒ4−ｒ1)／２, φ2-＝(φ3＋φ4＋φ1)／３
Ｓ3＝(ｒ4−ｒ2)×(ｒ1−ｒ2)／２, φ3-＝(φ4＋φ1＋φ2)／３
Ｓ4＝(ｒ3−ｒ1)×(ｒ2−ｒ1)／２, φ4-＝(φ1＋φ2＋φ3)／３
である。また、
ｆ1-＝(ｆ2＋ｆ3＋ｆ4)／３
ｆ2-＝(ｆ3＋ｆ4＋ｆ1)／３
ｆ3-＝(ｆ4＋ｆ1＋ｆ2)／３
ｆ4-＝(ｆ1＋ｆ2＋ｆ3)／３
とすると、
Σa∇aｆa＝(１／Ｖe){Ｓ1ｆ1-＋Ｓ2ｆ2-＋Ｓ3ｆ3-＋Ｓ4ｆ4-}
と書ける。同様に、
Σa∇a・ｖa＝(１／Ｖe){Ｓ1・ｖ1-＋Ｓ2・ｖ2-＋Ｓ3・ｖ3-＋Ｓ4・ｖ4-}
Σa∇a×ｖa＝(１／Ｖe){Ｓ1×ｖ1-＋Ｓ2×ｖ2-＋Ｓ3×ｖ3-＋Ｓ4×ｖ4-}
と書ける。
【００５１】
第４に、対称要素のこう配・発散・回転について説明する。１次元の線分、２次元の４角形、３次元の６面体について調べる。これらの要素は、写像関数
ｒ＝ΣaＮaｒa
により計算空間に写像すると、各座標軸に対して対称となるから、対称要素と呼ぶ。対称要素の形状関数は、それぞれの次元に対して、
Ｎa＝(1/2)(１＋ξaξ)  （１次元）
Ｎa＝(1/4)(１＋ξaξ)(１＋ηaη)  （２次元）
Ｎa＝(1/8)(１＋ξaξ)(１＋ηaη)(１＋ζaζ)  （３次元）
と書く事ができる。ξa,ηa,ζaは、節点により決まる＋１または−１である。１次元の場合、
Ｎ1＝(1/2)(１−ξ)；  ξa＝−１
Ｎ2＝(1/2)(１＋ξ)；  ξa＝＋１
となる。
【００５２】
よって、各次元の要素平均は、
<∂Ｎa／∂ξi>e＝(1/2)ξai  （１次元（ｉ＝１））
<∂Ｎa／∂ξi>e＝(1/4)ξai  （２次元（ｉ＝１，２））
<∂Ｎa／∂ξi>e＝(1/8)ξai  （３次元（ｉ＝１，２，３））
と表示できる。これより、
∇a＊φa＝Σiｇei＊(1/2)ξaiφ  （１次元）
∇a＊φa＝Σiｇei＊(1/4)ξaiφ  （２次元）
∇a＊φa＝Σiｇei＊(1/8)ξaiφ  （３次元）
を得る。１次元、２次元、３次元をまとめて、一般化公式、
∇a＊φa＝Σi(１／Ωe(ξ))ｇei＊(ξaiφa)
が導ける。ここで、Ωe(ξ)は、計算空間の各次元の体積である。１次元の場合の結果を、図４に示す。２次元の場合の結果を、図５に示す、３次元の場合の結果を、図６に示す。
【００５３】
特に、２次元の４角形要素に対して、各辺に対する面積Ｓi（ｉ＝１，２，３，４）を、次のように定義する。
(Δｒ)1×ｋ＝Ｓ1    ｋ×(Δｒ)3＝Ｓ3
(Δｒ)2×ｋ＝Ｓ2    ｋ×(Δｒ)4＝Ｓ4
このとき、面積分公式
<∇＊φ>e＝Σa∇a＊φa＝ΣsΔΓs(ｘ)＊φs／Ωe(ｘ)
を用いると、発散と回転は、
Σa∇a＊ｖa
＝（１／Ｓe）[Ｓ1・ｖ1*＋Ｓ2・ｖ2*＋Ｓ3・ｖ3*]
Σa(∇a×ｖa)・ｋ
＝（１／Ｓe）[(Δｒ)1・ｖ1*＋(Δｒ)2・ｖ2*
−(Δｒ)3・ｖ3*−(Δｒ)4・ｖ4*]
と書くこともできる。
【００５４】
平行４辺形の場合は、この結果が、体積分公式で求めた結果
Σa∇a＊ｖa
＝（１／Ｓe）[Ｓ1*・δｖ1*＋Ｓ2*・δｖ2*]
Σa(∇a×ｖa)・ｋ
＝（１／Ｓe）[Ｓ1*×δｖ1*＋Ｓ2*×δｖ2*]・ｋ
に一致することが証明できる。３次元の場合も同様である。すなわち、平行６面体の場合は、面積分公式と体積分公式の結果は一致する。そして、ゆがんだ要素に付しては、計算が少し複雑になるけれども、精度を重視すれば、面積分公式を用いるべきである。
【００５５】
第５に、係数行列の解析的表示と非記憶化について説明する。電磁熱流体の有限要素法解析には、多くの係数行列が表れる。これらの係数行列を、表４と表５に列挙しておく。２階微分の係数行列であるOpen行列、Cross行列、Dot行列の離散nabla演算子表示を、表４に示す。また、移流行列、拡散行列、BTD行列の離散nabla演算子を、表５に示す。ここで、curl curl演算子より生ずるCross行列の場合は、積の順序に注意しなければならない。また、Laplacianより生ずるDot行列は、拡散行列とも呼ばれる。BTD行列は、移流項の上流化より生ずる行列である。
【００５６】
ただし、係数行列の離散nabla演算子の積表示は、hourglass項も含めた定義になっている。例えば、拡散行列において
Ｄab＝Ωe∇a・∇b
＝Ｊeｇei・ｇej<(∂Ｎa／∂ξi)(∂Ｎb／∂ξj)>e
＝Ｊeｇei・ｇej<∂Ｎa／∂ξi>e<∂Ｎb／∂ξj>e＋(hourglass項)
と定義されている。すなわち、∇aと∇bを独立に内積すると、primary partにhourglass項（secondary part）を付加する必要がある。係数行列については、周知のものであるので、これ以上の詳しい説明は省略するが、必要であれば、文献１などを参照されたい。
【００５７】
［表４］  係数行列と離散nabla演算子（その１  ２階微分）
Open行列  非対称行列
Ｏab＝∫∇Ｎa∇ＮbdΩ(ｒ)
＝Ωe(ｒ)∇a∇b
＝ΣiΣjΩe(ｒ)ｇeiｇej<(∂Ｎa／∂ξi)(∂Ｎb／∂ξj)>
Cross行列  反対称行列
Ｒab＝∫∇Ｎa×∇ＮbdΩ(ｒ)
＝Ωe(ｒ)∇a×∇b
＝ΣiΣjΩe(ｒ)ｇei×ｇej<(∂Ｎa／∂ξi)(∂Ｎb／∂ξj)>
Dot行列  対称行列
Ｄab＝∫∇Ｎa・∇ＮbdΩ(ｒ)
＝Ωe(ｒ)∇a・∇b
＝ΣiΣjΩe(ｒ)ｇei・ｇej<(∂Ｎa／∂ξi)(∂Ｎb／∂ξj)>
【００５８】
［表５］  係数行列と離散nabla演算子（その２）
質量行列（対称行列）
Ｍab＝∫ＮaＮbdΩ(ｒ)
集中質量行列（ｎは要素の節点数）
Ｍ*ab＝Ｊeδab    Ｊe＝Ωe(ｒ)／ｎ
移流行列（非対称行列）
Ａab＝∫Ｎaｖ・∇ＮbdΩ(ｒ)
＝ΣjΩe(ｒ)ｖe・ｇej<Ｎa(∂Ｎb／∂ξj)>
拡散行列（Dot行列）
Ｄab＝∫∇Ｎa・∇ＮbdΩ(ｒ)
＝Ωe(ｒ)∇a・∇b
＝ΣiΣjΩe(ｒ)ｇei・ｇej<(∂Ｎa／∂ξi)(∂Ｎb／∂ξj)>
BTD行列（対称行列）
Ｂab＝∫ｖ・∇Ｎaｖ・∇ＮbdΩ(ｒ)
＝Ωe(ｒ)ｖe・∇aｖe・∇b
＝ΣiΣjΩe(ｒ)ｖe・ｇeiｖe・ｇej<(∂Ｎa／∂ξi)(∂Ｎb／∂ξj)>
計算に必要なものは、<∂ｒ／∂ξi>(ｉ＝１，２，３)及び<Ｎa(∂Ｎb／∂ξj)>,<(∂Ｎa／∂ξi)(∂Ｎb／∂ξj)>だけである。
【００５９】
この係数行列を記憶すると、莫大なメモリーが必要とされる。記憶しないで毎回すべて計算すると、計算時間が長くなる。よって、非記憶化と高速計算法が必要である。そのために、行列計算において、行列要素を２つの部分に分解し、積和で表示する。すなわち、
（係数行列要素）＝Σ（係数行列基底）×〈部分係数因子〉
とする。
【００６０】
〈部分係数因子〉は、各要素の形状によらない部分で、解析的に求める。計算空間で計算できる〈部分係数因子〉は、基本的には、<(∂Ｎa／∂ξi)(∂Ｎb／∂ξj)>e，<Ｎa(∂Ｎb／∂ξi)>eの２個のみである。これらは必要になる度に生成し非記憶化する。一方、残りの部分を(係数行列基底)と呼ぶ。(係数行列基底)は、要素形状に依存する部分である。これは、共変基底ベクトルｇi（≡∂ｒ／∂ξi(ｉ＝１，２，３)）を用いて容易に計算できる。(係数行列基底)は、要素の種類だけ計算する。例えば、全体が８種類の要素から構成されていれば、８個の要素について計算しておけばよい。１億個の要素があっても、皆同じ形状の要素であるならば、ただ１個の要素についての情報<∂ｒ／∂ξi>e(ｉ＝１，２，３)が記憶されていればよい。この部分を係数行列基底記憶装置に格納しておく。行列計算で係数行列要素が必要になる時、(係数行列基底)を読み出し、〈部分係数因子〉を生成し、積和演算で合成して係数行列要素を生成することにより、非記憶化する。
【００６１】
第６に、反変基底ベクトルのObject化について説明する。離散nabla演算子は、
∇e≡Σiｇei<∂Ｎa／∂ξi>e
で定義されている。ここで、反変基底ベクトルｇeiの計算が重要である。そして、反変基底ベクトルｇeiは、有限体積法と有限要素法に共通である。要素の種類が増加すると、計算量も増加する。かつ、記憶容量も増大する。そこで、要素の種類とその座標が与えられたら、この反変基底ベクトルｇeiを計算するルーチンをObject化し、かつLSI化して、メッセージとして節点座標を送ればｇeiを得られるように、組込み関数的に利用できるようにする。３次元の４／５／６面体要素に対して、チップ化すれば、混合要素での解析が可能となり、格子生成が極めて容易になる。これがデルチップである（図７参照）。名前デルは、nabla演算子の別名である。
【００６２】
図８と図９を参照して、離散nabla演算子を使った有限要素法の非記憶型計算手順を説明する。ステップ１で、偏微分方程式の種類を選択し、パラメータと境界条件を入力する。偏微分方程式情報記憶装置６に格納してある計算可能な偏微分方程式に関する情報を表示装置７に表示し、入力装置１から、必要な情報を入力する。入力された偏微分方程式のパラメータと境界条件のデータをパラメータ境界条件記憶装置２に格納しておく。ステップ２で、節点位置のベクトルを、入力装置１から入力する。これは、必ずしも個別データを入力することだけではなく、自動的に格子を生成する場合も含む。節点位置ベクトル記憶装置３に、すべての節点の座標値を格納しておく。
【００６３】
ステップ３で、既知ベクトル演算手段４により、パラメータと境界条件から既知ベクトルを計算して、既知ベクトル記憶装置５に、既知ベクトルを格納しておく。ステップ４で、係数行列基底演算手段８により、節点座標から反変基底ベクトルと要素体積などを計算し、係数行列基底記憶装置９に、係数行列基底を格納しておく。この係数行列基底演算手段８として、図７のデルチップ72を利用することができる。
【００６４】
ステップ５で、解析的に計算した形状関数のデータから部分係数因子生成手段11により生成した部分係数因子と、係数行列基底記憶装置９から読み出した係数行列基底から、係数行列要素演算手段12により、係数行列要素を計算する。係数行列要素を、行列要素記憶ワークエリア14の一時的に格納する。ステップ６で、行列演算手段13により、係数行列要素と既知ベクトルから解ベクトルを求める計算を行なう。このステップ５とステップ６を繰り返して、係数行列基底と部分係数因子から計算した係数行列要素を、ワークエリアを使って一時的に記憶して行列計算を行なうことを続けることにより、解を求める。ステップ７で、ＣＲＴディスプレイなどの表示装置７に解析結果を表示する。このように、行列計算において必要になるたびに、係数行列要素を計算するので、係数行列全体を記憶するメモリーは必要ない。係数行列要素は、（節点数）×（節点数）だけあるが、０でない値をとるのは、１行につき数十個程度で、対角成分付近に集中しているので、ワークエリアも少なくて済む。
【００６５】
上記のように、本発明の実施の形態では、非構造格子の有限体積法または有限要素法により計算する数値解析装置を、離散nabla演算子で記述された係数行列要素の要素形状に依存する因子を記憶し、要素形状に依存しない因子を解析的に計算しておき、必要になるたびに係数行列要素をこれから生成する構成としたので、記憶容量を少なくしながらも、計算を高速化できる。
【００６６】
【発明の効果】
以上の説明から明らかなように、本発明では、非構造格子の有限体積法または有限要素法により偏微分方程式の境界値問題の解を計算する数値解析装置を、離散nabla演算子で表わした係数行列要素のうちの要素形状に依存する因子を要素形状に対応させて記憶する係数行列基底記憶手段と、離散nabla演算子で表わした係数行列要素のうちの要素形状に依存しない因子を種類別に非記憶化し部分係数因子を生成する部分係数因子生成手段と、行列計算で必要になるたびに係数行列基底記憶手段を参照して得た係数行列基底と部分係数因子生成手段で生成された部分係数因子とから係数行列要素を生成する計算手段とを具備する構成としたので、記憶容量を少なくしながらも、計算を高速化できるという効果が得られる。
【図面の簡単な説明】
【図１】本発明の実施の形態における数値解析装置で利用する１次元単体要素のこう配・発散・回転の説明図、
【図２】本発明の実施の形態における数値解析装置で利用する２次元単体要素のこう配・発散・回転の説明図、
【図３】本発明の実施の形態における数値解析装置で利用する３次元単体要素のこう配・発散・回転の説明図、
【図４】本発明の実施の形態における数値解析装置で利用する１次元の対称要素のこう配・発散・回転の説明図、
【図５】本発明の実施の形態における数値解析装置で利用する２次元の対称要素のこう配・発散・回転の説明図、
【図６】本発明の実施の形態における数値解析装置で利用する３次元の対称要素のこう配・発散・回転の説明図、
【図７】本発明の実施の形態における数値解析装置のデルチップの概念図である。
【図８】本発明の実施の形態における数値解析装置の処理手順を示す流れ図、
【図９】本発明の実施の形態における数値解析装置の機能ブロック図である。
【符号の説明】
１  入力装置
２  パラメータ境界条件記憶装置
３  節点位置ベクトル記憶装置
４  既知ベクトル演算手段
５  既知ベクトル記憶装置
６  偏微分方程式情報記憶装置
７  表示装置
８  係数行列基底演算手段
９  係数行列基底記憶装置
10  形状関数記憶装置
11  部分係数因子生成手段
12  係数行列要素演算手段
13  行列演算手段
14  行列要素記憶ワークエリア
71  要素の節点座標生成手段
72  デルチップ
73  データ利用手段
【図１】
0004482646.tif 000002
【図２】
0004482646.tif 000003
【図３】
0004482646.tif 000004
【図４】
0004482646.tif 000005
【図５】
0004482646.tif 000006
【図６】
0004482646.tif 000007
【図７】
0004482646.tif 000008
【図８】
0004482646.tif 000009
【図９】
0004482646.tif 000010