(57)【特許請求の範囲】
【請求項１】
非圧縮流体に関して、ソレノイダル・ベクトル場の保存形Navier-Stokes運動方程式を、速度を陽的に圧力を陰的に離散化した方程式を用いて、サイクル誤差自己調整法を用いた同時緩和法で時間進行を行なって次の時間ステップの速度と圧力をコンピュータにより求める有限要素法の数値計算方法において、前記時間進行の処理では、前記コンピュータの時間進行手段により、データ再構築手段を用いて、質量行列の集中化を行ってEuler前進法により速度の予測子を求め、速度の発散が０という条件の修正速度ポテンシャルに関するPoisson方程式をNewton-Raphson法演算手段によりNewton-Raphson法で解いて得た修正速度ポテンシャルから速度と圧力を求める計算を収束するまで繰り返すように演算制御手段により前記時間進行手段と前記Newton-Raphson法演算手段と前記データ再構築手段とを制御することを特徴とする数値計算方法。
【請求項２】
磁場の誘導方程式を、磁場を陽的に離散化した方程式を用いて、サイクル誤差自己調整法を用いた同時緩和法で時間進行を行なって次の時間ステップの磁場と数値残差をコンピュータにより求める有限要素法の数値計算方法において、前記時間進行の処理では、前記コンピュータの時間進行手段により、データ再構築手段を用いてEuler前進法により磁場の予測子を求め、修正磁場ポテンシャルに関するPoisson方程式をNewton-Raphson法演算手段によりNewton-Raphson法で解いて得た修正磁場ポテンシャルから磁場と数値残差を求める計算を収束するまで繰り返すように演算制御手段により前記時間進行手段と前記Newton-Raphson法演算手段と前記データ再構築手段とを制御することを特徴とする数値計算方法。
【請求項３】
定常磁場の電位に関するPoisson方程式を、電流場を陽的に電場を陰的に離散化した方程式を用いて、サイクル誤差自己調整法を用いた同時緩和法で電流場と電場をコンピュータにより求める有限要素法の数値計算方法において、前記コンピュータの時間進行手段により、データ再構築手段を用いてEuler前進法により電流場と電場の予測子を求め、修正電位ポテンシャルに関するPoisson方程式を、Newton-Raphson法演算手段によりNewton-Raphson法により電流場の発散が０という条件で解き、修正電位ポテンシャルから電流場と電場を求める計算を、収束するまで繰り返すように演算制御手段により前記時間進行手段と前記Newton-Raphson法演算手段と前記データ再構築手段とを制御することを特徴とする数値計算方法。
【発明の詳細な説明】
【０００１】
【発明の属する技術分野】
本発明は、ソレノイダル・ベクトル場の有限要素法による数値計算方法に関し、特に、時間積分法に関して質量集中化とデータの再構築を適用して、Poisson方程式の解法に関してサイクル誤差自己調整法を用いた同時緩和法により数値計算を行なう方法に関する。
【０００２】
【従来の技術】
流体力学や構造力学などにおける問題に現れる偏微分方程式の境界値問題を、計算機による数値計算で解く方法の１つに有限要素法がある。有限要素法について、簡単に説明する。詳細については、文献１［棚橋隆彦著「流れの有限要素法解析Ｉ，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を第ｉ,ｊ成分とする行列の逆行列である。このようにして、連立一次方程式を解くことによって、偏微分方程式の境界値問題を解くことができる。
【００１１】
上記の例は、時間的に変化しない定常状態であり、連立一次方程式を解くことで解が求まる場合であるが、時間的に変化し、非線形方程式を解く必要がある場合について、非圧縮性流体のNavier-Stokesの方程式を解く場合を例にして簡単に説明する。詳細は文献１を参照されたい。速度をｖとし、圧力をｐとし、Reynolds数をＲeとすると、Navier-Stokesの方程式は、
∂ｖ／∂ｔ＋ｖ・∇ｖ＝−∇ｐ＋(１／Ｒe)∇2ｖ
∇・ｖ＝０
と記述できる。ｐは、Poisson方程式
∇2ｐ＝∇・(ｖ・∇ｖ)
を満たす。この方程式を弱形式にして空間離散化すると、非線形の連立方程式が得られる。解析対象領域を要素分割し、近似解を要素ごとに定義された関数によって表現することを離散化という。時間について離散化して、境界条件と初期条件を入れて、Newton-Raphson法により解くと、ｔ＝Δｔの解が得られる。この解から次の時刻の状態を求めるという時間進行を実行する。時間進行法については、例えば、前進Euler法を使う。これを繰り返すことにより、ｔ＝０〜ＮΔｔまでの状態を求めることができる。
【００１２】
この際、連立方程式の係数として質量行列がでるので、質量を集中化して質量行列を対角化すると、計算が簡単になる。要素形状が一様であれば、質量を集中化しても誤差は混入しない。しかし、大小の要素が混在する場合に質量を集中化すると誤差が入る。質量集中化の誤差を除去するのには、例えば、マルチパスアルゴリズムがある。
【００１３】
このように、従来の有限要素法では、微分形のNavier-Stokesの方程式から出発して、有限要素法で離散化し、時間進行法の際、質量を集中化している。単一格子を用いた有限要素法では、時間進行法においてnode by nodeのvertex centered法が採用され、有限体積法では、element by elementのcell centered法が採用されている。
【００１４】
ところで、ベクトルの発散が０になる場を、ソレノイダル・ベクトル場あるいは単にソレノイダル場という。このソレノイダル・ベクトル場から、一般にPoisson方程式が導かれる。複雑形状の大規模問題において、このPoisson方程式を、有限要素法で高速・高精度かつ短時間に解くことが大きな課題である。従来、差分法のHSMAC法などで解かれていた。
【００１５】
いくつかのソレノイダル・ベクトル場の例をあげる。非圧縮流体の場では、ｖを速度ベクトルとすると、
∇・ｖ＝０
となる。磁場では、Ｂを磁束密度ベクトルとすると、
∇・Ｂ＝０
となる。電場では、Ｅを電界ベクトルとすると、
∇・Ｅ＝０
となる。電流場では、ｊを電流密度ベクトルとすると、
∇・ｊ＝０
となる。Coulombゲージ場では、Ａを磁気ベクトル・ポテンシャルとすると、
∇・Ａ＝０
となる。
【００１６】
ソレノイダル・ベクトル場の方程式を有限要素法で計算する場合には、Poisson方程式を直接解かずに、時間発展方程式と連立させて、ソレノイダル・ベクトル場を求める。このとき、Poisson方程式を、同時緩和法の反復手段に用いる。以下、サイクル誤差自己調整法を用いた同時緩和法を利用して解く方法を説明する。詳細については、文献１を参照されたい。
【００１７】
非圧縮流体の流れ場を例にして、解析対象の基礎となる方程式を説明する。非圧縮流体の運動を記述する方程式は、Navier-Stokesの方程式
(∂ｖ／∂ｔ)＋ｖ・∇ｖ＝−∇ｐ＋ν∇2ｖ
と、連続の方程式
∇・ｖ＝０
からなる。ここで、ｖは速度ベクトルであり、ｐは密度で基準化された圧力であり、νは動粘性係数である。ただし、外力項は無視している。この解析手法では、サイクル誤差自己調整法を用いるから、保存形のNavier-Stokes方程式
(∂ｖ／∂ｔ)＋∇・(ｖｖ)＝−∇ｐ＋∇・τ
τ＝ν[∇ｖ＋(∇ｖ)T]
を用いる。これを原理式とよぶ。ここで、τは粘性応力テンソルである。
【００１８】
速度の発散が完全に０の場合、非保存形と保存形は等価であるが、数値解析上では誤差が発生し、
｜∇・ｖ｜＝ε
となる。例えば、単位量を１とした場合、ε＝10-3である。この誤差がある場合、非保存形と保存形の運動方程式の意味は異なる。保存形は、速度の誤差を含めて運動量がバランスするけれども、非保存形ではアンバランスになる。すなわち、非保存形方程式では数値的な不安定性が発生し易い。しかし、保存形方程式では安定化する。サイクル誤差自己調整法を用いる場合、運動方程式として保存形の原理式を用いることが必須である。このことは、他のソレノイダル場についても同様で、サイクル誤差自己調整法として時間発展方程式を用いる場合、そのソレノイダル場の保存形の原理式を使用することが重要である。
【００１９】
参考のために、粘性項と対流項に用いられるいろいろなFormを、表１に示しておく。
（表１）
0004649673.tif 000002
【００２０】
Poisson方程式の解法について説明する。MAC（Marker and Cell）法では、
∇・ｖn+1＝０
の条件より、圧力ｐに関するPoisson方程式が導かれる。肩付きの(ｎ＋１)は、時刻を示す添え字である。紛れるおそれはないので、時刻を示すｎと法線ベクトルを示すｎを同じ文字で表わす。ｖn+1は、時刻Δｔ（ｎ＋１）における速度である。SMAC（Simplified MAC）法では、同じ条件より、修正速度ポテンシャルφに関するPoisson方程式が導かれる。修正速度ポテンシャルφについては、
φ＝(ｐn+1−ｐn)Δｔ
である。そして、HSMAC（Hybrid-SMAC）法では、このPoisson方程式が、速度と圧力の同時緩和法で解かれる。同時緩和法については、後で説明する。GSMAC（Generalized SMAC）有限要素法では、修正速度ポテンシャルφに関するPoisson方程式が、
∇2φ＝∇・ｖ~
と導かれる。このPoisson方程式が、速度と圧力の同時緩和法で解かれる。速度の予測子ｖ~は、Navier-Stokesの方程式から、Euler前進法により、
（ｖ~−ｖn）／Δｔ＋∇・(ｖｖ)n＝−∇ｐn＋∇・τn
と求める。
【００２１】
ここで、Poisson方程式の具体的な解法を示す。Poisson方程式を、要素内で積分する。重み関数を１とした非構造格子の有限体積法を適用する。ガウスの定理を適用すると、Poisson方程式は、
∬(∂φ／∂ｎ)dＳ＝∬ｖ~・ｎdＳ  （積分範囲はＳe）
（ｎは面の外向き単位法線ベクトル）
となる。ここで、上式の右辺に着目して、
ε＝∬ｖ~・ｎdＳ  （積分範囲はＳe）
を定義する。εは、要素内から流出する流量である。そして、εはｐの関数である。ｐを修正した結果、速度場がソレノイダル場になれば、
ε(ｐ＋δｐ)＝０
を満足する。Taylor展開して修正量δｐを求めると、
δｐ＝−ε／(∂ε／∂ｐ)
となる。これが、Newton-Raphson法の原理である。
【００２２】
微係数∂ε／∂ｐを、力積の関係式
δｖ／Δｔ＝−∇δｐ
を用いて求める。その結果、
∂ε／∂ｐ＝(δ／δｐ)(∬ｖ・dＳ)（積分範囲はＳe）
＝(１／δｐ)(∬δｖ・ｎdＳ)（積分範囲はＳe）
＝(Δｔ／δｐ)(∬(δｐ／Δｎ)dＳ)（積分範囲はＳe）
＝Δｔ∬(１／Δｎ)dＳ（積分範囲はＳe）
(∵δｖn／Δｔ＝−∂δｐ／∂ｎ＝−(０−δｐ)／Δｎ)
を得る。よって、
λ＝(１／Ｖe)∬(１／Δｎ)dＳ（積分範囲はＳe）
と置けば、微係数は、
∂ε／∂ｐ＝λＶeΔｔ
となる。
【００２３】
つぎに、λの計算方法を説明する。図９に示すように、Ｋを対象要素の重心とし、Ｌを隣接要素の重心とし、ｎを面の外向き単位法線ベクトルとすれば、法線方向の増分は、
Δｎ＝(Ｋ→Ｌ)・ｎ
と評価できる。ここで、dＳ＝ΔＳは、面の微小面積である。（Ｋ→Ｌ）は、ＫからＬへのベクトルである。図面では、上付き矢印で示した。よって、
∬(１／Δｎ)dＳ＝Σ(ΔＳ／Δｎ)＝Σ(ΔＳ／((Ｋ→Ｌ)・ｎ))
となる。Σは、それぞれの面について総和をとるものとする。面の単位法線ベクトルは、要素の種類に応じて、図８のようになる。これらの結果は、∇2φの優対角近似に一致する。
【００２４】
優対角近似について説明する。GSMAC有限要素法は、HSMAC法の有限要素法への拡張である。したがって、要素の形状に依存しない形で優対角近似をすることが必要である。まず、Laplacianの定義式を説明する。Laplacianは、任意の体積要素Ｖに対して、
∇2φ＝limV→0(１／Ｖ)∬(∂φ／∂ｎ)dＳ（積分範囲はＳe）
で定義される。これを、要素Ｖeに適用し、修正速度ポテンシャルが０（要素外）からφ（要素内）に増加した場合を考えると、
∂φ／∂ｎ＝(０−φ)／Δｎ
であるから、各要素に対して、
−∇2φ＝(１／Ｖe)(∬(１／Δｎ)dＳ)φ（積分範囲はＳe）
＝λφ
となる。ここで、φは要素内で一定である。
【００２５】
図10に示す６面体要素Ｋに対して、次の積分公式が成立する。
∬(∂φ／∂ｎ)
＝Ａ1(φL−φK)＋Ａ2(φa−φc)＋Ａ3(φb−φd)
ここで係数は、
Ａ1＝(１／12Ｖ)((ａ→ｃ)×(ｂ→ｄ))・((ａ→ｃ)×(ｂ→ｄ))
Ａ2＝(１／12Ｖ)((ａ→ｃ)×(ｂ→ｄ))・((ｂ→ｄ)×(Ｌ→Ｋ))
Ａ3＝(１／12Ｖ)((ａ→ｃ)×(ｂ→ｄ))・((Ｌ→Ｋ)×(ａ→ｃ))
である。Ｖは８面体の体積で、
Ｖ＝(１／６)((ａ→ｃ)×(ｂ→ｄ))・(Ｋ→Ｌ)
である。この公式を用いて、点SOR(Successive Over Relaxation)法で、Poisson方程式を解くこともできるが、ここでは、優対角近似法を用いる。直交座標においては、優対角近似法と点SOR法が一致することが証明できる。
【００２６】
優対角近似の具体的な計算方法を説明する。前述の式の右辺の面積分を実行すると、
∬(∂φ／∂ｎ)dＳ
＝Σi=16{Ａ1i(φLi−φKi)＋Ａ2i(φai−φci)＋Ａ3i(φbi−φdi)}
の形となる。ここで、ｉ＝１〜６は、６面体の番号を表している。よって、Laplace演算子の優対角近似は、
∇2φ＝(１／Ｖe)(Σi=16Ａ1iφ)
と書ける（図10参照）。これより、Laplace演算子の優対角近似は、体積要素に対して、
−∇2＝(１／Ｖe)∬(１／Δｎ)dＳ（積分範囲はＳe）
＝(１／Ｖe)(Σi=16Ａ1i)
＝λ
となる。特に、２次元の直交要素に対して、
λ＝２((１／Δｘ2)＋(１／Δｙ2))
となる。これは、SOLA法に一致する。
【００２７】
したがって、修正速度ポテンシャルφに関するPoisson方程式
∇2φ＝∇・ｖ~
は、差分法で優対角近似を用いて解かれ、非構造格子を用いるGSMAC有限要素法において、Laplacianの逆演算は、
φ＝−λ-1Ｄ~
と近似できる。ただし、
λ＝(１／Ｖe)∬(１／Δｎ)dＳ（積分範囲はＳe）
Ｄ~＝(１／Ｖe)∬(ｖ~・ｎ)dＳ（積分範囲はＳe）
である。
【００２８】
ここで、同時緩和法について説明する。同時緩和法は、
∇・ｖ＝０，∇・Ｂ＝０，∇・ｊ＝０
を満足するソレノイダル場から導かれるPoisson方程式の解法に対して有効である。そして、速度場からは修正速度ポテンシャルに関するPoisson方程式が、磁場からは修正磁場ポテンシャルに関するPoisson方程式が、電流場からは電位ポテンシャルに関するPoisson方程式が、それぞれ導かれる。
【００２９】
修正速度ポテンシャルに関するPoisson方程式を説明する。まず、速度の予測値ｖ~を求める。つぎに、修正速度ポテンシャルφを導入して、速度と圧力を同時緩和法で、
∇・ｖn+1＝０
を満足させながら解く。すなわち、初期値（反復子ｋ＝０のときの値）を、
ｖ(0)＝ｖ~、ｐ(0)＝ｐn
とする。反復子ｋを０から１づつ増加させながら、つぎの計算を繰り返す。φ(k)に関するPoisson方程式
∇2φ(k)＝∇・ｖ(k)
を、Newton-Raphson法により解き、求めたφ(k)から、
ｖ(k+1)＝ｖ(k)−∇φ(k)
ｐ(k+1)＝ｐ(k)＋φ(k)／Δｔ
を求める。反復の結果、すべての要素で、
｜∇・ｖ(k+1)｜≦εv（例えば、εv＝0.001）
を満足したとき、
ｖn+1＝ｖ(k+1)、ｐn+1＝ｐ(k+1)
として、つぎの時間ステップに進む。
【００３０】
この同時緩和法の計算手順を簡略化したアルゴリズムで示すと、次のようになる。図11にPAD形式で示す。
0004649673.tif 000003
【００３１】
修正磁場ポテンシャルに関するPoisson方程式を説明する。まず、磁場の予測値Ｂ~を求める。つぎに、修正磁場ポテンシャルψを導入して、磁場とスカラー関数Ｒを同時緩和法で、
∇・Ｂn+1＝０
を満足させながら解く。すなわち、初期値（反復子ｋ＝０のときの値）を、
Ｂ(0)＝Ｂ~、Ｒ(0)＝Ｒn
とする。反復子ｋを０から１づつ増加させながら、つぎの計算を繰り返す。ψ(k)に関するPoisson方程式
∇2ψ(k)＝∇・Ｂ(k)
を、Newton-Raphson法により解き、求めたψ(k)から、
Ｂ(k+1)＝Ｂ(k)−∇ψ(k)
Ｒ(k+1)＝Ｒ(k)＋ψ(k)／Δｔ
を求める。反復の結果、すべての要素で、
｜∇・Ｂ(k+1)｜≦εB（例えば、εB＝0.001）
を満足したとき、
Ｂn+1＝Ｂ(k+1)、Ｒn+1＝Ｒ(k+1)
として、つぎの時間ステップに進む。すなわち、この解法により、従来困難とされていた磁場のソレノイダル性が保証される。この場合、記号の対応（表２）を用いて、速度場と同じサブルーチンが使用できる（図12参照）。
（表２）
0004649673.tif 000004
【００３２】
電磁場が定常と見なされる場合の電位に関するPoisson方程式を解く方法を説明する。最初に、電位に関するPoisson方程式を近似的に解いて、電位の予測子Ｖ~を求める。つぎに、Ohmの法則
ｊ~＝−∇Ｖ~＋ｖ×Ｂ
により、電流の予測子ｊ~を求める。この予測子ｊ~は当然、
∇・ｊ~≠０
である。そこで、修正電位ポテンシャルχを導入して、修正電位ポテンシャルχを、次の反復法により求める。初期値（反復子ｋ＝０のときの値）を、
ｊ(0)＝ｊ~，Ｖ(0)＝Ｖ~
とする。反復子ｋを０から１づつ増加させながら、つぎの計算を繰り返す。χ(k)に関するPoisson方程式
∇2χ(k)＝∇・ｊ(k)
を、Newton-Raphson法により解き、求めたχ(k)から、
ｊ(k+1)＝ｊ(k)−∇χ(k)
Ｖ(k+1)＝Ｖ(k)＋χ(k)
を求める。反復の結果、すべての要素で、
｜∇・ｊ(k+1)｜≦εj（例えば、εj＝0.001）
を満足したとき、
ｊ＝ｊ(k+1)、Ｖ＝Ｖ(k+1)
として、反復を終了する。この結果、ソレノイダル性を満足する電流場が求められる。解法のアルゴリズムは、速度場の場合と同じである。
【００３３】
【発明が解決しようとする課題】
しかし、従来の有限要素法による数値計算方法では、要素の節点に質量を集中化するという離散化をしていたので、時間進行法による計算の際、誤差が混入して正確な計算ができないという問題があった。
【００３４】
本発明は、上記従来の問題を解決して、ソレノイダル・ベクトル場のNavier-Stokesの方程式を有限要素法で解く場合に、計算精度を上げるとともに解の収束を安定化することを目的とする。
【００３５】
【課題を解決するための手段】
上記の課題を解決するために、本発明では、ソレノイダル・ベクトル場の保存形Navier-Stokes運動方程式を、速度を陽的に圧力を陰的に離散化し、サイクル誤差自己調整法を用いた同時緩和法で時間進行を行なって次の時間ステップの速度と圧力を求める有限要素法の数値計算方法の時間進行の処理で、速度の予測子を質量行列の集中化とデータの再構築法を用いてEuler前進法により求め、速度の発散が０という条件の修正速度ポテンシャルに関するPoisson方程式をNewton-Raphson法で解いて得た修正速度ポテンシャルから速度と圧力を求める計算を収束するまで繰り返す構成とした。
【００３６】
このように構成したことにより、ソレノイダル・ベクトル場のNavier-Stokesの方程式を有限要素法で解く場合に、計算精度が上がるとともに、解が安定に収束するようになる。この数値計算方法は、すべてのソレノイダル・ベクトル場に適用できる。
【００３７】
【発明の実施の形態】
以下、本発明の実施の形態について、図１〜図７を参照しながら詳細に説明する。
【００３８】
（第１の実施の形態）
本発明の第１の実施の形態は、ソレノイダル・ベクトル場の保存形Navier-Stokes運動方程式を、速度を陽的に、圧力を陰的に離散化し、サイクル誤差自己調整法を用いた同時緩和法で、速度の予測子を、質量行列の集中化とデータの再構築法を用いて、Euler前進法により求め、速度の発散が０という条件の修正速度ポテンシャルに関するPoisson方程式をNewton-Raphson法で解き、修正速度ポテンシャルから速度と圧力を求めることを、収束するまで繰り返して、次の時間ステップの速度と圧力を求める時間進行を行なう数値計算方法である。
【００３９】
図１は、本発明の第１の実施の形態における数値計算方法の概略の流れ図である。図２は、時間進行の流れ図である。図３は、数値計算方法を実行する数値計算装置の機能ブロック図である。図３において、入力装置１は、パラメータなどを入力するキーボードなどの手段である。記憶装置２は、パラメータや境界条件などを記憶するメモリである。出力装置３は、計算結果を出力する表示装置などの出力手段である。演算制御手段４、全体の演算動作を制御する手段である。Newton-Raphson法演算手段５は、Poisson方程式を解く手段である。時間進行手段６は、Δｔ時間後の状態を計算する演算手段である。データ再構築手段７は、節点のデータから重心または図心のデータを求めるとともに、重心または図心のデータから節点のデータを再構成する手段である。
【００４０】
上記のように構成された本発明の第１の実施の形態における数値計算方法の手順を説明する。最初に、計算対象の基礎となる方程式を説明する。非圧縮流体の運動を記述する方程式は、Navier-Stokesの方程式
(∂ｖ／∂ｔ)＋ｖ・∇ｖ＝−∇ｐ＋ν∇2ｖ
と、連続の方程式
∇・ｖ＝０
からなる。ここで、ｐは密度で基準化された圧力であり、νは動粘性係数である。ただし、外力項は無視している。サイクル誤差自己調整法を用いるから、粘性項に応力発散形を用いた保存形のNavier-Stokes方程式
(∂ｖ／∂ｔ)＋∇・(ｖｖ)＝−∇ｐ＋∇・τ
τ＝ν[∇ｖ＋(∇ｖ)T]
を用いる。τは粘性応力テンソルである。
【００４１】
次に、図１を参照しながら、数値計算方法の概略手順を説明する。ステップ１で、保存形のNavier-Stokes運動方程式を、速度を陽的に、圧力を陰的に離散化する。ステップ２で、ｔ＝０の初期値で初期化する。ステップ３で、ｔ＝ｎΔｔの状態から、ｔ＝（ｎ＋１）Δｔの状態を求める時間進行を行なう。時間進行の具体的な処理は、図２の流れ図に示す。ステップ４で、所定回数Ｎだけ実行したか調べ、まだであればステップ３を実行し、Ｎ回実行したら終了する。
【００４２】
図２を参照しながら、時間進行の手順を説明する。ステップ11で、速度の予測子ｖ~を、時刻ｔ＝ｎΔｔにおける速度ｖnと圧力ｐnから求める。このとき、質量行列の集中化とデータの再構築法を用いる。質量行列の集中化とデータの再構築法については、後で説明する。ステップ12で、同時緩和法の繰返計算の初期値の設定を行なう。ステップ13で、反復子ｋを１つ進める。予めｋを−１に初期化しておくので、最初のｋは０となる。ステップ14で、Newton-Raphson法によりPoisson方程式
∇2φ(k)＝∇・ｖ(k)
を解いて、修正速度ポテンシャルφ(k)を求める。ステップ15で、修正速度ポテンシャルφ(k)から、速度ｖ(k+1)と圧力ｐ(k+1)を求める。ステップ16で、速度の発散が十分小さくなったかどうか調べる。小さくなければ、収束していないと判断して、ステップ13に戻る。収束していれば、ステップ17で、速度ｖ(k+1)と圧力ｐ(k+1)を、時刻ｔ＝（ｎ＋１）Δｔの値ｖn+1,ｐn+1とする。
【００４３】
運動方程式の時間進行法について説明する。保存形の運動方程式
（∂ｖ／∂ｔ）＋∇・（ｖｖ）＝−∇ｐ＋∇・τ
を、速度を陽的に、圧力を陰的に離散化すると、
（ｖn+1−ｖn）／Δｔ＋∇・(ｖｖ)n＝−∇ｐn+1＋∇・τn
となる。肩付きのｎは、時刻を示す添え字である。ｖnは、時刻ｎΔｔにおける速度ベクトルである。時刻ｎΔｔにおける速度ｖnと圧力ｐnと粘性応力テンソルτnとから、時刻（ｎ＋１）Δｔにおける速度ｖn+1と圧力ｐn+1を計算する。その際、質量行列の集中化とデータの再構築法を用いて速度の予測子ｖ~を求め、修正速度ポテンシャルφに関するPoisson方程式を解いて、時刻（ｎ＋１）Δｔにおける値を求める。
【００４４】
速度の予測子について説明する。速度の予測子ｖ~は、Euler前進法により、方程式
（ｖ~−ｖn）／Δｔ＋∇・(ｖｖ)n＝−∇ｐn＋∇・τn
から求める。速度の予測子ｖ~は、時刻（ｎ＋１）Δｔにおける速度ｖn+1の第１近似である。時刻ｎΔｔにおける速度ベクトルｖnと、圧力ｐnと、粘性応力テンソルτnとから計算する。この計算には、質量行列の集中化とデータの再構築法を用いる。速度の予測子ｖ~は、発散０の条件を満たさないので、
∇・ｖ~≠０
となる。
【００４５】
修正速度ポテンシャルのφに関するPoisson方程式の導出方法を説明する。離散化した運動方程式と速度の予測子の方程式との差を作ると、
（ｖn+1−ｖ~）／Δｔ＝−∇（ｐn+1−ｐn）
となる。ここで、
ｖ'＝−(ｖn+1−ｖ~)
とすると、
∇×ｖ'＝Δｔ(∇×∇)(ｐn+1−ｐn)＝０
であるから、スカラーポテンシャルφが存在して、
∇φ＝−ｖ'
となる。φを修正速度ポテンシャルとよぶ。修正速度ポテンシャルφは、
∇φ／Δｔ＝∇（ｐn+1−ｐn）
であるから、
φ＝(ｐn+1−ｐn)Δｔ
となる。ところが、
∇・ｖn+1＝０
であるから、
∇・（ｖn+1−ｖ~）＝∇・（−ｖ~）＝−∇2φ
となる。これで、修正速度ポテンシャルφに関するPoisson方程式
∇2φ＝∇・ｖ~
が導かれた。
【００４６】
このPoisson方程式
∇2φ＝∇・ｖ~
の解は、Newton-Raphson法により、
φ＝−λ-1Ｄ~
と求まる。ただし、
λ＝(１／Ｖe)∬(１／Δｎ)dＳ（積分範囲はＳe）
Ｄ~＝(１／Ｖe)∬(ｖ~・ｎ)dＳ（積分範囲はＳe）
である。
【００４７】
同時緩和法について説明する。時刻ｎΔｔまでの状態が求まっているとする。
反復計算の初期値（反復子ｋ＝０のときの値）を、
ｖ(0)＝ｖ~
ｐ(0)＝ｐn
とする。φ(0)に関するPoisson方程式
∇2φ(0)＝∇・ｖ(0)
を、Newton-Raphson法により解き、
λ＝(１／Ｖe)∬(１／Δｎ)dＳ（積分範囲はＳe）
Ｄ(0)＝(１／Ｖe)∬(ｖ(0)・ｎ)dＳ（積分範囲はＳe）
φ(0)＝−λ-1Ｄ(0)
と求まて、φ(0)を得る。φ(0)から、
ｖ(1)＝ｖ(0)−∇φ(0)
ｐ(1)＝ｐ(0)＋φ(0)／Δｔ
を計算する。以下、同様にして、
∇2φ(k)＝∇・ｖ(k)
を解いてφ(k)を求め、
ｖ(k+1)＝ｖ(k)−∇φ(k)
ｐ(k+1)＝ｐ(k)＋φ(k)／Δｔ
を求める計算を繰り返す。すべての要素で、
｜∇・ｖ(k+1)｜≦εv（例えば、εv＝0.001）
を満足したとき、
ｖn+1＝ｖ(k+1)
ｐn+1＝ｐ(k+1)＝ｐ(k)＋φ(k)／Δｔ
として、つぎの時間ステップに進む。
【００４８】
ここで、質量の集中化とデータの再構築法について説明する。データの再構築の説明の前に、場の関数の時空平均について説明しておく。時空４次元空間の平均操作を、
φ*cn≡(１／ＶΔｔ)∬φdＶdｔ
（ｔの積分範囲は−Δｔ/2〜Δｔ/2、Ｖの積分範囲はＶ）
で定義し、φ*cnを、体積Ｖの図心点に割り当てる。これを、データの再構築と呼ぶ。時空平均は、時間平均φ*cと空間平均φ*nに分割できる。時間平均と空間平均は、順序が交換可能で、どちらから先に平均してもよい。また、平均操作を２回した場合の値は、１回の平均操作の値に等しい。
【００４９】
質量平均は、重みを密度ρに選んだ平均で、
φ*g≡(∫ρφdΩ)／（∫ρdΩ）  on Ω
で定義される。積分領域は、領域Ωである。そして、重心に、質量平均φ*gを割り当てる。
【００５０】
要素平均は、重み関数が１の平均で、
φ*c≡(∫φdΩ)／（∫dΩ）  on Ωe
で定義される。積分領域は、Ωeである。要素平均は、混同の恐れがない限り、上付き星印（＊）を省略して、φeと書く。要素平均は、空間平均と同じ意味を持つので、要素平均値φ*eを、要素の図心に割り当てる。
【００５１】
節点平均は、形状関数Ｎiを重みとした平均で、
φ*i≡<φ>i≡(∫ＮiφdΩ)／（∫ＮidΩ）  on Ωe
で定義する。積分領域は、形状関数Ｎiの台である。この場合は、要素の大きさを、体積加重平均したものである。そして、平均値φ*iは、節点ｉに割り当てられる。よって、上付き星印を省略して、φiと書いてよい。
【００５２】
質量の集中化について説明する。質量の集中化とは、重心に全質量を集めることを意味する。データの再構築とは、φ*gを、重心に割り当てることを意味する。特に、密度ρが一定の場合は、質量平均と空間平均は一致するから、図心に全体積を集め、φ*cを、図心に割り当てることと同じである。
【００５３】
離散系について、質量の集中化とデータの再構築を説明する。要素Ｖiを考える。節点iには、ｎ個のサブドメインからなる要素Ｖk（ｋ＝１，２，…，ｎ）が隣接している。要素Ｖkの質量をｍk、密度をρkとすれば、平均値の意味で、
ｍk＝ρkＶk
が成立し、要素Ｖiの全質量Ｍiと全体積Ｖiは、
Ｍi＝∫dｍ＝ｍ1＋ｍ2＋・・・＋ｍn  （全質量）
Ｖi＝∫dＶ＝Ｖ1＋Ｖ2＋・・・＋Ｖn  （全体積）
と表せる。このとき、要素Ｖiの質量平均は、
φ*g＝(∫ρφdＶ)／(∫ρdＶ)  on Ｖi
＝(ｍ1φ1＋ｍ2φ2＋・・・＋ｍnφn)／(ｍ1＋ｍ2＋・・・＋ｍn)
＝(ρ1Ｖ1φ1＋ρ2Ｖ2φ2＋・・・＋ρnＶnφn)
／(ρ1Ｖ1＋ρ2Ｖ2＋・・・＋ρnＶn)
と表せる。ここで、φk（ｋ＝１，２，…，ｎ）は、要素Ｖkのφの平均値である。特に、密度が一定で、
ρ1＝ρ2＝・・・＝ρn＝一定
ならば、φ*gは、φ*cに一致する。ただし、
φ*c＝(∫φdＶ)／(∫dＶ)
＝(Ｖ1φ1＋Ｖ2φ2＋・・・＋Ｖnφn)／(Ｖ1＋Ｖ2＋・・・＋Ｖn)
である。このφ*gと、φ*cを、重心と図心に割り当てるのが、データの再構築である。すなわち、連続要素Ｖiをｎ個のサブ要素に分割し、サブ要素を質点系と考えて、全体を再構築する。
【００５４】
サイクル誤差自己調整法を採用する場合、保存則を表現する原理式を用いて離散化する。原理式には、Lagrange表示とEuler表示の２通りの表示方法がある。Lagrange表示は、物質検査体積Ｖm(t)を用いるものであり、Euler表示は、固定検査体積Ｖを用いるものである。物質検査体積を用いるLagrange表示の場合は、その検査体積の重心点が重要な意味を持つ。よって、時間進行法は、この重心点に関して定式化される。一方、Euler表示においては、図心点が重要になり、時間進行法は、図心点に関して定式化される。密度ρが一定の非圧縮性流体においては、重心点と図心点は一致する。
【００５５】
重心と図心の定義
ｒg≡(∫ρｒdΩ)／(∫ρdΩ)  （重心）
ｒc≡(∫ｒdΩ)／(∫dΩ)      （図心）
【００５６】
一方、有限要素法の質量集中化は、節点で実行される。節点は、重心点とも図心点とも異なっている。よって、大小さまざまな要素を組み合わせた場合、質量の集中化は誤差を発生し、非物理的となる。ゆえに、重心または図心にデータを再構築し、発生した誤差を除去し、物理的に正しい時間進行を行なう必要がある。
【００５７】
Euler的時間進行法について説明する。固定検査体積をＶとする。このとき、質量保存則と運動量のつり合い方程式は、
(∂／∂ｔ)(∫ρdＶ)  (積分範囲はＶ)
＋∫∇・(ρｖｖ)dＶ＝０  (積分範囲はＶ)
(∂／∂ｔ)(∫ρｖdＶ)  (積分範囲はＶ)
＋∫∇・(ρｖｖ)dＶ  (積分範囲はＶ)
＝(∫∇・ＴdＶ)＋(∫ρｂdＶ)  (積分範囲はＶ)
と書ける。ここでも、体積Ｖの重心をｒgとする。体積Ｖ内の全質量は、時間の関数として、
Ｍ(t)＝∫ρdＶ  (積分範囲はＶ)
と書ける。また、
ｖ*g≡∫ρｖdＶ／Ｍ  （積分範囲はＶ）
ｂ*g≡∫ρｂdＶ／Ｍ  （積分範囲はＶ）
Π≡ρｖｖ−Ｔ：運動量テンソル
とする。ここで、運動量のつり合い方程式の時間積分を実行すると、
(Ｍn+1ｖ*gn+1−Ｍnｖ*gn)／Δｔ＋(∇・Π)*gn+1/2
＝Ｍt+1/2ｂ*gn+1/2
となる。この場合、
Ｍn+1≠Ｍn
であるから、このような分離は、時間進行に適さない。非圧縮で、ρ＝const.ならば、
Ｍn+1＝Ｍn+1/2＝Ｍn＝Ｍ＝const.
であるから、両辺をＭで割り、時間進行可能となる。
【００５８】
非圧縮流体場を有限要素法で定式化する場合、質量保存則に対しては、重み関数として１を用いる。各要素で発散の要素平均値を計算する。すなわち、
θe＝<∇・ｖ>e
＝(１／Ωe(ｘ))∫∇・ｖdΩ(ｘ)  （積分範囲はΩe）
＝∇a・ｖa
とする。ここで、∇aは、離散ナブラ演算子である。また、
ｖ＝ΣＮaｖa＝Ｎaｖa
で、Ｎaは形状関数である。
【００５９】
運動方程式のEuler前進にともなう半離散化式は、速度を陽的、圧力を陰的に離散化して、
(ｖn+1−ｖn)／Δｔ＋∇・(ｖｖ)en
＝−∇ｐen+1＋∇・τen＋∇・τｖenｖen・∇ｖn
となる。ただし、拡散項に対しては、安定化のためアワーグラス(hourglass)行列を付加する。最後の項は上流化項で、BTD項と呼ぶ。τは、安定化に伴う時間パラメータで、普通、△ｔ／２である。最適化パラメータτoptその他を採用することもできる。
【００６０】
上記の運動方程式の節点平均操作< >iによって得られる節点速度ｖin+1は、物理的に図心速度ｖcn+1を意味する。すなわち、図心点ｃにデータを再構築する必要がある。大小の要素が混在する場合、このデータの再構築が大きな意味を持つ。
【００６１】
データの再構築法について説明する。データの再構築法は、Taylor展開によって実行される。関数ｆ(ｒ)をTaylor展開すると、
ｆ(ｒ＋Δｒ)＝ｆ(ｒ)＋Δｒ・∇ｆ＋(１／２)(Δｒ・∇)2ｆ＋・・・
と書ける。ここで、重心をｒg、節点をｒiとし、
ｒ＝ｒg
ｒ＋Δｒ＝ｒi
Δｒ＝ｒi−ｒg＝(Δｒ)g
と置く。すなわち、重心まわりのTaylor展開は、２次以上を省略して、
ｆ(ｒi)＝ｆ(ｒg)＋(Δｒ・∇)gｆ＋(１／２)(Δｒ・∇)g2ｆ
となる。右辺第２項と第３項が補正項である。同様に、図心をｒc、節点をｒiとし、
ｒ＝ｒc
ｒ＋Δｒ＝ｒi
Δｒ＝ｒi−ｒc＝(Δｒ)c
と置けば、図心に対する補正式は、
ｆ(ｒi)＝ｆ(ｒc)＋(Δｒ・∇)cｆ＋(１／２)(Δｒ・∇)c2ｆ
と求まる。ここで、関数ｆの１階微分と２階微分の重心と図心での値が必要になるが、誤差の２乗が最小になるように決定するのがよい。
【００６２】
図３を参照しながら、数値計算装置の動作を説明する。最初に、パラメータと境界条件と初期条件とメッシュデータを、図３に示す入力装置１から入力する。
入力したデータは、記憶装置２に格納しておく。時間進行を行なう際に、節点データを重心データまたは図心データに変換する。重心で時間進行する場合では、データ再構築手段７により節点のデータを重心のデータに変換する。図心で時間進行する場合は、節点のデータを図心のデータに変換する。
【００６３】
初期条件を使って、ｔ＝Δｔの予測子、すなわち、ｎ＝１の速度の予測子ｖ~の計算を行なう。速度の予測子ｖ~と圧力ｐ0を初期値として、同時緩和法の繰返計算を行なう。すなわち、Poisson方程式
∇2φ(k)＝∇・ｖ(k)
を、Newton-Raphson法演算手段５で、Newton-Raphson法により解いて、修正速度ポテンシャルφ(k)を得る。時間進行手段６により、修正速度ポテンシャルφ(k)から、次の近似値ｖ(k+1),ｐ(k+1)を求める。これを、収束するまで繰り返す。収束したら、ｖ(k+1),ｐ(k+1)を、データ再構築手段７により、重心または図心のデータから節点のデータを再構成し、新しい状態ｖn+1,ｐn+1とする。この計算結果を出力装置３に出力する。この計算を、演算制御手段４の制御により、所定回数（Ｎ回）繰り返す。このようにして、Navier-Stokesの方程式を解くことができる。
【００６４】
上記のように、本発明の第１の実施の形態では、数値計算方法を、ソレノイダル・ベクトル場の保存形Navier-Stokes運動方程式を、速度を陽的に、圧力を陰的に離散化し、サイクル誤差自己調整法を用いた同時緩和法で、速度の予測子を、質量行列の集中化とデータの再構築法を用いて、Euler前進法により求め、速度の発散が０という条件の修正速度ポテンシャルに関するPoisson方程式をNewton-Raphson法で解き、修正速度ポテンシャルから速度と圧力を求めることを、収束するまで繰り返して、次の時間ステップの速度と圧力を求める時間進行を行なう構成としたので、計算精度が上がり、解が安定に収束する。
【００６５】
（第２の実施の形態）
本発明の第２の実施の形態は、磁場の誘導方程式を離散化し、磁場の予測子をデータの再構築法を用いて求め、修正磁場ポテンシャルに関するPoisson方程式を解いて磁場を求めるサイクル誤差自己調整法を用いた同時緩和法で時間進行を行なう数値計算方法である。
【００６６】
図４は、本発明の第２の実施の形態における数値計算方法の時間進行の手順を示す流れ図である。基本的な演算手順は、図１に示した第１の実施の形態と同じである。数値計算装置の基本的な構成も、図３に示した第１の実施の形態と同じである。
【００６７】
上記のように構成された本発明の第２の実施の形態における数値計算方法の処理手順を説明する。最初に、修正磁場ポテンシャルに関するPoisson方程式の導出を説明する。磁場の誘導方程式を、
∂Ｂ／∂ｔ＋ｖ・∇Ｂ＝−∇Ｒ＋Ｂ・∇ｖ＋γ∇2Ｂ
と変形する。ここで、
Ｒ＝γ∇・Ｂ
は、理論的には０であるが、数値残差として残しておく。
速度ｖとＲを陰的に定式化すると、
(Ｂn+1−Ｂn)／Δｔ＋ｖen+1・∇Ｂn
＝−∇Ｒen+1＋Ｂen・∇ｖn+1＋γ∇2Ｂn＋BTDn
となる。ここで、
BTD＝∇・Ｋ・∇Ｂ
である。
【００６８】
前進Euler法で、磁束密度の予測値を求めると、
(Ｂ~−Ｂn)／Δｔ＋ｖen+1・∇Ｂn
＝−∇Ｒen＋Ｂen・∇ｖn+1＋γ∇2Ｂ＋BTDn
となる。もとの離散式から、この予測子の式を差し引くと、差式
(Ｂn+1−Ｂ~)／Δｔ＝−∇(Ｒen+1−Ｒen)
が導かれる。修正磁場ポテンシャルψを、
Ｂn+1＝Ｂ~−∇ψ
で定義する。
∇・Ｂn+1＝０
の条件より、修正磁場ポテンシャルψに関するPoisson方程式
∇2ψ＝∇・Ｂ~
が導かれる。
【００６９】
方程式の解法の手順を説明する。方程式
∂Ｂ／∂ｔ＋ｖ・∇Ｂ＝−∇Ｒ＋Ｂ・∇ｖ＋γ∇2Ｂ
を、
∇・Ｂ＝０
という条件のもとで、磁束密度Ｂと数値残差Ｒの同時緩和法で解く。解法は、運動方程式と同様である。まず、前進Euler法で、磁場の予測値Ｂ~を求める。速度ｖen+1はすでに求められているので、
(Ｂ~−Ｂn)／Δｔ＋ｖen+1・∇Ｂn
＝−∇Ｒen＋Ｂen・∇ｖn+1＋γ∇2Ｂ＋BTDn
から、予測値Ｂ~を求める。このとき、質量行列の集中化とデータの再構築法を用いる。
【００７０】
反復子ｋが０の場合の初期値を、
Ｂ(0)＝Ｂ~、Ｒ(0)＝Ｒn
とする。修正磁場ポテンシャルψ(k)に関するPoisson方程式
∇2ψ(k)＝∇・Ｂ(k)
を、Newton-Raphson法により、
∇・Ｂn+1＝０
を満足させながら解く。つぎの漸化式
Ｂ(k+1)＝Ｂ(k)−∇ψ(k)
Ｒ(k+1)＝Ｒ(k)＋ψ(k)／Δｔ
で、近似を進める。この計算を繰り返す。反復の結果、すべての要素で、
｜∇・Ｂ(k+1)｜≦εB（例えば、εB＝0.001）
を満足したとき、
Ｂn+1＝Ｂ(k+1)、Ｒn+1＝Ｒ(k+1)
として、つぎの時間ステップに進む。
【００７１】
図１と図３と図４を参照しながら、数値計算装置の動作を説明する。最初に、パラメータと境界条件と初期条件とメッシュデータを、図３に示す入力装置１から入力する。入力したデータは、記憶装置２に格納しておく。節点データを重心データまたは図心データに変換する。重心で時間進行する場合では、データ再構築手段７により節点のデータを重心のデータに変換する。図心で時間進行する場合は、節点のデータを図心のデータに変換する。
【００７２】
初期条件を使って、ｔ＝Δｔの予測子、すなわち、ｎ＝１の磁場の予測子Ｂ~の計算を行なう。磁場の予測子Ｂ~とＲ0を初期値Ｂ(0)とＲ(0)として、同時緩和法の繰返計算を行う。すなわち、Poisson方程式
∇2ψ(0)＝∇・Ｂ(0)
をNewton-Raphson法演算手段５でNewton-Raphson法により解く。時間進行手段６により、Newton-Raphson法の解から修正磁場ポテンシャルψを計算し、修正磁場ポテンシャルψから、次の近似値Ｂ(k+1),Ｒ(k+1)を求める。これを、収束するまで繰り返す。収束したら、Ｂ(k+1),Ｒ(k+1)を新しい状態Ｂ(n+1),Ｒ(n+1)とする。
【００７３】
データ再構築手段７により、重心または図心のデータから節点のデータを再構成する。計算結果を出力装置３に出力する。演算制御手段４の制御により、所定回数（Ｎ回）繰り返す。このようにして、Navier-Stokesの方程式を解くことができる。
【００７４】
上記のように、本発明の第２の実施の形態では、数値計算方法を、磁場の誘導方程式を離散化し、磁場の予測子をデータの再構築法を用いて求め、修正磁場ポテンシャルに関するPoisson方程式を解いて磁場を求めるサイクル誤差自己調整法を用いた同時緩和法で時間進行を行なう構成としたので、計算精度が上がり、解が安定に収束する。
【００７５】
（第３の実施の形態）
本発明の第３の実施の形態は、定常電磁場の電位に関するPoisson方程式を離散化し、電流場と電場の予測子を、データの再構築法を用いて求め、サイクル誤差自己調整法を用いた同時緩和法で、修正電位ポテンシャルに関するPoisson方程式を解き、電流場と電場を求める数値計算方法である。
【００７６】
図５は、本発明の第３の実施の形態の数値計算方法の概略手順を示す流れ図である。図６は、同時緩和処理の手順を示す流れ図である。図７は、数値計算方法を実行する数値計算装置の機能ブロック図である。図７において、入力装置１は、パラメータなどを入力するキーボードなどの手段である。記憶装置２は、パラメータや境界条件などを記憶するメモリである。出力装置３は、計算結果を出力する表示装置などの出力手段である。演算制御手段４、全体の演算動作を制御する手段である。Newton-Raphson法演算手段５は、Poisson方程式を解く手段である。データ再構築手段７は、節点のデータから重心または図心のデータを求めるとともに、重心または図心のデータから節点のデータを再構成する手段である。
同時緩和手段８は、繰返演算を行なう手段である。
【００７７】
上記のように構成された本発明の第３の実施の形態における数値計算方法の処理手順を説明する。最初に、電位に関するPoisson方程式を説明する。直流電流の場合、誘起磁場は存在しないから電場の回転は０となる。よって、電場Ｅは、電位Ｖの勾配として、
Ｅ＝−∇Ｖ
と表せる。このとき、Ohmの法則の発散をとると、電位に関するPoisson方程式
∇2Ｖ＝∇・(ｖ×Ｂ)
が導かれる。この方程式は、
∇・ｊ＝０
の条件で解かれる。Poisson方程式を解くことを目的とせず、物理的に重要な、発散が０となるソレノイダル性を満足する電流場を求めることに重点を置く。
【００７８】
電位に関するPoisson方程式
∇2Ｖ＝∇・(ｖ×Ｂ)
を近似的に解いて、電位の予測子Ｖ~を求める。つぎに、Ohmの法則
ｊ~＝−∇Ｖ~＋ｖ×Ｂ
により、電流の予測子ｊ~を求める。この予測子ｊ~は当然、
∇・ｊ~＝０
を満足しない。そこで、修正電位ポテンシャルχを導入して、修正電位ポテンシャルχを、次の反復法により求める。
【００７９】
すなわち、の初期値（反復子ｋ＝０の値）を、
ｊ(0)＝ｊ~，Ｖ(0)＝Ｖ~
とする。χ(k)に関するPoisson方程式
∇2χ(k)＝∇・ｊ(k)
を、Newton-Raphson法により解く。得られたχ(k)から、
ｊ(k+1)＝ｊ(k)−∇χ(k)
Ｖ(k+1)＝Ｖ(k)＋χ(k)
を計算する。この計算を繰り返した結果、すべての要素で、
｜∇・ｊ(k+1)｜≦εj（例えば、εj＝0.001）
を満足したとき、
ｊ＝ｊ(k+1)、Ｖ＝Ｖ(k+1)
として、反復を終了する。このようにして、ソレノイダル性を満足する電流場を求める。
【００８０】
図５と図６と図７を参照しながら、数値計算装置の動作を説明する。ステップ31で、Poisson方程式
∇2Ｖ＝∇・(ｖ×Ｂ)
を離散化する。離散化したPoisson方程式のパラメータを、図７の記憶装置２に、境界条件とメッシュデータとともに格納する。ここまでは準備段階である。ステップ32で、電位の予測子Ｖ~と電流の予測子ｊ~を求める。ステップ33で、予測子を初期値として、同時緩和法の計算を行なう。その詳細を、図６に示す。
【００８１】
図６に示す同時緩和処理のステップ41で、初期値を設定する。ステップ42で、反復子ｋを１つ進める。ｋの初期値は−１である。ステップ43で、Newton-Raphson法演算手段５により、修正電位ポテンシャルχ(k)を求める。ステップ44で、修正電位ポテンシャルχ(k)から、電流ｊ(k+1)と電位Ｖ(k+1)を求める。ステップ45で、収束を検査する。収束していなければ、ステップ42から繰り返す。収束していれば、電流ｊ(k+1)と電位Ｖ(k+1)の値を解ｊ、Ｖとする。
【００８２】
上記のように、本発明の第３の実施の形態では、数値計算方法を、定常電磁場の電位に関するPoisson方程式を離散化し、電流場と電場の予測子を、データの再構築法を用いて求め、サイクル誤差自己調整法を用いた同時緩和法で、修正電位ポテンシャルに関するPoisson方程式を解き、電流場と電場を求める構成としたので、定常状態の計算でも計算精度が上がり、解が安定に収束する。
【００８３】
【発明の効果】
以上の説明から明らかなように、本発明では、ソレノイダル・ベクトル場の保存形Navier-Stokes運動方程式を、速度を陽的に圧力を陰的に離散化し、サイクル誤差自己調整法を用いた同時緩和法で時間進行を行なって次の時間ステップの速度と圧力を求める有限要素法の数値計算方法の時間進行の処理で、速度の予測子を質量行列の集中化とデータの再構築法を用いてEuler前進法により求め、速度の発散が０という条件の修正速度ポテンシャルに関するPoisson方程式をNewton-Raphson法で解いて得た修正速度ポテンシャルから速度と圧力を求める計算を収束するまで繰り返す構成としたので、ソレノイダル・ベクトル場のNavier-Stokesの方程式を有限要素法で解く場合に、計算精度が上がるとともに、解が安定に収束するようになるという効果が得られる。
【図面の簡単な説明】
【図１】本発明の第１の実施の形態における数値計算方法の概略の流れ図、
【図２】本発明の第１の実施の形態における数値計算方法の時間進行の流れ図、
【図３】本発明の第１の実施の形態における数値計算方法を実行する数値計算装置の機能ブロック図、
【図４】本発明の第２の実施の形態における数値計算方法の時間進行の流れ図、
【図５】本発明の第３の実施の形態における数値計算方法の概略の流れ図、
【図６】本発明の第３の実施の形態における数値計算方法の同時緩和の流れ図、
【図７】本発明の第３の実施の形態における数値計算方法を実行する数値計算装置の機能ブロック図、
【図８】要素の面に応じた積分∬(１／Δｎ)dＳの評価方法を説明する図、
【図９】∬(１／Δｎ)dＳの積分方法を説明する図、
【図１０】要素Ｋと要素Ｌの境界面上での∂φ／∂ｎの計算方法を説明する図、
【図１１】圧力ｐと速度ｖの同時緩和法の手順を示すPAD図、
【図１２】数値残差Ｒと磁束密度Ｂの同時緩和法の手順を示すPAD図である。
【符号の説明】
１  入力装置
２  記憶装置
３  出力装置
４  演算制御手段
５  Newton-Raphson法演算手段
６  時間進行手段
７  データ再構築手段
８  同時緩和手段
【図１】
0004649673.tif 000005
【図２】
0004649673.tif 000006
【図３】
0004649673.tif 000007
【図４】
0004649673.tif 000008
【図５】
0004649673.tif 000009
【図６】
0004649673.tif 000010
【図７】
0004649673.tif 000011
【図８】
0004649673.tif 000012
【図９】
0004649673.tif 000013
【図１０】
0004649673.tif 000014
【図１１】
0004649673.tif 000015
【図１２】
0004649673.tif 000016