(57)【要約】【課題】保存形式の流体方程式を解く際、特にファンデルワールス状態方程式を用いて相変化計算を行う際には、特性速度が虚数となる領域があるため、ＦＤＳ、ＦＶＳといった解法は使えない。またセントラルスキームでは主格子、双対格子と保存量の補間（もしくは積分）のような操作を行うため誤差が大きくなる。またＡＵＬＦスキームは時間伴にずれた定義点を持つ主格子と双対格子の初期値を辻褄が合うように設定することが困難であった。【解決手段】主格子から双対格子の値を求める工程Ａと、流速に基づく風上差分法によって流束の移流部分の時間発展を解く工程Ｂと、工程Ａによって求めた物理量の中心差分によって流束の圧力が陽に現れる部分の時間発展を解く工程Ｃとを備えることを特徴とする。よって双対格子の値は仮の値としての役割しかないため誤差の発生が抑えられるし、双対格子の初期値を設定する必要もない。【選択図】図１
2009110082.tif 000002
【特許請求の範囲】
【請求項１】
保存形式の流体方程式を、空間を構造格子で離散化し、格子の格子線で囲まれる最小の空間をセルと定義し、時間発展の時間刻み幅をΔｔとし、該構造格子の保存量の定義点に対して時間軸ではΔｔ／２だけ進み、空間座標では保存量の定義点がセル中心の時には格子点を、格子点の時にはセル中心を双対格子の保存量の定義点として計算する数値計算方法であって、双対格子の保存量を求める工程Ａと、流速に基づく風上差分法によって流束の移流部分の時間発展を解く工程Ｂと、工程Ａによって求めた物理量の中心差分によって流束の圧力が陽に現れる部分の時間発展を解く工程Ｃとを備えることを特徴とする流体の数値計算方法。
【請求項２】
請求項１に記載の流体の数値計算方法であって、上記工程Ｂの移流させる速度として、上記工程Ａで求めた流速を使うことを特徴とする流体の数値計算方法。
【請求項３】
請求項１に記載の流体の数値計算方法であって、上記工程ＡにおいてＬａｘ−Ｆｒｉｅｄｒｉｃｈｓ法によって保存量を求めることを特徴とする流体の数値計算方法。
【請求項４】
請求項１に記載の流体の数値計算方法であって、上記工程Ｂの移流させる速度として、上記工程Ａで求めた流速を使い、上記工程ＡにおいてＬａｘ−Ｆｒｉｅｄｒｉｃｈｓ法によって保存量を求めることを特徴とする流体の数値計算方法。
【請求項５】
コンピュータに請求項１乃至４の何れか１項に記載の流体の数値計算方法を実行させることを特徴とするプログラム。
【請求項６】
請求項５に記載のプログラムを記録したことを特徴とするコンピュータ読み取り可能な記録媒体。
【請求項７】
保存形式の流体方程式を、空間を構造格子で離散化し、格子の格子線で囲まれる最小の空間をセルと定義し、時間発展の時間刻み幅をΔｔとして計算する流体の計算装置であって、時間軸ではΔｔ／２だけ進み、空間座標では保存量の定義点がセル中心の時には格子点、格子点の時にはセル中心の保存量を求める手段Ａと、流速に基づく風上差分法によって流束の移流部分の時間発展を解く手段Ｂと、手段Ａによって求めた物理量の中心差分によって流束の圧力が陽に現れる部分の時間発展を解く手段Ｃとを備えることを特徴とする流体の計算装置。
【発明の詳細な説明】
【技術分野】
【０００１】
本発明は流体の数値計算方法、プログラム、プログラムを記憶した記録媒体及び装置に関するものである。
【背景技術】
【０００２】
流体の流れを扱う学問、流体力学は古くから計算されているが、その基礎方程式は非線形な連立方程式であり、限られた場合を除いて解析解を得ることは困難である。そのためコンピュータによる数値計算によって流れ場の近似解を得ようとする試みがなされてきた。
【０００３】
図２にコンピュータによって数値計算プログラムを実行する際の装置の構成を示す。
【０００４】
２０２はＣＰＵ群であり、ここで演算を行う。
【０００５】
２０３は入力装置群であり、マウス、キーボード、スキャナなどの情報入力装置を指す。
【０００６】
２０４は一次記憶装置群であり、メモリなどの一時的な記憶装置を指す。
【０００７】
２０５は出力装置群であり、ディスプレイ、プリンタなどの情報出力装置を指す。
【０００８】
２０６は通信装置群であり、ＬＡＮボートなど通信装置をさす。
【０００９】
２０７は二次記憶装置群であり、ハードディスク、フレキシブルディスク、ＣＤ―ＲＯＭ、ＤＶＤ―ＲＯＭなどのデータ保存のための装置を指す。
【００１０】
２０８はバスラインであり、コンピュータ内の装置群間の情報を伝達する役割を持つ。
【００１１】
２０９はＬＡＮケーブルであり、異なるコンピュータ間の通信装置に接続し、情報を伝達する。
【００１２】
以下で述べる流体の数値計算方法はプログラムとして、二次記憶装置群２０７に保存されている。実行時には一次記憶装置群２０４に読み込まれ、２０３、２０９の外部からの指令によって、２０２のＣＰＵ群で実行される。
【００１３】
解きたい問題の種類によって解くべき基礎方程式の形は異なる。気体−液体間の相変化について解きたい場合には状態方程式としてファンデルワールス状態方程式を用いた上で、質量保存、運動量保存及びエネルギー保存の方程式を基礎方程式として選択する方法がある。ファンデルワールス状態方程式はその密度の高い部分と低い部分がそれぞれ液体と気体を表現している。従って、ファンデルワールス状態方程式に従って密度を変化させていけば自然と相変化を表現することができる。
【００１４】
しかし気体から液体、液体から気体へと状態が遷移する部分では情報の伝播速度とも言うべき音速が虚数となってしまう領域が存在する。そのため音速を直接用いるＦＤＳ(Flux Difference Splitting)、ＦＶＳ(Flux Vector Splitting)といった手法を適用することは難しい。
【００１５】
よって音速を直接用いない手法を用いることが考えられる。そのうち最も近接する技術は非特許文献１に示されるセントラルスキームであり、特許文献１に示されるスキーム（ＡＵＬＦスキーム）である。
【００１６】
上記２種のスキーム及び本発明の説明のために、一次元の格子幅がΔｘで一様な格子を用いて説明を行う。また質量保存、運動量保存及びエネルギー保存の方程式として以下の（式１）、（式２）、（式３）を用いて説明を行う。ここでｔは時間、ｘは空間座標、Ｑは保存量、Ｅは流束、ρは密度、ｕは速度、ｅは単位質量あたりの全エネルギー、ｐは圧力である。
【００１７】
【数１】
2009110082.tif 000003
【００１８】
【数２】
2009110082.tif 000004
【００１９】
【数３】
2009110082.tif 000005
上記の例は本発明の説明を簡単にするためであり、本発明を制限するようなものではない。二次元、三次元空間への拡張は容易である。また簡単な変換で格子幅が一様な直交格子に変換できるような格子、例えば不等間隔格子、斜交格子、境界適合構造格子への適用も容易である。方程式についてもその他の項、例えば粘性項、表面張力項、熱伝導項が加わっても本発明は適用可能である。（式１）の右辺にソース項、例えば発熱項が加わっても本発明が適用可能である。また保存量Ｑの選び方には任意性がある。（式４）から（式７）に示されるような内部エネルギーを用いる手法や（式４）及び、（式８）から（式１０）のように温度を用いる手法もあり、いずれも本発明が適用可能である。尚これらの式の中でεは単位質量あたりの内部エネルギー、Ｔは温度、Ｃｖは定積比熱である。
【００２０】
【数４】
2009110082.tif 000006
【００２１】
【数５】
2009110082.tif 000007
【００２２】
【数６】
2009110082.tif 000008
【００２３】
【数７】
2009110082.tif 000009
【００２４】
【数８】
2009110082.tif 000010
【００２５】
【数９】
2009110082.tif 000011
【００２６】
【数１０】
2009110082.tif 000012
図３に説明のために作成した格子を示す。格子の格子線が囲む最小の空間をセルと定義する。図中の３０１で示す空間がセルである。以下の説明ではセル中心に保存量Ｑを定義する。３０２がセル中心である。図３で作られた格子を以下主格子と呼ぶ。
【００２７】
セントラルスキーム及び特許文献１に示されるＡＵＬＦスキームでは双対格子を用いる。主格子における時間刻み幅をΔｔとする時、双対格子は主格子に対して空間座標でΔｘ／２だけずれており、時間軸でΔｔ／２だけずれるように作成する。図４の４０１が双対格子のセル中心である。なお４０２が主格子のセル中心である。
【００２８】
また以下の説明を簡単にするために、主格子のセル中心をｘ軸の負の方向から順番に整数で番号をつける。整数ｉを用いて、主格子のｉ番目のセル中心の空間座標をｘｉと下付き添え字を用いて示すこととする。また物理量がｉ番目の主格子のセル中心に定義されていることも下付き添え字で示すこととする。例えばＱがｉ番目の主格子のセル中心に定義されていればＱｉと書く事とする。時間方向も初期時刻ｔinitからの経過時間ｔが整数ｎを用いてｔ＝ｎΔｔと書けることから、ｎを上付き添え字にして時刻を表現することとする。つまりｔｎ＝ｎΔｔである。同様に時刻ｔｎの物理量であることを上付き添え字を用いて示すこととし、例えば時刻ｔｎのＱの値はＱｎと書く。
【００２９】
また添え字が半整数の時にはｘｉ＋１／２：＝ｘｉ＋Δｘ／２及びｔｎ＋１／２：＝ｔｎ＋Δｔ／２のように定義される。物理量の添え字も同様である。
【００３０】
セントラルスキームの最も低次な形は次のような形である。
【００３１】
【数１１】
2009110082.tif 000013
【００３２】
【数１２】
2009110082.tif 000014
そして、Ｌａｘ−Ｆｒｉｅｄｒｉｃｈｓ法を用いて上式を解くことができる。
【００３３】
図５にセントラルスキームで求める保存量Ｑの場所を矢印でつないだ図を示す。矢印の元の値を用いて矢印の先の値を求める。
【００３４】
またＡＵＬＦスキームは保存量Ｑを主格子と双対格子で別々に定義し、
【００３５】
【数１３】
2009110082.tif 000015
【００３６】
【数１４】
2009110082.tif 000016
とし、Ｅａｄｖによる時間発展は風上差分を、Ｅｐによる時間発展は双対格子の値を用いた蛙飛び差分を用いている。図６にＡＵＬＦスキームを用いて求められる保存量Ｑの場所を矢印でつないだ図を示す。主に矢印の元の値を用いて矢印の先の値を求める。
【非特許文献１】H. Nessyahu, E.Tadmor,“Non-oscillatory Central Differencing for Hyperbolic Conservation Laws”Journal of Computational Physics 87 (1990) 408-463.
【特許文献１】特開２００５−０５００６３号公報
【発明の開示】
【発明が解決しようとする課題】
【００３７】
セントラルスキームでは図５から理解できるように、保存量Ｑの求められる場所が主格子、双対格子と移り変わっていく過程で補間のような操作を繰り返すため大きな誤差を生む。もう少し具体的な例をあげよう。仮想的にｕ＝０の状態を考える。この時、流束が０となり、本来時間発展をしても保存量Ｑの値は変わらないはずである。セントラルスキームは（式１５）、（式１６）のように時間発展を行う。
【００３８】
【数１５】
2009110082.tif 000017
【００３９】
【数１６】
2009110082.tif 000018
（式１６）からわかるように、数値拡散が生じてしまう。この精度はセル内分布の再構築精度に依存するが、セントラルスキームを用いた場合、誤差の発生を避けることは難しい。しかもこの誤差は累積するために時間進行に従って誤差はどんどんと大きくなる。
【００４０】
ＡＵＬＦスキームは主格子と双対格子とで別々に保存量Ｑを定義している。従って図６からも理解できるように、主格子と双対格子のそれぞれ、即ち時間軸も空間座標も異なった定義点の初期条件を辻褄が合うように与えなければならない。通常初期条件は人為的に与えるが、時間軸が異なっている定義点の初期条件を辻褄が合うように決めるというのはもはや基礎方程式を解く必要があり、人為的に決めることは困難である。
【００４１】
また主格子と双対格子とで別々に保存量Ｑを定義している上に、互いの値を用いる機会が流束を介してのみ行われる。従って、両者の格子の保存量Ｑが空間的には近接しているにも関わらず乖離してしまう可能性がある。
【００４２】
また数値計算を実行する際には、計算の途中から再開できるように保存量Ｑの値を含む変数の値を二次記憶装置に書き込んで保存することがよく行われるが、ＡＵＬＦスキームは主格子と双対格子の両者の保存量Ｑを保存しておかなければならない。
【課題を解決するための手段】
【００４３】
上記の課題を解決するための数値計算方法は、
保存形式の流体方程式を、空間を構造格子で離散化し、格子の格子線で囲まれる最小の空間をセルと定義し、時間発展の時間刻み幅をΔｔとし、該構造格子の保存量の定義点に対して時間軸ではΔｔ／２だけ進み、空間座標では保存量の定義点がセル中心の時には格子点を、格子点の時にはセル中心を双対格子の保存量の定義点として計算する数値計算方法であって、双対格子の保存量を求める工程Ａと、流速に基づく風上差分法によって流束の移流部分の時間発展を解く工程Ｂと、工程Ａによって求めた物理量の中心差分によって流束の圧力が陽に現れる部分の時間発展を解く工程Ｃとを備えることを特徴とする流体の数値計算方法である。
【００４４】
また、上記課題を解決するための計算装置は、保存形式の流体方程式を、空間を構造格子で離散化し、格子の格子線で囲まれる最小の空間をセルと定義し、時間発展の時間刻み幅をΔｔとして計算する流体の計算装置であって、時間軸ではΔｔ／２だけ進み、空間座標では保存量の定義点がセル中心の時には格子点、格子点の時にはセル中心の保存量を求める手段Ａと、流速に基づく風上差分法によって流束の移流部分の時間発展を解く手段Ｂと、手段Ａによって求めた物理量の中心差分によって流束の圧力が陽に現れる部分の時間発展を解く手段Ｃとを備えることを特徴とする流体の計算装置である。
【発明の効果】
【００４５】
本発明の双対格子の値は仮の値としての役割しかない。よってセントラルスキームとは異なり、あくまでも主格子の保存量Ｑから主格子の保存量Ｑへと時間発展を行う。よって数値拡散の累積が抑えられ、誤差の発生を抑えることができる。
【００４６】
また本発明ではＡＵＬＦスキームとは異なり、双対格子の値を常に主格子から求めるため、双対格子の初期条件を与える必要がない。よって双対格子との辻褄について気にすることなく、主格子の初期条件を設定することができる。また双対格子の値を常に主格子から求めるため、空間的に近接する主格子と双対格子の保存量の値が乖離することがない。また計算の途中から再開できるように保存量Ｑの値を含む変数の値を保存する際に、保存量Ｑについては主格子の値のみ保存すればよいためＡＵＬＦスキームに比べて少ない二次記憶装置容量で良い。
【００４７】
図１に本発明を一次元の格子幅がΔｘで一様な格子へと適用した場合の、求められる保存量Ｑの場所を矢印でつないだ図を示す。主に矢印の元の値を用いて矢印の先の値を求める。
【００４８】
繰り返しになるが、二次元、三次元空間への拡張は容易である。また簡単な変換で格子幅が一様な直交格子に変換できるような格子、例えば不等間隔格子、直交格子、一般座標格子への適用も容易である。方程式についても（式１）、（式２）、（式３）のような方程式にも適用が可能な他、その他の項、例えば粘性項、表面張力項、熱伝導項が加わっても本発明が適用可能である。（式１）の右辺にソース項、例えば発熱項が加わっても本発明が適用可能である。また保存量Ｑの選び方には任意性がある。（式４）から（式７）に示されるような内部エネルギーを用いる手法や（式４）及び、（式８）から（式１０）のように温度を用いる手法もあり、いずれも本発明が適用可能である。
【００４９】
また図７に示すような直交格子をAdaptive Mesh Refinementした格子（以下ＡＭＲ格子）への適用を考えよう。ＡＭＲ格子では保存量Ｑをセル中心以外（例えば格子点）に定義する場合、保存性を保つためには複雑な処理が必要となる。それはコントロールボリュームの作り方が難しかったり、流束の受け渡しが難しかったりすることが原因である。本発明ではセントラルスキームやＡＵＬＦスキームとは異なり、保存量Ｑが主格子のセル中心でのみ定義されており、流束の受け渡しも常に主格子のセル同士で行われる。よってＡＭＲ格子へ保存量Ｑの保存性を確保したまま、簡単に適用することができる。
【００５０】
また以下の実施例１でも述べるが、高次化は容易である。工程Ａを高次化するための手順はセントラルスキーム関連の手法を用いれば良い。工程Ｂを高次化するためには高次の風上移流スキームを用いれば良い。
【発明を実施するための最良の形態】
【００５１】
以下、図面を参照して本発明を実施するための最良の形態に関する実施例について詳細に説明する。
【実施例１】
【００５２】
説明の簡単のために最も低次な形の説明を行う。（式１）、（式２）、（式３）で示されるような保存形式の流体方程式を解くこととし、流束を（式１３）、（式１４）のように流束の移流部分Ｅａｄｖと、流束の圧力が陽に現れる部分Ｅｐとに分離する。
【００５３】
まずはＥｐによる保存量Ｑの変化分を考える。最も低次な形では工程Ａとして（式１１）のように、双対格子のセル中心でＱｉ＋１／２ｎ＋１／２を計算する。Ｑｉ＋１／２ｎ＋１／２と状態方程式からＥｐｉ＋１／２ｎ＋１／２を求めることができるので、工程Ｃとして（式１７）のように計算することによってＥｐによる保存量Ｑの変化分ΔＱｐを計算することができる。
【００５４】
【数１７】
2009110082.tif 000019
次に流速に基づく風上差分法によって流束の移流部分の時間発展を解く工程ＢのＥａｄｖによる保存量Ｑの変化分を考える。最も低次な形は一次の風上差分であり、（式１８）、（式１９）のようにＥａｄｖを求める。但しここで用いたｕｉ＋１／２ｎ＋１／２は上記Ｑｉ＋１／２ｎ＋１／２から求める。
【００５５】
【数１８】
2009110082.tif 000020
【００５６】
【数１９】
2009110082.tif 000021
これを用いて、（式２０）のようにＥａｄｖによる保存量Ｑの変化分ΔＱａｄｖを計算することができる。
【００５７】
【数２０】
2009110082.tif 000022
以上よりＱｉｎ＋１は（式２３）のように求めることができる。
【００５８】
【数２１】
2009110082.tif 000023
以上の説明では最も低次な形で説明を行った。工程Ａを高次化するためには（式１１）を高次化すれば良い。その手順はセントラルスキーム関連の手法を用いれば良い。工程Ｂを高次化するためには（式１９）の部分を高次化すれば良い。（式１９）では一次の風上差分を用いているため、セル内の保存量Ｑの分布を０次の関数（一様分布）で再構築した。このセル内の分布をより高い精度の関数を用いて再構築すれば高次化が可能である。
【００５９】
理想気体の衝撃波管問題を用いて計算結果を示す。計算空間ｘ［０：１］であり、ｘ＜０．５にρ＝１．０、ｐ＝１．０なる理想気体を配置し、ｘ≧０．５にρ＝０．１２５、ｐ＝０．１なる理想気体を配置した。１００個の格子で分割し、Δｔ＝２．０２２×１０−３として、ｔ＝０からｔ＝０．１４１５４まで計算を行った。
【００６０】
図８に本発明による計算結果を示す。横軸が座標ｘ、縦軸がρである。衝撃波面に少しオーバーシュートが見られるものの正しく計算できていることがわかる。図９に同様の計算条件でセントラルスキーム、即ち（式１１）、（式１２）を用いた結果を点線で示す。横軸が座標ｘ、縦軸がρである。実線で示される本発明のスキームに比べて数値拡散が大きいことが分かり、本発明を用いることで、数値拡散を抑えることができていることがわかる。
【実施例２】
【００６１】
実施例２では本発明による数値計算手法を用いて、気液相変化計算を行う。一次元空間に液体を満たし、端部を熱することで気体へと状態変化する計算を行う。
【００６２】
状態方程式としてファンデルワールス状態方程式を用い、その定数には一般的な水の値を用いた。
【００６３】
質量保存、運動量保存及びエネルギー保存の方程式は相変化計算のために（式２２）、（式２３）、（式２４）のように変更を行った。保存量Ｑは（式２）と同様。αは表面張力の大きさをコントロールする係数を示し、Ｑｅは加熱によるソース項を示す。
【００６４】
【数２２】
2009110082.tif 000024
【００６５】
【数２３】
2009110082.tif 000025
【００６６】
【数２４】
2009110082.tif 000026
またΠｘｘ、Ｓｘｘ及びｑは（式２５）、（式２６）、（式２７）のように定義される。
【００６７】
【数２５】
2009110082.tif 000027
【００６８】
【数２６】
2009110082.tif 000028
【００６９】
【数２７】
2009110082.tif 000029
ここで、ηはせん断粘性係数、ζは体積粘性係数、λは熱伝導率を示す。
【００７０】
上記方程式を本発明を用いて解く。流束を（式２８）、（式２９）のように流束の移流部分Ｅａｄｖと、流束の圧力が陽に現れる部分Ｅｐと、その他の部分Ｅｏｔｈｅｒｓに分離する。
【００７１】
【数２８】
2009110082.tif 000030
【００７２】
【数２９】
2009110082.tif 000031
（式２９）を見れば分かるようにＥａｄｖとＥｐは（式１４）と変わらない。よって本発明の主眼たる部分は実施例１と同様の方法で解く事ができ、ＥｏｔｈｅｒｓとＳｏについての扱いのみが異なる。ＥｏｔｈｅｒｓとＳｏは陽解法で扱うものとし、中心差分を用いて離散化する。
【００７３】
計算空間はｘ［０：０．１］［μｍ］であり、２００個の格子で分割した。初期はρ＝４４５．４５［ｋｇ／ｍ３］、Ｔ＝４００［Ｋ］なる水で満たされている。またその他の物性値は、η＝１．８×１０−５［Ｐａ  ｓｅｃ］、ζ＝０、λ＝５．０×１０−２［J／（ｍ  ｓｅｃ  Ｋ）］、α＝１．０×１０−１５［Ｊ  ｍ５／ｋｇ２］で定数とする。
【００７４】
Δｔ＝１．０×１０−１４［ｓｅｃ］として、ｔ＝０からｔ＝２．５×１０−９［ｓｅｃ］まで計算を行った。Ｑｅは時刻ｔ＝０からｔ＝１．５×１０−９［ｓｅｃ］まで、ｘ［０：１．５×１０−３］［μｍ］の範囲で、Ｑｅ＝７．６×１０１７［Ｊ／（ｍ３  ｓｅｃ）］の大きさで与えるものとした。
【００７５】
境界条件はｘ＝０では対称境界条件、ｘ＝０．１付近ではρとＴは初期条件で固定し、速度は０次の外挿によって求めた。
【００７６】
図１０に本発明による計算結果を示す。横軸が座標ｘ、縦軸がρである。時刻ｔ＝０が実線で示されており、時間が経過するにつれて点線から長破線へと分布が変化する。このように加熱によって密度が小さくなり、やがて気体へと相変化を起こしていることが分かる。一方図１１に同様の計算条件でセントラルスキーム、即ち（式１１）、（式１２）を用いた結果を示す。セントラルスキームでは数値拡散が大きく発泡しない。このように本発明を用いることで気液相変化の計算が可能である。
【産業上の利用可能性】
【００７７】
流体の数値計算を行うことができる。とりわけ相変化計算を精度良く行うことができる。
【図面の簡単な説明】
【００７８】
【図１】本発明を一次元の格子幅がΔｘで一様な格子へと適用した場合の、求められる保存量Ｑの場所を矢印でつないだ図
【図２】本発明を実施可能な装置の図
【図３】一次元の一様な格子幅の格子
【図４】主格子と双対格子
【図５】セントラルスキームで求める保存量Ｑの場所を矢印でつないだ図
【図６】ＡＵＬＦスキームを用いて求める保存量Ｑの場所を矢印でつないだ図
【図７】直交格子をAdaptive Mesh Refinementした格子
【図８】実施例１における本発明を用いた計算結果
【図９】実施例１におけるセントラルスキームとの比較
【図１０】実施例２における本発明を用いた結果
【図１１】実施例２におけるセントラルスキームを用いた結果
【符号の説明】
【００７９】
２０２  ＣＰＵ
２０３  入力装置
２０４  一次記憶装置
２０５  出力装置
２０６  通信装置
２０７  二次記憶装置
２０８  バスライン
２０９  ＬＡＮケーブル
３０１  主格子のセル
３０２  主格子のセル中心
４０１  双対格子のセル中心
４０２  主格子のセル中心
【図１】
2009110082.tif 000032
【図２】
2009110082.tif 000033
【図３】
2009110082.tif 000034
【図４】
2009110082.tif 000035
【図５】
2009110082.tif 000036
【図６】
2009110082.tif 000037
【図７】
2009110082.tif 000038
【図８】
2009110082.tif 000039
【図９】
2009110082.tif 000040
【図１０】
2009110082.tif 000041
【図１１】
2009110082.tif 000042