(57)【要約】【課題】電磁場・温度場・流体場が同時に解けて、高速高精度で汎用性が高い解析システムを実現する。【解決手段】GSMAC-FEM流体解析手段で、電磁力にMaxwellの応力テンソルの発散型を用いた運動方程式と連続の方程式に高Reynolds数と高Hartmann数に安定なGSMAC有限要素を適用して流体場を計算する。出力データを、他の解析手段への入力データに変換する。CIP-FEM熱流体解析手段で、エネルギー方程式に高Prandtl数に安定なCIP有限要素を適用して温度場の影響を計算する。出力データを、他の解析手段への入力データに変換する。A-φFEM電磁流体解析手段で、高磁気Reynolds数に安定な方程式に辺有限要素を適用して電磁場の影響を計算する。出力データを、他の解析手段への入力データに変換する。このように、各解析手段を連続して順次動作させてそれぞれの解を求め、電磁場・温度場・流体場を同時に解く。【選択図】図１
2012074016.tif 000002
【特許請求の範囲】
【請求項１】
構造データとパラメータと初期値を入力する入力手段と、入力した構造データとパラメータと初期値を格納する記憶装置と、電磁力にMaxwellの応力テンソルの発散型を用いた運動方程式と連続の方程式に高Reynolds数と高Hartmann数に安定なGSMAC有限要素法を適用して流体場を計算するGSMAC-FEM流体解析手段と、エネルギー方程式に高Prandtl数に安定なCIP有限要素法を適用して温度場を計算するCIP-FEM熱流体解析手段と、高磁気Reynolds数に安定な方程式に辺有限要素法を適用して電磁場を計算するA-φFEM電磁流体解析手段と、GSMAC-FEM流体解析手段の出力データをCIP-FEM熱流体解析手段とA-φFEM電磁流体解析手段への入力データに変換する第１データ変換手段と、CIP-FEM熱流体解析手段の出力データをA-φFEM電磁流体解析手段とGSMAC-FEM流体解析手段への入力データに変換する第２データ変換手段と、A-φFEM電磁流体解析手段の出力データをGSMAC-FEM流体解析手段とCIP-FEM熱流体解析手段への入力データに変換する第３データ変換手段と、ある時刻においてGSMAC-FEM流体解析手段とCIP-FEM熱流体解析手段とA-φFEM電磁流体解析手段を連続して順次動作させてそれぞれの解を求めるように制御する計算制御手段とを具備することを特徴とする電磁熱流体解析装置。
【請求項２】
流体場・温度場・電磁場が密接に関連している電磁熱流体の運動をコンピュータにより数値解析する電磁熱流体解析方法であって、入力手段を介して入力された構造データと初期値とパラメータとを記憶手段に格納し、電磁力にMaxwellの応力テンソルの発散型を用いた運動方程式と連続の方程式に高Reynolds数と高Hartmann数に安定なGSMAC有限要素を適用して流体を解析するGSMAC-FEM流体解析手段で流体場を計算し、GSMAC-FEM流体解析手段の出力データをCIP-FEM熱流体解析手段とA-φFEM電磁流体解析手段への入力データに第１データ変換手段で変換し、エネルギー方程式に高Prandtl数に安定なCIP有限要素を適用して流体の温度場を解析するCIP-FEM熱流体解析手段で温度場の影響を計算し、CIP-FEM熱流体解析手段の出力データをA-φFEM電磁流体解析手段とGSMAC-FEM流体解析手段への入力データに第２データ変換手段で変換し、高磁気Reynolds数に安定な方程式に辺有限要素を適用して流体の電磁場を解析するA-φFEM電磁流体解析手段で電磁場の影響を計算し、A-φFEM電磁流体解析手段の出力データをGSMAC-FEM流体解析手段とCIP-FEM熱流体解析手段への入力データに第３データ変換手段で変換し、ある時刻においてGSMAC-FEM流体解析手段とCIP-FEM熱流体解析手段とA-φFEM電磁流体解析手段を連続して順次実施してそれぞれの解を求めることにより電磁熱流体の運動を解析することを特徴とする電磁熱流体解析方法。
【請求項３】
請求項２に記載の電磁熱流体解析方法をコンピュータに実行させるための処理手順を記述したコンピュータプログラム。
【請求項４】
請求項３記載のコンピュータプログラムを記録したコンピュータ読取可能な記録媒体。
【発明の詳細な説明】
【技術分野】
【０００１】
本発明は、有限要素法を用いる電磁熱流体解析装置に関し、特に、高磁気Reynolds数で高Prandtl数の電磁熱流体でも安定に解析できる汎用性のある電磁熱流体解析装置に関する。
【背景技術】
【０００２】
従来、有限要素法による流体解析システムとしては、圧力などの力学的な外場のみを考慮した流体解析システムと、電磁場を含む外場を考慮した電磁流体解析システムと、熱的な外場を考慮した熱流体解析システムがある。解析対象としては、特許文献１に開示された「ガラス溶融炉」などがある。以下に、これに関連する従来技術の例をいくつかあげる。
【０００３】
特許文献２に開示された「離散ナブラ演算子と非記憶型係数行列による数値解析装置」は、記憶容量を少なく計算精度を維持しながらも、計算を高速化する非構造格子の有限体積法または有限要素法により計算する数値解析装置である。図10（ａ）に示すように構成されている。nabla演算子∇がベクトル量であり、座標系に依存しない量であることに着目し、離散nabla演算子を定義する。離散nabla演算子で記述した係数行列要素を、要素形状に依存する部分と、要素形状に依存しない部分の積に分解する。要素形状に依存する部分のみを記憶する。他の部分は、解析的に求めておく。行列計算の際に必要になるたびに、これらから係数行列要素を生成する。非構造格子の有限体積法または有限要素法において、この離散nabla演算子を用いると、非記憶型の係数行列を作ることができる。したがって、係数行列要素をすべて記憶することなく係数行列を解くことができ、記憶装置の低容量化と計算の高速化ができる。
【０００４】
特許文献３に開示された「質量集中化とデータの再構築および双対格子による数値計算方法」は、有限要素法によるNavier-Stokes方程式の数値計算において、計算を安定化し精度を上げる方法である。図10（ｂ）に示すように構成されている。Navier-Stokes方程式の積分形の原理式を、双対格子による双対要素を用いて要素平均と節点平均を組み合わせて、時空４次元空間での平均操作により離散化する。構成方程式には要素平均を適用し、運動方程式には節点平均を適用して離散化する。質量を要素の重心に集中化して、要素の重心点または図心点にデータを再構築して、node by element法で時間進行する。このようにして、質量集中化による誤差を除いて安定に高精度で計算できる。
【０００５】
特許文献４に開示された「サイクル誤差自己調整法を用いた同時緩和法による有限要素法の数値計算方法」は、ソレノイダル・ベクトル場のNavier-Stokesの方程式を有限要素法で解く場合に、計算精度を上げるとともに解の収束を高速化かつ安定化する方法である。図10（ｃ）に示すように計算する。ソレノイダル・ベクトル場の保存形Navier-Stokes運動方程式を、速度を陽的に、圧力を陰的に離散化する。サイクル誤差自己調整法を用いた同時緩和法で時間進行を行って、次の時間ステップの速度と圧力を求める。時間進行の処理では、速度の予測子を、質量行列の集中化とデータの再構築法を用いて、Euler前進法により求める。速度の発散が０という条件の修正速度ポテンシャルに関するPoisson方程式をNewton-Raphson法で解き、修正速度ポテンシャルから速度と圧力を求めることを、収束するまで繰り返して、次の時間ステップの速度と圧力を求める。
【０００６】
特許文献５に開示された「連続体の数値解析方法」は、誘導磁場下でローレンツ力およびジュール熱を受ける連続体の物理量を少ない計算時間で数値解析する方法である。図10（ｄ）に示すように計算する。誘導磁場下で運動する連続体に生じるローレンツ力およびジュール熱が連続体の速度により移流する移流方程式を設定し、連続体の運動量保存式は毎時間段階で解き、マクスウェル方程式はとびとびの時間段階において解き、このとびとびの時間段階の間の時間段階において移流方程式を解くことにより、連続体内のローレンツ力およびジュール熱を時間的に更新する連続体の数値解析方法である。とびとびの時間段階の間の省略した時間段階におけるローレンツ力およびジュール熱の解析精度低下を最小限に食い止めながら、連続体の物理量についての計算時間を大幅に短縮することが可能となる。
【０００７】
特許文献６に開示された「鋳造シミュレーションによる最適鋳造条件策定方法」は、鋳造シミュレーションの精度を向上させ、空隙の発生を的確に予測できるようにすることで、鋳造シミュレーションにより常に最適な鋳造条件を容易に策定し得る方法である。図10（ｅ）に示すように、鋳型およびそのキャビティ領域を、座標系上に位置づけられた複数の微小セルからなる計算要素に分割する。鋳型と、鋳造時にそのキャビティ領域から流出する気体と、キャビティ領域に流入する溶融金属について、夫々物性値、初期条件、境界条件を入力し、各計算要素の経時的変動を算出する。
【０００８】
特許文献７に開示された「加熱装置の設定方法」は、目標となる温度分布を先に決定して、最初に流動解析を１回行うだけで、後は電磁場解析のみで理想的な発熱分布に近いコイルレイアウトと電圧設定を行える加熱装置の設計方法である。図10（ｆ）に示すように加熱する。誘導加熱により中空状の被加熱体を加熱する。被加熱体の中空内に流体が存在する加熱装置である。被加熱体の内表面の目標となる温度分布を設定する。温度分布を１つの条件として流動解析を行い、熱流束分布を算出する。熱流束分布のプロファイルを発熱密度分布のプロファイルに変換する。電磁場解析により発熱密度分布を算出する。前者の発熱密度分布と後者の発熱密度分布とを比較し、比較結果に応じて条件を変更し、誘導加熱の条件決定を行う。
【０００９】
非特許文献１において、「CIP有限要素法を用いた立方体キャビティ内自然対流解析」が報告されている。３次元問題の検証例として、立方体キャビティ内自然対流を解析したものである。図10（ｇ）に示すようにして計算する。立方体キャビティ内自然対流は、正方形キャビティ内自然対流と同様に、スキーム検証のベンチマーク問題として数多く研究が行われ、良好なデータが入手しやすい。ただし、３次元では、Ｒa≧107の高Rayleigh数流れは、依然として計算負荷が大きい。さらに、エネルギー方程式の対流項の強い非線形性が原因で、数値的な発散や振動が生じ、解析が困難になるため、解析例は少ない。３次元で103≦Ｒa≦108の流れを解析し、他者の解析結果と比較することで、CIP有限要素法の有効性を示した。
【００１０】
非特許文献２において、「電磁場と速度場の相互作用を考慮した新しい電磁流体解析手法の構築」が報告されている。速度起電力を考慮した、電磁場と速度場の相互作用を解析可能な手法の構築を行ったものである。図10（ｈ）に示すようにして計算する。電磁場の解法には、速度起電力を考慮したA法辺要素有限要素法を用いる。速度場の解法として、GSMAC有限要素法を用いる。この二つの解法を弱連成させることにより、電磁流体解析手法の構築を行う。３次元立方体cavity内MHD強制対流問題を解析できる。
【先行技術文献】
【特許文献】
【００１１】
【特許文献１】特開2002-328197号公報
【特許文献２】特開2002-328913号公報
【特許文献３】特開2002-328955号公報
【特許文献４】特開2003-067426号公報
【特許文献５】特開2005-182150号公報
【特許文献６】特開2005-246439号公報
【特許文献７】特開2007-059292号公報
【非特許文献】
【００１２】
【非特許文献１】槇原孝文、三好市朗、棚橋隆彦：「CIP有限要素法を用いた立方体キャビティ内自然対流解析」（論文No.00-0120）、日本機械学会論文集（Ｂ編）66巻650号(2000-10)、pp.2513-2520。
【非特許文献２】近藤陽介、棚橋隆彦：「電磁場と速度場の相互作用を考慮した新しい電磁流体解析手法の構築」、日本計算工学会論文集、NO.20060018、(2006-5-12)、pp.1-8。
【発明の開示】
【発明が解決しようとする課題】
【００１３】
しかし、上記従来の方法では、温度場と流体場を同時に解くことはできるし、電磁場と流体場も同時に解くことができるが、電磁場と温度場と流体場を同時に解くことはできない。実用的な問題の多くは、流体場・温度場・電磁場が密接に関連している。そこで、それぞれの数値解析システムを使って、工夫をこらして解析している。しかし、この３つの場を統一的に解ける単一の数値解析システムはない。
【００１４】
本発明の目的は、上記従来の問題を解決して、電磁場・温度場・流体場が同時に解けて、高速高精度で汎用性が高い解析システムを実現することである。
【課題を解決するための手段】
【００１５】
上記の課題を解決するために、本発明では、電磁熱流体解析装置を、構造データとパラメータと初期値を入力する入力手段と、入力した構造データとパラメータと初期値を格納する記憶装置と、電磁力にMaxwellの応力テンソルの発散型を用いた運動方程式と連続の方程式に高Reynolds数と高Hartmann数に安定なGSMAC有限要素を適用して流体場を計算するGSMAC-FEM流体解析手段と、エネルギー方程式に高Prandtl数に安定なCIP有限要素を適用して温度場の影響を計算するCIP-FEM熱流体解析手段と、高磁気Reynolds数に安定な方程式に辺有限要素を適用して電磁場の影響を計算するA-φFEM電磁流体解析手段と、GSMAC-FEM流体解析手段の出力データをCIP-FEM熱流体解析手段とA-φFEM電磁流体解析手段への入力データに変換する第１データ変換手段と、CIP-FEM熱流体解析手段の出力データをA-φFEM電磁流体解析手段とGSMAC-FEM流体解析手段への入力データに変換する第２データ変換手段と、A-φFEM電磁流体解析手段の出力データをGSMAC-FEM流体解析手段とCIP-FEM熱流体解析手段への入力データに変換する第３データ変換手段と、ある時刻においてGSMAC-FEM流体解析手段とCIP-FEM熱流体解析手段とA-φFEM電磁流体解析手段を連続して順次動作させてそれぞれの解を求めるように制御する計算制御手段とを具備する構成とした。
【発明の効果】
【００１６】
上記のように構成したことにより、電磁場・温度場・流体場が同時に解けて、高速高精度で汎用性が高い解析システムを実現でき、３次元でも非定常でも非構造格子でも、誘導磁場や熱対流などを考慮した電磁熱流体の流れが解析できる。
【図面の簡単な説明】
【００１７】
【図１】本発明の実施例における電磁熱流体解析装置の動作手順を示す流れ図である。
【図２】本発明の実施例における電磁熱流体解析装置の概念図である。
【図３】本発明の実施例における電磁熱流体解析装置の支配方程式を説明する図である。
【図４】本発明の実施例における電磁熱流体解析装置での接続行列による変数の変換と形状関数の変換を示す表と、関数空間の関係を示す表である。
【図５】本発明の実施例における電磁熱流体解析装置に関する双対空間と接続行列の関係を示す図である。
【図６】本発明の実施例における電磁熱流体解析装置に関する木・補木ゲージによる正則化を示す図と、A-V法の有限要素方程式を示す図である。
【図７】本発明の実施例における電磁熱流体解析装置に関する積分の説明図である。
【図８】本発明の実施例における電磁熱流体解析装置による計算例の説明図である。
【図９】本発明の実施例における電磁熱流体解析装置による計算例の説明図である。
【図１０】従来の流体解析装置の概念図である。
【発明を実施するための形態】
【００１８】
以下、本発明を実施するための最良の形態について、図１〜図９を参照しながら詳細に説明する。
【実施例】
【００１９】
本発明の実施例は、GSMAC-FEM流体解析手段で微小時間の流体の力学的変化を解析してデータ変換し、CIP-FEM熱流体解析手段で流体の熱力学的変化を解析してデータ変換し、A-φFEM電磁流体解析手段で流体の電磁気学的変化を解析してデータ変換し、これを順次実行して電磁熱流体の動きを解析する電磁熱流体解析装置である。
【００２０】
図１に、電磁熱流体解析装置の動作を示す流れ図を示す。図２に、電磁熱流体解析装置の概念図を示す。主なデータ経路のみを示してあり、当然必要な自明なデータ経路は省略してある。図３に、支配方程式を示す。図４に、電磁熱流体解析における接続行列による変数の変換と形状関数の変換と、電磁熱流体解析における関数空間の関係を示す。図５に、電磁熱流体解析における双対空間と接続行列の関係説明図を示す。図６に、電磁熱流体解析における木・補木ゲージによる正則化説明図と、電磁熱流体解析におけるA-V法の有限要素方程式を示す。図７に、電磁熱流体解析における積分の説明図を示す。図８、９に、電磁熱流体解析の計算例を示す。
【００２１】
図２において、入力手段１は、構造データやパラメータや初期値などを入力する手段である。構造データ記憶手段２は、解析対象の構造データを記憶保持する手段である。パラメータ記憶手段３は、解析対象のパラメータを記憶保持する手段である。初期値記憶手段４は、解析対象の初期値を記憶保持する手段である。GSMAC-FEM計算手段５は、運動方程式と連続の方程式に高Reynolds数に安定なGSMAC有限要素を適用して流体場を計算する流体解析手段である。CIP-FEM計算手段６は、エネルギー方程式に高Prandtl数に安定なCIP有限要素を適用して温度場の影響を計算する熱流体解析手段である。A-φFEM計算手段７は、電磁力にMaxwellの応力テンソルの発散型を用いた高Hartmann数と高磁気Reynolds数に安定な方程式に辺有限要素を適用して電磁場の影響を計算する電磁流体解析手段である。
【００２２】
中間データ記憶手段８は、１ステップ分時間進行した中間結果を記憶保持する手段である。第１データ変換手段９は、GSMAC-FEM流体解析手段の出力データをCIP-FEM熱流体解析手段とA-φFEM電磁流体解析手段への入力データに変換する手段である。第２データ変換手段10は、CIP-FEM熱流体解析手段の出力データをA-φFEM電磁流体解析手段とGSMAC-FEM流体解析手段への入力データに変換する手段である。第３データ変換手段11は、A-φFEM電磁流体解析手段の出力データをGSMAC-FEM流体解析手段とCIP-FEM熱流体解析手段への入力データに変換する手段である。出力手段12は、中間結果や最終結果を出力する手段である。演算制御手段13は、各手段の演算順序などを制御する手段である。
【００２３】
上記のように構成された本発明の実施例における電磁熱流体解析装置の機能と動作を説明する。最初に、図１を参照しながら、電磁熱流体解析装置の動作手順の概要を説明する。入力手段から、解析対象の構造データやパラメータや初期値などを入力し、構造データ記憶手段とパラメータ記憶手段と初期値記憶手段に保持する。まず(1)電磁力にMaxwellの応力テンソルの発散型を用いた運動方程式と連続の方程式に高Reynolds数と高Hartmann数に安定なGSMAC有限要素法を適用して流体場を計算するGSMAC-FEM流体解析手段で、微小時間の流体の力学的変化を解析して流体の速度を求め、その結果を次の計算用にデータ変換する。次に(2)エネルギー方程式にCIP有限要素を適用した高Prandtl数に安定なCIP-FEM流体解析手段で、流体の熱力学的変化を解析して流体の温度を求め、次の計算用にデータ変換する。さらに(3)電磁力にMaxwellの応力テンソルの発散型を用いて辺有限要素を適用した高磁気Reynolds数に安定なA-φFEM流体解析手段で、流体の電磁気学的変化を解析してベクトルポテンシャルと電位を求めて、次の計算用にデータ変換する。それぞれの変換後のデータは、変換直後に記憶手段に保管される。(1)から(3)までの一連の計算が終了したら、時間を1ステップ進めて、(1)から(3)までの一連の計算を所定の時間になるまで繰り返し、電磁場・温度場・流体場を同時に解く。計算が終了したら結果を出力する。
【００２４】
GSMAC-FEMの時間進行法では、速度の予測子を求めてから、離散化Poisson方程式を反復法で解く。求めた結果から速度・圧力を修正して、新しい時間の物理量を求める方法で、時間進行する。次に、時間進行後の物理量を温度場の離散化式に適用して、収束するまで計算を繰り返して解く方法で、流体場に温度場を弱連成させる。さらに、電磁場の離散化式の速度起電力項に、時間進行後の速度を適用して、収束するまで計算を繰り返して解く方法で、流体場に電磁場を弱連成させる。このようにして、流体場と温度場と電磁場を順次に解きながら時間進行を行う。
【００２５】
レイノルズ数（Reynolds number）は、慣性力と粘性力との比で定義される無次元数である。磁気Reynolds数は、磁気流体力学での層流の安定性を示す無次元パラメータである。Prandtl数は、流体の熱伝導性を示す無次元パラメータであり、流体の動粘度と熱拡散率の比である。Hartmann数は、磁気誘導によって生じる抵抗力が、粘性力に対してどの程度の大きさかを示す無次元数である。高Reynolds数で安定になるようにするには、例えば、特許文献２〜４、特開2002-328954号公報に開示された方法で、GSMAC-FEM流体解析手段を構成する。高Prandtl数で安定になるようにするには、例えば、非特許文献１に開示された方法で、CIP-FEM流体解析手段を構成する。高磁気Reynolds数と高Hartmann数で安定になるようにするには、例えば、非特許文献２に開示された方法で、A-φFEM流体解析手段に、速度起電力を考慮に入れた方程式を用いる。
【００２６】
次に、図２を参照しながら、電磁熱流体解析装置の動作の詳細を説明する。入力手段１を介して、構造データと、パラメータと、初期値を入力し、記憶手段２〜４に格納する。まず(1)GSMAC-FEM流体解析手段５により流体場を計算して速度を求め、その結果をCIP-FEM熱流体解析手段６とA-φFEM電磁流体解析手段７に入力できる形式にデータ変換手段９で変換する。次に(2)CIP-FEM熱流体解析手段６により温度場を計算して温度を求め、計算結果をA-φFEM計算手段７とGSMAC-FEM流体解析手段５に入力できる形式にデータ変換手段１０で変換する。次に(3)A-φFEM電磁流体解析手段７により電磁場を計算して、ベクトルポテンシャルと磁束を求めて、その結果をGSMAC-FEM流体解析手段５とCIP-FEM熱流体解析手段６に入力できる形式にデータ変換手段１１で変換する。データ変換手段９、１０および１１で変換されたデータは、中間データとして中間データ記憶手段８に保存される。中間データは、ある時刻における計算結果であるとともに、その次の時間ステップでの計算にも用いられる。計算結果は、出力手段１２を用いて適宜出力される。構造データは、計算の対象となる場の大きさ、形状、材質などが該当し、パラメータには流体の性質(比熱、密度、Prandtl数、レイノルズ数など)や、解析時のメッシュの数と大きさなどが該当する。また時刻t=0における温度や速度、電磁場などは初期値として与えられる。また、必要に応じて、計算結果や中間データをデータ記憶手段１４に保存するようにしてもよく、またデータ記憶手段１４から出力できるようにしてもよい。
【００２７】
次に、図３を参照しながら、電磁熱流体解析における支配方程式を説明する。流体場は、連続の式と運動方程式で表される。これを解いて、速度vが求められる。速度vは各要素の節点のデータとして求められる。速度vを温度場・電磁場に繰り込むときに、それぞれの場に適した形に変換する。温度場は、熱エネルギー方程式で表される。これを解いて、温度Ｔが求められる。温度Ｔは各要素の節点のデータとして求められる。温度Ｔを電磁場・流体場に繰り込むときに、各場に適した形に変換する。電磁場は、修正Ampereの法則と電荷保存則で表される。これを解いて、ベクトルポテンシャルＡと磁束φが求められる。ベクトルポテンシャルＡと磁束φを変換して、速度vの計算と温度Ｔの計算に利用する。
【００２８】
次に、図４〜図８を参照しながら、計算方法の細部を説明する。流体場にはGSMAC-FEMを適用する。温度場にはCIP-FEMを適用する。渦電流場解析には色々な解析手法がある。有限積分法のFDTD法は電磁場の物理構造をうまく利用した時間領域の解を得る良い手法である。一方、複雑形状物体まわりの電磁場解析には汎用性のある有限要素法が適している。ベクトル形状関数を用いた辺有限要素法はこれに適している。この手法も電磁場の物理構造をうまく取り込んでいる。基本的にはすべての領域をA-φ法で定式化する。空気領域に境界要素法を適用する場合も多いが、ここでは、全領域を有限要素法で解く。これにより、単純化と高速化が期待できる。導体が運動する場合、誘導方程式の移流項を考慮する必要がある。移流項は場の不安定性と密接に関係している。この原因は、磁気Reynolds数が高くなると、磁力線が流線に凍結されることによる。以下、主に高磁気Reynolds数に安定な渦電流の時間領域での解析手法について説明する。
【００２９】
流体場・温度場・電磁場の方程式をGalerkin有限要素法で定式化する。流体場の計算方法は、線形のスカラー形状関数を基底関数としたGSMAC法である。温度場の計算方法は、関数と微係数が連続となるスカラー形状関数を基底関数としたCIP法である。電磁場の計算方法は、微分形式に基づいたベクトル形状関数を基底関数とした辺有限要素法（A-φ法）である。流体場・温度場・電磁場よりなる場の支配方程式について、以下に説明する。
【００３０】
場は流体場・温度場・電磁場よりなる。流体場の支配方程式は、
ρ0(dｖ/dt)＝−gradｐ＋div(2ηｄ)＋Ｊ×Ｂ−ρ0(βΔＴ)ｇ
divｖ＝０
である。温度場の支配方程式は、
ρｃv(dＴ/dt)＝div(ｋgradＴ)＋ｓ[Ｔ]
ｓ[Ｔ]＝2ηｄ:ｄ＋Ｊ2/σ
である。電磁場の支配方程式は、
rotＥ＝−(∂Ｂ/∂t)
rotＨ＝Ｊ
divＢ＝０
divＪ＝０
である。
【００３１】
流体場は、Boussinesq近似された運動方程式
ρ0(dｖ/dt)＝−gradｐ＋div(２ηｄ)＋Ｆem-ρ0(βΔＴ)ｇ
と非圧縮連続の方程式
divｖ＝０
により表現される。ここで、ρ0は基準密度（スカラー）、ｖは速度ベクトル、ｐは圧力（スカラー）、ηは粘性係数（スカラー）、ｄは変形速度テンソル、Ｆemは電磁力（ベクトル）で、Maxwellの応力テンソルＴemの発散として計算される。即ち、
Ｆem＝divＴem＝div(ＢＨ−(1/2)(Ｂ・Ｈ)Ｉ)
である。ここで、Ｂは磁束密度（ベクトル）、Ｈは磁場（ベクトル）、Ｉは恒等テンソルである。ＢＨは、Ｂi・Ｈjを（ｉ,ｊ）要素とするテンソルである。テンソルの発散は、テンソル場の偏微分係数を縮約して得られるベクトルである。ｄ:ｄは、tr(ｄｄt)＝Σijｄij2である。Boussinesq近似において密度は、
ρ＝ρ0(１−βΔＴ)
と表せる。ここで、βは体膨張率（スカラー）、ΔＴ＝Ｔ−Ｔ0は基準温度Ｔ0からの温度差（スカラー）、ｇは重力加速度（ベクトル）である。
【００３２】
温度場は、エネルギー方程式
ρｃv(dＴ/dt)＝div(ｋgradＴ)+ｓ[Ｔ]
により表現される。ここで、ｃvは定容比熱（スカラー）、ｋは熱伝導率（スカラー）、ｓ[Ｔ]は発熱項である。発熱項は、粘性発熱項２ηｄ:ｄとJoule発熱項Ｊ2/σの和である。ここで、Ｊは電流密度（ベクトル）、σは電気伝導度（スカラー）である。
【００３３】
電磁場は、準定常近似のMaxwellの方程式
rotＥ＝−(∂Ｂ/∂t)
rotＨ＝Ｊ
divＢ＝０
により表現される。電流の発散が零となる両立性の条件は
divＪ＝０
である。A-φ法で定式化すると
Ｂ＝rotＡ
Ｅ＝−(∂Ａ/∂t)−gradφ＝−(∂/∂t)(Ａ＋gradＶ)
となる。ここで、φは電気スカラーポテンシャル、Ａは磁気ベクトルポテンシャルである。Ｖは対称化ポテンシャルで、
φ＝∂Ｖ/∂t
である。このとき、
divＢ＝０
rotＥ＝−(∂Ｂ/∂t)
は自動的に満たされる。
【００３４】
よって、解くべき方程式は
rotＨ＝Ｊ
のみである。これより、誘導方程式
rot[(1/σ)(1/μ)rotＢ]＝−(∂Ｂ/∂t)−rot(Ｂ×ｖ)
が導ける。μは透磁率（スカラー）である。これをA-V法で表すと
rot(νrotＡ)＝−σ((∂/∂t)(Ａ＋gradＶ)＋(rotＡ)×ｖ)
となる。また、両立性の条件は
div(σ[(∂/∂t)(Ａ＋gradＶ)＋(rotＡ)×ｖ])＝０
となる。ここで、Ｖ＝０のゲージを選択すればＡ法となる。静止導体の場合、速度はｖ＝０である。境界条件は
ｎ×Ａ＝０
Ｖ＝const.
である。これらの式を辺有限要素法で解く。
【００３５】
次に、接続行列について説明する。接続行列は特異行列で、勾配行列・回転行列・発散行列の３種類のものがある。
{∂ｅ}＝[Ｇ]{ｎ}
{∂ｆ}＝[Ｒ]{ｅ}
{∂ｖ}＝[Ｄ]{ｆ}
これらは、辺の境界が接点、面の境界が辺、体の境界が面であることを表している。ここで、接続行列の成分表示は
[Ｇ]＝Ｇ＝Ｇen＝Ｇen
[Ｒ]＝Ｒ＝Ｒfe＝Ｒfe
[Ｄ]＝Ｄ＝Ｄvf＝Ｄvf
であり、これらの転置行列は
[Ｇ]t＝Ｇt＝Ｇne＝Ｇne
[Ｒ]t＝Ｒef
[Ｄ]t＝Ｄfv
である。このように、行列の表示方法には色々ある。
【００３６】
[Ｇ]、[Ｒ]、[Ｄ]は成分が−１、０、＋１からなる接続行列である。このとき、
∂∂＝０  →  [Ｒ][Ｇ]＝０、[Ｄ][Ｒ]＝０
が成立する。これは、外微分演算
dd＝０  →  (rot)(grad)＝０、(div)(rot)＝０
に対応する。このとき、次の対応関係
Ｂ＝rotＡ  ←→  {Ｂ}＝[Ｒ]{Ａ}  ←→  Ｂf＝ＲfeＡe
Ｅ＝−gradφ−(∂Ａ/∂t)  ←→  {Ｅ}＝−[Ｇ]{φ}−(∂{Ａ}/∂t)
が成立する。
【００３７】
[Ｒ][Ｇ]＝０の証明
∂{∂ｆ}＝∂[Ｒ]{ｅ}＝[Ｒ]{∂ｅ}＝[Ｒ][Ｇ]{ｎ}＝０
∴[Ｒ][Ｇ]＝０
【００３８】
{Ｂ}＝[Ｒ]{Ａ}の証明
面ｆを通過する磁束を計算する。磁束密度を、面基底要素を用いて、
Ｂ＝ΣfＭfＢf
と表す。このとき、直交性
∫f'Ｍf・ｎdＳ＝δff'
が成立する。よって、Ｂの法線成分
∫fＢ・ｎdＳ＝Ｂf
は境界面上で連続である。Ｂの接線成分は一般に不連続である。つぎに、Ｂの法線成分をＡの接線成分で表現する。
Ｂf＝∫fＢ・ｎdＳ＝∫f(rotＡ)・ｎdＳ＝∫∂fＡ・dｓ
＝∫Ａ・dｓ（積分範囲はΣeＲfeｅ）
＝ΣeＲfe∫eＡ・dｓ＝ΣeＲfeＡe
ゆえに、
Ｂf＝ΣeＲfeＡe  →  {Ｂ}＝[Ｒ]{Ａ}  ←→  Ｂ＝rotＡ
を得る。
【００３９】
微分演算子と接続行列の関係をまとめて示す。微分演算子と接続行列は次の通りである。
{∂ｅ}＝Ｇ{ｎ}    Ｅ＝gradφ    {Ｅe}＝Ｇ{φn}
{∂ｆ}＝Ｒ{ｅ}    Ｂ＝rotＡ     {Ｂf}＝Ｒ{Ａe}
{∂ｖ}＝Ｄ{ｆ}    ρ＝divＢ     {ρv}＝Ｄ{Ｂf}
【００４０】
Stokesの定理を示す。
∫Ωdω＝∫∂Ωω
dと∂は双対である。
〈Ω, dω〉＝〈∂Ω, ω〉
Ｅe＝∫eＥ・dｒ＝∫∂eφ＝Ｇen∫nφ＝Ｇenφn
∴Ｅe＝Ｇenφn
Ｂf＝∫fＢ・ｎdＳ＝∫∂fＡ・dｒ＝Ｒfe∫eＡ・dｒ＝ＲfeＡe
∴Ｂf＝ＲfeＡe
ρv＝∫vρdＶ＝∫∂vＢ・ｎdＳ＝Ｄvf∫fＢ・ｎdＳ＝ＤvfＢf
∴ρv＝ＤvfＢf
【００４１】
Maxwellの方程式と接続行列の関係をまとめて示す。基本場と物質場のMaxwellの方程式と接続行列は次の通りである。
基本場（基本空間）              物質場（双対空間）
divＢ＝０, [Ｄ]{Ｂ}＝０         divＤ＝ρe, [Ｇ]t{Ｄ}＝{ρe}
rotＥ＝−(∂Ｂ/∂t),            rotＨ＝Ｊ＋(∂Ｄ/∂t)
[Ｒ]{Ｅ}＝−(∂{Ｂ}/∂t)        [Ｒ]t{Ｈ}＝{Ｊ}＋(∂{Ｄ}/∂t)
Ｂ＝rotＡ, {Ｂ}＝[Ｒ]{Ａ}       Ｈ＝−gradφm, {Ｈ}＝−[Ｄ]t{φm}
Ｅ＝−(∂Ａ/∂t)−gradφ
{Ｅ}＝−(∂{Ａ}/∂t)−[Ｇ]{φ}
【００４２】
次に、基底間の関係式について説明する。基底の間につぎの関係式が成立する。
{rotＮ}＝[Ｒ]t{Ｍ}  →  ＶrotN⊂ＶM
{gradＮ}＝[Ｇ]t{Ｎ}  →  ＶgradN⊂ＶM
ここで、ＶrotNは、辺基底の回転rotＮeが張る空間である。ＶMは、面基底Ｍfが張る空間である。ＶgradNは、節点基底の勾配gradＮnが張る空間である。ＶNは、辺基底Ｎeが張る空間である。
【００４３】
{rotＮ}＝[Ｒ]t{Ｍ}を証明する。まず
Ｂ＝rotＡ＝ΣeＡerotＮe
である。一方、
Ｂ＝ΣfＭfＢf＝ΣfＭf∫fＢ・ｎdＳ＝ΣfＭf∫∂fＡ・dｓ
＝ΣfＭf(ΣeＲfeＡe)＝ΣeＡe(ΣfＲfeＭf)
である。両式を比較して
rotＮe＝ΣfＲfeＭf  ←→  {rotＮ}＝[Ｒ]t{Ｍ}
を得る。これは、ＶrotN⊂ＶMを意味する（証明終わり）。
【００４４】
次に、接続行列のランクについて説明する。有限要素の総節点数をNd、総辺数をEd、総面数をFc、総要素数をVmとする。独立な辺の数は木の枝の数に等しく、Nd−１である。よって、
{∂ｅ}＝[Ｇ]{ｎ}
rank[Ｇ]＝Nd−１
である。独立な面の数は補木の弦の数に等しく、（Fc−１）である。よって、
{∂ｆ}＝[Ｒ]{ｅ}
rank[Ｒ]＝Fc−１
である。そして、辺数の総数Edは、
（木の枝の数）＋（補木の弦の数）＝（Nd−１）＋（Fc−１）＝Ed
となる。ここで、単連結領域の場合、FcとEdとNdの間にEulerの関係式
Fc−Ed＋Nd＝２
が成立する。
【００４５】
図４（ａ）に、接続行列による変数の変換と形状関数の変換をまとめて示す。図４（ｂ）に、関数空間の関係をまとめて示す。ここで、節点、辺、面、体の基底形状関数は、
Ｎn＝ｗn, Ｎe＝ｗe, Ｍf＝ｗf, Ｎv＝ｗv
であり、基底空間は
ＶgradN＝gradＷ0⊂ＶN＝Ｗ1, ＶrotN＝rotＷ1⊂ＶM＝Ｗ2
である。
【００４６】
図５に、双対空間（基本場と物質場）と接続行列の関係式をまとめて示す。双対空間の接続行列は基本空間の転置行列である。ここで、物質場（双対空間）におけるgrad rot divの接続行列は、
Ｇ~＝Ｄt
Ｒ~＝Ｒt
Ｄ~＝Ｇt
となる。
【００４７】
次に、静磁場と辺有限要素法について説明する。磁場には静磁場と動電磁場の問題がある。静磁場は動電磁場の特別な場合であるが、動電磁場を理解するのに静磁場解析を理解しておくのが先決である。ここでは、静磁場の解析を、辺有限要素法を用いて行う。すなわち、ベクトル形状関数を用いた辺有限要素法で定式化する。磁場の支配方程式
rot(νrotＡ)＝Ｊ  (ν＝1/μ)
に重み付き残さ法を適用する。ベクトル重み関数ｗを内積して領域積分すると
∫Ωｗ・rot(νrotＡ)dΩ＝∫Ωｗ・ＪdΩ
となる。ここで、微分公式
div[(rotＡ)×ｗ]＝ｗ・[rot(rotＡ)]−(rotｗ)・(rotＡ)
を用いる。部分積分して発散項が零になることを用いると、弱形式
∫Ω(rotｗ)・(νrotＡ) dΩ＝∫Ωｗ・ＪdΩ
が導ける。辺要素基底（ベクトル形状関数）Ｎeを適用し、Ａを
Ａ＝ΣeＮeＡe
と補間する。重み関数ｗに同じ辺要素基底を適用し、
ｗ＝Ｎe
とすれば、Galerkin有限要素法となる。辺要素基底は直交性
∫ee'Ｎe・dｓ＝δee'
を満足し、要素境界上でＡの接線成分が連続である。すなわち、
∫eＡ・dｓ＝Ａe
となる。Ａの法線成分は一般に不連続である。
【００４８】
有限要素方程式について説明する。弱形式で定式化された方程式にＡとｗの辺要素基底で表された式を代入すると有限要素法の離散式
Σe'Ａe'∫Ω(rotＮe)・(νrotＮe')dΩ＝∫ΩＮe・ＪdΩ
が導ける。ここで、和はすべての辺ｅ＝１,２,…,Edの和である。Edは全辺数である。ここで、
{rotＮ}＝[Ｒ]t{Ｍ}
を用いる。上式より、行列方程式
[Ｋ]{Ａ}＝{Ｊ}
を得る。ただし、
Ｋee'＝∫Ω(rotＮe)・(νrotＮe')dΩ
Ｊe＝∫ΩＮe・ＪdΩ
である。
【００４９】
次に、有限要素行列の性質について説明する。有限要素行列は特異行列で正則行列でない。すなわち、
[Ｋ]＝[Ｒ]t[ν][Ｒ]：ｅ×ｅ行列
[Ｒ]t：ｅ×ｆ行列  [Ｒ]：ｆ×ｅ行列
[ν]：対称正則ｆ×ｆ行列
νff'＝∫ΩνＭf・Ｍf'dΩ
である。行列[Ｋ]のdimとrankは異なる。すなわち、
dim[Ｋ]＝Ed
rank[Ｋ]＝rank[Ｒ]＝Ed−Nd＋１
である。よって、行列[Ｋ]は特異行列である。
【００５０】
次に、木・補木ゲージによる正則化について説明する。木の自由度を消去して、[Ｋ]を正則化して[Ｋ’]とする(図６（ａ）参照)。補木の弦の数は独立なループの個数に等しい。木の自由度を消去して正則化すると、最小固有値が減少し、条件数が悪化する場合が極めて多い。よって、特異行列をそのままICCG法で解いた方が良い。
【００５１】
誘導方程式と辺有限要素法について説明する。動電磁場の支配方程式は誘導方程式である。この誘導方程式をA-φ法で記述する。この場合、電流は励振電流（ソース電流）Ｊsと誘導電流Ｊiの和として
Ｊ＝Ｊs＋Ｊi
と表せる。そして、誘導電流は
Ｊi＝σ(Ｅ＋ｖ×Ｂ)＝−σ((∂Ａ/∂t)＋(rotＡ)×ｖ＋gradφ)
となる。
【００５２】
A-φ法の基礎方程式について説明する。A-φ法の基礎方程式と両立性の条件はつぎのとおりである。
rot(νrotＡ)＝Ｊs＋Ｊi
＝Ｊs−σ((∂Ａ/∂t)＋(rotＡ)×ｖ＋gradφ)
div(σ[gradφ＋(∂Ａ/∂t)＋(rotＡ)×ｖ])＝０
この方程式を基礎に色々な変形が可能である。
【００５３】
A-V法の基礎方程式について説明する。A-φ法の基礎方程式において、添え字を省略して
Ｊs＝Ｊ  φ＝∂Ｖ/∂t
と置くと、A-V法の基礎方程式
rot(νrotＡ)＋σ((∂(Ａ＋gradＶ)/∂t)＋(rotＡ)×ｖ)＝Ｊ
div(σ[(∂(Ａ＋gradＶ)/∂t)＋(rotＡ)×ｖ])＝０
が得られる。ここで、Ａを辺要素で補間し、Ｖを節点要素で補間すると
Ａ＝ΣeＮeＡe
Ｖ＝ΣnＮnＶn
となる。速度ｖ＝０の場合を考える。これより弱形式
∫Ω(rotＮe)・ν(rotＡ)dΩ＋ｊω∫ΩσＮe・(Ａ＋gradＶ)dΩ
＝∫ΩＮe・ＪdΩ
ｊω∫ΩσgradＮn・(Ａ＋gradＶ)dΩ＝０
が導ける。よって、A-V法の有限要素方程式は図６（ｂ）に示すようになる。さらに、[σ]、[ν]は正定値行列で、その成分は、
σee'＝∫ΩσＮe・Ｎe'dΩ
νff'＝∫ΩνＭf・Ｍf'dΩ
である。
【００５４】
Ａ法の基礎方程式について説明する。A-V法の基礎方程式で、Ｖ＝０のゲージを採用すればＡ法の行列方程式
[Ｋ]{Ａ}＝{Ｊ}
が得られる。ここで、
[Ｋ]＝[Ｒ]t[ν][Ｒ]＋ｊω[σ]
である。この式は単純化されているが、収束性はＡ法よりもA-V法の方が一般に勝っている。
【００５５】
時間領域積分について説明する。Ａ法の基礎方程式として
rot(νrotＡ)＋σ((∂Ａ/∂t)＋(rotＡ)×ｖ)＝０
div(σ[(∂Ａ/∂t)＋(rotＡ)×ｖ])＝０
を考える。弱形式で定式化すると
∫Ω(rotＮe)・ν(rotＡ)dΩ＋∫ΩσＮe・((∂Ａ/∂t)＋(rotＡ)×ｖ)dΩ
＝０
となる。よって、離散化式
【数１】
2012074016.tif 000003
を得る。ここで、
[Ｋ]＝[Ｒ]t[ν][Ｒ]−[Ａd]
[Ａd]＝∫ΩσＮe・(ｖ×rotＮe')dΩ
[σ]＝∫ΩσＮe・Ｎe'dΩ
である。時間積分には一般化台形公式
{θ[Ｋn+1]＋[σ]/Δt}{Ａn+1}＝−{(１−θ)[Ｋn]−[σ]/Δt}{Ａn}
を用いる。ここで、θは混合パラメータで、この方程式は陰解法である。両立性の条件は
−∫ΩσgradＮn・((∂Ａ/∂t)＋(rotＡ)×ｖ)dΩ
＝∫ΩσgradＮn・ＪidΩ＝０
となる。この式は、divＪi＝０と同等である。よって、大域的には境界において電流が流出しない条件
ｎ・Ｊi＝０
と同等となる。よって、この条件を満足しないと解は不定となる。
【００５６】
流体場には温度差が浮力として作用し、Lorentz力として電磁力が働く。温度場には熱源項として、粘性散逸の発熱項とJoule発熱項がある。電磁場には導体の運動による誘導起電力がある。すなわち、流体場と温度場と電磁場が互いに連成する場合を考える。そして、これら３つの場に、次に示す３種類の有限要素法を適用する。流体場には、スカラー形状関数を用い、GSMAC-FEMを適用する。温度場には、スカラー形状関数を用い、CIP−FEMを適用する。電磁場には、ベクトル形状関数を用い、辺有限要素法（A-φ法）を適用する。これにより、高Reynolds数、高Prandtl数、高磁気Reynolds数に安定な解が得られる。有限要素法に用いられる基底関数（形状関数）には、節点基底関数、辺基底関数、面基底関数、体基底関数の４種類がある。節点基底関数と体基底関数はスカラー関数であるが、辺基底関数と面基底関数はベクトル関数である。これらの基底関数は、図７に示すように、微分形式と密接の関係している。そして、積分値が未知変数となっている。特に電磁場においては、Gaussの定理やStokesの定理が重要で、これらの定理を辺有限要素法はうまく表現している。
【００５７】
次に、図８を参照しながら、電磁熱流体が空気中におかれた場合の、電磁熱流体解析の計算例を説明する。電磁熱流体としては、金属微粒子を含有する溶融ガラス(MHD)とする。また空気中には電磁場が存在するとする。図８（ａ）に解析モデルを示す。大立方体の中心の小立方体が電磁熱流体の領域である。ｘ方向に電位差が生じるように、小立方体の端面の電位を設定している。−ｚ方向に重力がかかり、ｚ方向に磁場がある。図８（ｂ）に主な解析パラメータを示す。図８（ｃ）に解析メッシュを示す。解析メッシュは６面体要素に分割されている。
【００５８】
図９（ａ）、（ｂ）、（ｃ）に計算開始から約2000秒後の計算結果を示す。図９（ａ）は流体場の、図９（ｂ）は温度場の解析結果である。また、電磁場の解析の結果得られた例として磁束密度の分布を図９（ｃ）に示す。白黒でかつ簡易的に表示してあるので明瞭ではないが、小立方体上部でｘ断面において流体の速度が大きくなっているところがあることが見て取れる。また温度は、小立方体最上部よりも内部のほうが高温になっていることがわかる。磁束密度は小立方体の中心部にはほとんど無く周辺部で高くなっていることがわかる。
【００５９】
上記のように、本発明の実施例では、電磁熱流体解析装置を、GSMAC-FEM流体解析手段で微小時間の流体の力学的変化を解析してデータ変換し、CIP-FEM流体解析手段で流体の熱力学的変化を解析してデータ変換し、A-φFEM流体解析手段で流体の電磁気学的変化を解析してデータ変換し、これを繰り返して電磁熱流体の動きを解析する構成としたので、電磁場・温度場・流体場が同時に解けて、高速高精度で汎用性が高い解析システムが実現できる。
【産業上の利用可能性】
【００６０】
本発明の電磁熱流体解析装置は、流体場・温度場・電磁場が密接に関連している実用的な多くの問題を解析する装置として最適である。例えば、シリコン単結晶装置の開発や、超電導船の推進装置の開発や、プラズマ・ジェット装置の開発や、核融合の研究などに利用できる。
【符号の説明】
【００６１】
１  入力手段
２  構造データ記憶手段
３  パラメータ記憶手段
４  初期値記憶手段
５  GSMAC-FEM計算手段
６  CIP-FEM計算手段
７  A-φFEM計算手段
８  中間データ記憶手段
９  第１データ変換手段
10  第２データ変換手段
11  第３データ変換手段
12  出力手段
13  演算制御手段
【図１】
2012074016.tif 000004
【図２】
2012074016.tif 000005
【図３】
2012074016.tif 000006
【図４】
2012074016.tif 000007
【図５】
2012074016.tif 000008
【図６】
2012074016.tif 000009
【図７】
2012074016.tif 000010
【図８】
2012074016.tif 000011
【図９】
2012074016.tif 000012
【図１０】
2012074016.tif 000013