構造データとパラメータと初期値を入力する入力手段と、入力した構造データとパラメータと初期値を格納する記憶装置と、電磁力にMaxwellの応力テンソルの発散型を用いた運動方程式と連続の方程式に高Reynolds数と高Hartmann数に安定なGSMAC有限要素法を適用して流体場を計算するGSMACFEM流体解析手段と、エネルギー方程式に高Prandtl数に安定なCIP有限要素法を適用して温度場を計算するCIPFEM熱流体解析手段と、高磁気Reynolds数に安定な方程式に辺有限要素法を適用して電磁場を計算するAφFEM電磁流体解析手段と、GSMACFEM流体解析手段の出力データをCIPFEM熱流体解析手段とAφFEM電磁流体解析手段への入力データに変換する第1データ変換手段と、CIPFEM熱流体解析手段の出力データをAφFEM電磁流体解析手段とGSMACFEM流体解析手段への入力データに変換する第2データ変換手段と、AφFEM電磁流体解析手段の出力データをGSMACFEM流体解析手段とCIPFEM熱流体解析手段への入力データに変換する第3データ変換手段と、ある時刻においてGSMACFEM流体解析手段とCIPFEM熱流体解析手段とAφFEM電磁流体解析手段を連続して順次動作させてそれぞれの解を求めるように制御する計算制御手段とを具備することを特徴とする電磁熱流体解析装置。流体場温度場電磁場が密接に関連している電磁熱流体の運動をコンピュータにより数値解析する電磁熱流体解析方法であって、入力手段を介して入力された構造データと初期値とパラメータとを記憶手段に格納し、電磁力にMaxwellの応力テンソルの発散型を用いた運動方程式と連続の方程式に高Reynolds数と高Hartmann数に安定なGSMAC有限要素を適用して流体を解析するGSMACFEM流体解析手段で流体場を計算し、GSMACFEM流体解析手段の出力データをCIPFEM熱流体解析手段とAφFEM電磁流体解析手段への入力データに第1データ変換手段で変換し、エネルギー方程式に高Prandtl数に安定なCIP有限要素を適用して流体の温度場を解析するCIPFEM熱流体解析手段で温度場の影響を計算し、CIPFEM熱流体解析手段の出力データをAφFEM電磁流体解析手段とGSMACFEM流体解析手段への入力データに第2データ変換手段で変換し、高磁気Reynolds数に安定な方程式に辺有限要素を適用して流体の電磁場を解析するAφFEM電磁流体解析手段で電磁場の影響を計算し、AφFEM電磁流体解析手段の出力データをGSMACFEM流体解析手段とCIPFEM熱流体解析手段への入力データに第3データ変換手段で変換し、ある時刻においてGSMACFEM流体解析手段とCIPFEM熱流体解析手段とAφFEM電磁流体解析手段を連続して順次実施してそれぞれの解を求めることにより電磁熱流体の運動を解析することを特徴とする電磁熱流体解析方法。請求項2に記載の電磁熱流体解析方法をコンピュータに実行させるための処理手順を記述したコンピュータプログラム。請求項3記載のコンピュータプログラムを記録したコンピュータ読取可能な記録媒体。しかし、上記従来の方法では、温度場と流体場を同時に解くことはできるし、電磁場と流体場も同時に解くことができるが、電磁場と温度場と流体場を同時に解くことはできない。実用的な問題の多くは、流体場温度場電磁場が密接に関連している。そこで、それぞれの数値解析システムを使って、工夫をこらして解析している。しかし、この3つの場を統一的に解ける単一の数値解析システムはない。本発明の目的は、上記従来の問題を解決して、電磁場温度場流体場が同時に解けて、高速高精度で汎用性が高い解析システムを実現することである。上記の課題を解決するために、本発明では、電磁熱流体解析装置を、構造データとパラメータと初期値を入力する入力手段と、入力した構造データとパラメータと初期値を格納する記憶装置と、電磁力にMaxwellの応力テンソルの発散型を用いた運動方程式と連続の方程式に高Reynolds数と高Hartmann数に安定なGSMAC有限要素を適用して流体場を計算するGSMACFEM流体解析手段と、エネルギー方程式に高Prandtl数に安定なCIP有限要素を適用して温度場の影響を計算するCIPFEM熱流体解析手段と、高磁気Reynolds数に安定な方程式に辺有限要素を適用して電磁場の影響を計算するAφFEM電磁流体解析手段と、GSMACFEM流体解析手段の出力データをCIPFEM熱流体解析手段とAφFEM電磁流体解析手段への入力データに変換する第1データ変換手段と、CIPFEM熱流体解析手段の出力データをAφFEM電磁流体解析手段とGSMACFEM流体解析手段への入力データに変換する第2データ変換手段と、AφFEM電磁流体解析手段の出力データをGSMACFEM流体解析手段とCIPFEM熱流体解析手段への入力データに変換する第3データ変換手段と、ある時刻においてGSMACFEM流体解析手段とCIPFEM熱流体解析手段とAφFEM電磁流体解析手段を連続して順次動作させてそれぞれの解を求めるように制御する計算制御手段とを具備する構成とした。上記のように構成したことにより、電磁場温度場流体場が同時に解けて、高速高精度で汎用性が高い解析システムを実現でき、3次元でも非定常でも非構造格子でも、誘導磁場や熱対流などを考慮した電磁熱流体の流れが解析できる。以下、本発明を実施するための最良の形態について、図1~図9を参照しながら詳細に説明する。本発明の実施例は、GSMACFEM流体解析手段で微小時間の流体の力学的変化を解析してデータ変換し、CIPFEM熱流体解析手段で流体の熱力学的変化を解析してデータ変換し、AφFEM電磁流体解析手段で流体の電磁気学的変化を解析してデータ変換し、これを順次実行して電磁熱流体の動きを解析する電磁熱流体解析装置である。図1に、電磁熱流体解析装置の動作を示す流れ図を示す。図2に、電磁熱流体解析装置の概念図を示す。主なデータ経路のみを示してあり、当然必要な自明なデータ経路は省略してある。図3に、支配方程式を示す。図4に、電磁熱流体解析における接続行列による変数の変換と形状関数の変換と、電磁熱流体解析における関数空間の関係を示す。図5に、電磁熱流体解析における双対空間と接続行列の関係説明図を示す。図6に、電磁熱流体解析における木補木ゲージによる正則化説明図と、電磁熱流体解析におけるAV法の有限要素方程式を示す。図7に、電磁熱流体解析における積分の説明図を示す。図8、9に、電磁熱流体解析の計算例を示す。図2において、入力手段1は、構造データやパラメータや初期値などを入力する手段である。構造データ記憶手段2は、解析対象の構造データを記憶保持する手段である。パラメータ記憶手段3は、解析対象のパラメータを記憶保持する手段である。初期値記憶手段4は、解析対象の初期値を記憶保持する手段である。GSMACFEM計算手段5は、運動方程式と連続の方程式に高Reynolds数に安定なGSMAC有限要素を適用して流体場を計算する流体解析手段である。CIPFEM計算手段6は、エネルギー方程式に高Prandtl数に安定なCIP有限要素を適用して温度場の影響を計算する熱流体解析手段である。AφFEM計算手段7は、電磁力にMaxwellの応力テンソルの発散型を用いた高Hartmann数と高磁気Reynolds数に安定な方程式に辺有限要素を適用して電磁場の影響を計算する電磁流体解析手段である。中間データ記憶手段8は、1ステップ分時間進行した中間結果を記憶保持する手段である。第1データ変換手段9は、GSMACFEM流体解析手段の出力データをCIPFEM熱流体解析手段とAφFEM電磁流体解析手段への入力データに変換する手段である。第2データ変換手段10は、CIPFEM熱流体解析手段の出力データをAφFEM電磁流体解析手段とGSMACFEM流体解析手段への入力データに変換する手段である。第3データ変換手段11は、AφFEM電磁流体解析手段の出力データをGSMACFEM流体解析手段とCIPFEM熱流体解析手段への入力データに変換する手段である。出力手段12は、中間結果や最終結果を出力する手段である。演算制御手段13は、各手段の演算順序などを制御する手段である。上記のように構成された本発明の実施例における電磁熱流体解析装置の機能と動作を説明する。最初に、図1を参照しながら、電磁熱流体解析装置の動作手順の概要を説明する。入力手段から、解析対象の構造データやパラメータや初期値などを入力し、構造データ記憶手段とパラメータ記憶手段と初期値記憶手段に保持する。まず電磁力にMaxwellの応力テンソルの発散型を用いた運動方程式と連続の方程式に高Reynolds数と高Hartmann数に安定なGSMAC有限要素法を適用して流体場を計算するGSMACFEM流体解析手段で、微小時間の流体の力学的変化を解析して流体の速度を求め、その結果を次の計算用にデータ変換する。次にエネルギー方程式にCIP有限要素を適用した高Prandtl数に安定なCIPFEM流体解析手段で、流体の熱力学的変化を解析して流体の温度を求め、次の計算用にデータ変換する。さらに電磁力にMaxwellの応力テンソルの発散型を用いて辺有限要素を適用した高磁気Reynolds数に安定なAφFEM流体解析手段で、流体の電磁気学的変化を解析してベクトルポテンシャルと電位を求めて、次の計算用にデータ変換する。それぞれの変換後のデータは、変換直後に記憶手段に保管される。からまでの一連の計算が終了したら、時間を1ステップ進めて、からまでの一連の計算を所定の時間になるまで繰り返し、電磁場温度場流体場を同時に解く。計算が終了したら結果を出力する。GSMACFEMの時間進行法では、速度の予測子を求めてから、離散化Poisson方程式を反復法で解く。求めた結果から速度圧力を修正して、新しい時間の物理量を求める方法で、時間進行する。次に、時間進行後の物理量を温度場の離散化式に適用して、収束するまで計算を繰り返して解く方法で、流体場に温度場を弱連成させる。さらに、電磁場の離散化式の速度起電力項に、時間進行後の速度を適用して、収束するまで計算を繰り返して解く方法で、流体場に電磁場を弱連成させる。このようにして、流体場と温度場と電磁場を順次に解きながら時間進行を行う。レイノルズ数Reynoldsnumberは、慣性力と粘性力との比で定義される無次元数である。磁気Reynolds数は、磁気流体力学での層流の安定性を示す無次元パラメータである。Prandtl数は、流体の熱伝導性を示す無次元パラメータであり、流体の動粘度と熱拡散率の比である。Hartmann数は、磁気誘導によって生じる抵抗力が、粘性力に対してどの程度の大きさかを示す無次元数である。高Reynolds数で安定になるようにするには、例えば、特許文献2~4、特開2002328954号公報に開示された方法で、GSMACFEM流体解析手段を構成する。高Prandtl数で安定になるようにするには、例えば、非特許文献1に開示された方法で、CIPFEM流体解析手段を構成する。高磁気Reynolds数と高Hartmann数で安定になるようにするには、例えば、非特許文献2に開示された方法で、AφFEM流体解析手段に、速度起電力を考慮に入れた方程式を用いる。次に、図2を参照しながら、電磁熱流体解析装置の動作の詳細を説明する。入力手段1を介して、構造データと、パラメータと、初期値を入力し、記憶手段2~4に格納する。まずGSMACFEM流体解析手段5により流体場を計算して速度を求め、その結果をCIPFEM熱流体解析手段6とAφFEM電磁流体解析手段7に入力できる形式にデータ変換手段9で変換する。次にCIPFEM熱流体解析手段6により温度場を計算して温度を求め、計算結果をAφFEM計算手段7とGSMACFEM流体解析手段5に入力できる形式にデータ変換手段10で変換する。次にAφFEM電磁流体解析手段7により電磁場を計算して、ベクトルポテンシャルと磁束を求めて、その結果をGSMACFEM流体解析手段5とCIPFEM熱流体解析手段6に入力できる形式にデータ変換手段11で変換する。データ変換手段9、10および11で変換されたデータは、中間データとして中間データ記憶手段8に保存される。中間データは、ある時刻における計算結果であるとともに、その次の時間ステップでの計算にも用いられる。計算結果は、出力手段12を用いて適宜出力される。構造データは、計算の対象となる場の大きさ、形状、材質などが該当し、パラメータには流体の性質や、解析時のメッシュの数と大きさなどが該当する。また時刻t0における温度や速度、電磁場などは初期値として与えられる。また、必要に応じて、計算結果や中間データをデータ記憶手段14に保存するようにしてもよく、またデータ記憶手段14から出力できるようにしてもよい。次に、図3を参照しながら、電磁熱流体解析における支配方程式を説明する。流体場は、連続の式と運動方程式で表される。これを解いて、速度vが求められる。速度vは各要素の節点のデータとして求められる。速度vを温度場電磁場に繰り込むときに、それぞれの場に適した形に変換する。温度場は、熱エネルギー方程式で表される。これを解いて、温度Tが求められる。温度Tは各要素の節点のデータとして求められる。温度Tを電磁場流体場に繰り込むときに、各場に適した形に変換する。電磁場は、修正Ampereの法則と電荷保存則で表される。これを解いて、ベクトルポテンシャルAと磁束φが求められる。ベクトルポテンシャルAと磁束φを変換して、速度vの計算と温度Tの計算に利用する。次に、図4~図8を参照しながら、計算方法の細部を説明する。流体場にはGSMACFEMを適用する。温度場にはCIPFEMを適用する。渦電流場解析には色々な解析手法がある。有限積分法のFDTD法は電磁場の物理構造をうまく利用した時間領域の解を得る良い手法である。一方、複雑形状物体まわりの電磁場解析には汎用性のある有限要素法が適している。ベクトル形状関数を用いた辺有限要素法はこれに適している。この手法も電磁場の物理構造をうまく取り込んでいる。基本的にはすべての領域をAφ法で定式化する。空気領域に境界要素法を適用する場合も多いが、ここでは、全領域を有限要素法で解く。これにより、単純化と高速化が期待できる。導体が運動する場合、誘導方程式の移流項を考慮する必要がある。移流項は場の不安定性と密接に関係している。この原因は、磁気Reynolds数が高くなると、磁力線が流線に凍結されることによる。以下、主に高磁気Reynolds数に安定な渦電流の時間領域での解析手法について説明する。流体場温度場電磁場の方程式をGalerkin有限要素法で定式化する。流体場の計算方法は、線形のスカラー形状関数を基底関数としたGSMAC法である。温度場の計算方法は、関数と微係数が連続となるスカラー形状関数を基底関数としたCIP法である。電磁場の計算方法は、微分形式に基づいたベクトル形状関数を基底関数とした辺有限要素法Aφ法である。流体場温度場電磁場よりなる場の支配方程式について、以下に説明する。場は流体場温度場電磁場よりなる。流体場の支配方程式は、ρ0=-gradpdivJ×B-ρ0gdivv=0である。温度場の支配方程式は、ρcv=divsTsT=2ηddJ2σである。電磁場の支配方程式は、rotE=-rotH=JdivB=0divJ=0である。流体場は、Boussinesq近似された運動方程式ρ0=-gradpdivFemρ0gと非圧縮連続の方程式divv=0により表現される。ここで、ρ0は基準密度スカラー、vは速度ベクトル、pは圧力スカラー、ηは粘性係数スカラー、dは変形速度テンソル、Femは電磁力ベクトルで、Maxwellの応力テンソルTemの発散として計算される。即ち、Fem=divTem=divIである。ここで、Bは磁束密度ベクトル、Hは磁場ベクトル、Iは恒等テンソルである。BHは、BiHjをij要素とするテンソルである。テンソルの発散は、テンソル場の偏微分係数を縮約して得られるベクトルである。ddは、tr=Σijdij2である。Boussinesq近似において密度は、ρ=ρ0と表せる。ここで、βは体膨張率スカラー、ΔT=T-T0は基準温度T0からの温度差スカラー、gは重力加速度ベクトルである。温度場は、エネルギー方程式ρcv=divsTにより表現される。ここで、cvは定容比熱スカラー、kは熱伝導率スカラー、sTは発熱項である。発熱項は、粘性発熱項2ηddとJoule発熱項J2σの和である。ここで、Jは電流密度ベクトル、σは電気伝導度スカラーである。電磁場は、準定常近似のMaxwellの方程式rotE=-rotH=JdivB=0により表現される。電流の発散が零となる両立性の条件はdivJ=0である。Aφ法で定式化するとB=rotAE=--gradφ=-となる。ここで、φは電気スカラーポテンシャル、Aは磁気ベクトルポテンシャルである。Vは対称化ポテンシャルで、φ=∂V∂tである。このとき、divB=0rotE=-は自動的に満たされる。よって、解くべき方程式はrotH=Jのみである。これより、誘導方程式rotrotB=--rotが導ける。μは透磁率スカラーである。これをAV法で表すとrot=-σ×vとなる。また、両立性の条件はdiv×v=0となる。ここで、V=0のゲージを選択すればA法となる。静止導体の場合、速度はv=0である。境界条件はn×A=0V=constである。これらの式を辺有限要素法で解く。次に、接続行列について説明する。接続行列は特異行列で、勾配行列回転行列発散行列の3種類のものがある。∂e=Gn∂f=Re∂v=Dfこれらは、辺の境界が接点、面の境界が辺、体の境界が面であることを表している。ここで、接続行列の成分表示はG=G=Gen=GenR=R=Rfe=RfeD=D=Dvf=Dvfであり、これらの転置行列はGt=Gt=Gne=GneRt=RefDt=Dfvである。このように、行列の表示方法には色々ある。G、R、Dは成分が-1、0、1からなる接続行列である。このとき、∂∂=0→RG=0、DR=0が成立する。これは、外微分演算dd=0→=0、=0に対応する。このとき、次の対応関係B=rotA←→B=RA←→Bf=RfeAeE=-gradφ-←→E=-Gφ-が成立する。RG=0の証明∂∂f=∂Re=R∂e=RGn=0∴RG=0B=RAの証明面fを通過する磁束を計算する。磁束密度を、面基底要素を用いて、B=ΣfMfBfと表す。このとき、直交性∫fMfndS=δffが成立する。よって、Bの法線成分∫fBndS=Bfは境界面上で連続である。Bの接線成分は一般に不連続である。つぎに、Bの法線成分をAの接線成分で表現する。Bf=∫fBndS=∫fndS=∫∂fAds=∫Ads積分範囲はΣeRfee=ΣeRfe∫eAds=ΣeRfeAeゆえに、Bf=ΣeRfeAe→B=RA←→B=rotAを得る。微分演算子と接続行列の関係をまとめて示す。微分演算子と接続行列は次の通りである。∂e=GnE=gradφEe=Gφn∂f=ReB=rotABf=RAe∂v=Dfρ=divBρv=DBfStokesの定理を示す。∫Ωdω=∫∂Ωωdと∂は双対である。Ωdω=∂ΩωEe=∫eEdr=∫∂eφ=Gen∫nφ=Genφn∴Ee=GenφnBf=∫fBndS=∫∂fAdr=Rfe∫eAdr=RfeAe∴Bf=RfeAeρv=∫vρdV=∫∂vBndS=Dvf∫fBndS=DvfBf∴ρv=DvfBfMaxwellの方程式と接続行列の関係をまとめて示す。基本場と物質場のMaxwellの方程式と接続行列は次の通りである。基本場基本空間物質場双対空間divB=0DB=0divD=ρeGtD=ρerotE=-rotH=JRE=-RtH=JB=rotAB=RAH=-gradφmH=-DtφmE=--gradφE=--Gφ次に、基底間の関係式について説明する。基底の間につぎの関係式が成立する。rotN=RtM→VrotN⊂VMgradN=GtN→VgradN⊂VMここで、VrotNは、辺基底の回転rotNeが張る空間である。VMは、面基底Mfが張る空間である。VgradNは、節点基底の勾配gradNnが張る空間である。VNは、辺基底Neが張る空間である。rotN=RtMを証明する。まずB=rotA=ΣeAerotNeである。一方、B=ΣfMfBf=ΣfMf∫fBndS=ΣfMf∫∂fAds=ΣfMf=ΣeAeである。両式を比較してrotNe=ΣfRfeMf←→rotN=RtMを得る。これは、VrotN⊂VMを意味する証明終わり。次に、接続行列のランクについて説明する。有限要素の総節点数をNd、総辺数をEd、総面数をFc、総要素数をVmとする。独立な辺の数は木の枝の数に等しく、Nd-1である。よって、∂e=GnrankG=Nd-1である。独立な面の数は補木の弦の数に等しく、Fc-1である。よって、∂f=RerankR=Fc-1である。そして、辺数の総数Edは、木の枝の数補木の弦の数=Nd-1Fc-1=Edとなる。ここで、単連結領域の場合、FcとEdとNdの間にEulerの関係式Fc-EdNd=2が成立する。図4aに、接続行列による変数の変換と形状関数の変換をまとめて示す。図4bに、関数空間の関係をまとめて示す。ここで、節点、辺、面、体の基底形状関数は、Nn=wnNe=weMf=wfNv=wvであり、基底空間はVgradN=gradW0⊂VN=W1VrotN=rotW1⊂VM=W2である。図5に、双対空間基本場と物質場と接続行列の関係式をまとめて示す。双対空間の接続行列は基本空間の転置行列である。ここで、物質場双対空間におけるgradrotdivの接続行列は、G=DtR=RtD=Gtとなる。次に、静磁場と辺有限要素法について説明する。磁場には静磁場と動電磁場の問題がある。静磁場は動電磁場の特別な場合であるが、動電磁場を理解するのに静磁場解析を理解しておくのが先決である。ここでは、静磁場の解析を、辺有限要素法を用いて行う。すなわち、ベクトル形状関数を用いた辺有限要素法で定式化する。磁場の支配方程式rot=Jに重み付き残さ法を適用する。ベクトル重み関数wを内積して領域積分すると∫ΩwrotdΩ=∫ΩwJdΩとなる。ここで、微分公式div×w=wrot-を用いる。部分積分して発散項が零になることを用いると、弱形式∫ΩdΩ=∫ΩwJdΩが導ける。辺要素基底ベクトル形状関数Neを適用し、AをA=ΣeNeAeと補間する。重み関数wに同じ辺要素基底を適用し、w=Neとすれば、Galerkin有限要素法となる。辺要素基底は直交性∫eeNeds=δeeを満足し、要素境界上でAの接線成分が連続である。すなわち、∫eAds=Aeとなる。Aの法線成分は一般に不連続である。有限要素方程式について説明する。弱形式で定式化された方程式にAとwの辺要素基底で表された式を代入すると有限要素法の離散式ΣeAe∫ΩdΩ=∫ΩNeJdΩが導ける。ここで、和はすべての辺e=12...Edの和である。Edは全辺数である。ここで、rotN=RtMを用いる。上式より、行列方程式KA=Jを得る。ただし、Kee=∫ΩdΩJe=∫ΩNeJdΩである。次に、有限要素行列の性質について説明する。有限要素行列は特異行列で正則行列でない。すなわち、K=RtνR:e×e行列Rt:e×f行列R:f×e行列ν:対称正則f×f行列νff=∫ΩνMfMfdΩである。行列Kのdimとrankは異なる。すなわち、dimK=EdrankK=rankR=Ed-Nd1である。よって、行列Kは特異行列である。次に、木補木ゲージによる正則化について説明する。木の自由度を消去して、Kを正則化してK’とする。補木の弦の数は独立なループの個数に等しい。木の自由度を消去して正則化すると、最小固有値が減少し、条件数が悪化する場合が極めて多い。よって、特異行列をそのままICCG法で解いた方が良い。誘導方程式と辺有限要素法について説明する。動電磁場の支配方程式は誘導方程式である。この誘導方程式をAφ法で記述する。この場合、電流は励振電流ソース電流Jsと誘導電流Jiの和としてJ=JsJiと表せる。そして、誘導電流はJi=σ=-σ×vgradφとなる。Aφ法の基礎方程式について説明する。Aφ法の基礎方程式と両立性の条件はつぎのとおりである。rot=JsJi=Js-σ×vgradφdiv×v=0この方程式を基礎に色々な変形が可能である。AV法の基礎方程式について説明する。Aφ法の基礎方程式において、添え字を省略してJs=Jφ=∂V∂tと置くと、AV法の基礎方程式rotσ∂t×v=Jdiv∂t×v=0が得られる。ここで、Aを辺要素で補間し、Vを節点要素で補間するとA=ΣeNeAeV=ΣnNnVnとなる。速度v=0の場合を考える。これより弱形式∫ΩνdΩjω∫ΩσNedΩ=∫ΩNeJdΩjω∫ΩσgradNndΩ=0が導ける。よって、AV法の有限要素方程式は図6bに示すようになる。さらに、σ、νは正定値行列で、その成分は、σee=∫ΩσNeNedΩνff=∫ΩνMfMfdΩである。A法の基礎方程式について説明する。AV法の基礎方程式で、V=0のゲージを採用すればA法の行列方程式KA=Jが得られる。ここで、K=RtνRjωσである。この式は単純化されているが、収束性はA法よりもAV法の方が一般に勝っている。時間領域積分について説明する。A法の基礎方程式としてrotσ×v=0div×v=0を考える。弱形式で定式化すると∫ΩνdΩ∫ΩσNe×vdΩ=0となる。よって、離散化式を得る。ここで、K=RtνR-AdAd=∫ΩσNedΩσ=∫ΩσNeNedΩである。時間積分には一般化台形公式θKn1σΔtAn1=-Kn-σΔtAnを用いる。ここで、θは混合パラメータで、この方程式は陰解法である。両立性の条件は-∫ΩσgradNn×vdΩ=∫ΩσgradNnJidΩ=0となる。この式は、divJi=0と同等である。よって、大域的には境界において電流が流出しない条件nJi=0と同等となる。よって、この条件を満足しないと解は不定となる。流体場には温度差が浮力として作用し、Lorentz力として電磁力が働く。温度場には熱源項として、粘性散逸の発熱項とJoule発熱項がある。電磁場には導体の運動による誘導起電力がある。すなわち、流体場と温度場と電磁場が互いに連成する場合を考える。そして、これら3つの場に、次に示す3種類の有限要素法を適用する。流体場には、スカラー形状関数を用い、GSMACFEMを適用する。温度場には、スカラー形状関数を用い、CIP-FEMを適用する。電磁場には、ベクトル形状関数を用い、辺有限要素法Aφ法を適用する。これにより、高Reynolds数、高Prandtl数、高磁気Reynolds数に安定な解が得られる。有限要素法に用いられる基底関数形状関数には、節点基底関数、辺基底関数、面基底関数、体基底関数の4種類がある。節点基底関数と体基底関数はスカラー関数であるが、辺基底関数と面基底関数はベクトル関数である。これらの基底関数は、図7に示すように、微分形式と密接の関係している。そして、積分値が未知変数となっている。特に電磁場においては、Gaussの定理やStokesの定理が重要で、これらの定理を辺有限要素法はうまく表現している。次に、図8を参照しながら、電磁熱流体が空気中におかれた場合の、電磁熱流体解析の計算例を説明する。電磁熱流体としては、金属微粒子を含有する溶融ガラスとする。また空気中には電磁場が存在するとする。図8aに解析モデルを示す。大立方体の中心の小立方体が電磁熱流体の領域である。x方向に電位差が生じるように、小立方体の端面の電位を設定している。-z方向に重力がかかり、z方向に磁場がある。図8bに主な解析パラメータを示す。図8cに解析メッシュを示す。解析メッシュは6面体要素に分割されている。図9a、b、cに計算開始から約2000秒後の計算結果を示す。図9aは流体場の、図9bは温度場の解析結果である。また、電磁場の解析の結果得られた例として磁束密度の分布を図9cに示す。白黒でかつ簡易的に表示してあるので明瞭ではないが、小立方体上部でx断面において流体の速度が大きくなっているところがあることが見て取れる。また温度は、小立方体最上部よりも内部のほうが高温になっていることがわかる。磁束密度は小立方体の中心部にはほとんど無く周辺部で高くなっていることがわかる。上記のように、本発明の実施例では、電磁熱流体解析装置を、GSMACFEM流体解析手段で微小時間の流体の力学的変化を解析してデータ変換し、CIPFEM流体解析手段で流体の熱力学的変化を解析してデータ変換し、AφFEM流体解析手段で流体の電磁気学的変化を解析してデータ変換し、これを繰り返して電磁熱流体の動きを解析する構成としたので、電磁場温度場流体場が同時に解けて、高速高精度で汎用性が高い解析システムが実現できる。