(57)【要約】【課題】流体シミュレータの作成に当たって、気液の界面形状をシャープに保ったまま移動計算を行うために、格子点上に離散的に定義された離散関数から実一次元上の（不連続）関数を構成する補間アルゴリズムを提供する。【解決手段】１次元実領域上に構成された１次元構造格子と、該１次元構造格子の各胞中心において値が定まる事によって定義される該１次元構造格子上の関数ｆとが用意されたとき、該関数ｆと該格子の幾何を用いることにより、前記１次元実領域上の関数Ｆを定義する補間方法において、着目する格子に対して、該格子上の点を前進差分と後進差分の符号が異なるときには傾きをゼロとし、符号が同じ場合は、その絶対値の小さい方の２倍を傾きに選び、中心での値として前記着目する胞の中心の値ｆを持ち、該傾きを持つ１次関数によって、該胞で定まる１次元実領域の部分領域上での関数Ｆの値を定める。【選択図】　　　　図１
2004062749.tif 000002
【特許請求の範囲】
【請求項１】
１次元実領域上に構成された１次元構造格子と、該１次元構造格子の各胞中心において値が定まる事によって定義される該１次元構造格子上の関数ｆとが用意されたとき、該関数ｆと該格子の幾何を用いることにより、前記１次元実領域上の関数Ｆを定義する補間方法において、
着目する胞において、該関数ｆの前進差分と後進差分の符号が異なるときには傾きをゼロとし、該前進差分と後進差分の符号が同じ場合は、その絶対値の小さい方の２倍を傾きに選び、
中心での値として前記着目する胞の中心の値ｆを持ち、該傾きを持つ１次関数によって、該胞で定まる１次元実領域の部分領域上での該関数Ｆの値を定めることを特徴とする補間方法。
【請求項２】
移流型微分方程式の数値解法における補間方法として用いることを特徴とする請求項１記載の補間方法。
【請求項３】
２次元実領域上の構成された２次元構造格子と、該２次元構造格子の各胞中心において値が定まる事によって定義される該２次元構造格子上の関数ｆとが用意されたとき、該関数ｆと該格子の幾何を用いることにより、前記２次元実領域上の関数Ｆを定義する補間方法において、
着目する胞を胞甲とし、
胞甲の各辺を東西あるいは南北にそれぞれ合わせた際に東西南北方向に隣合う４胞に対して、前記胞甲の中心点と各４胞の中心点との前記関数ｆの片側差分値の２倍の値をそれぞれＤｆ東、Ｄｆ西、Ｄｆ南、Ｄｆ北とし、前記Ｄｆ東と前記Ｄｆ西が異符号の場合はゼロ、同符号の場合はその絶対値の小さいものをＤｆ東西とし、前記Ｄｆ南と前記Ｄｆ北が異符号の場合はゼロ、同符号の場合はその絶対値の小さいものをＤｆ南北とし、前記Ｄｆ東西と前記Ｄｆ南北の傾きを有し前記胞甲の中心において胞甲の中心での前記関数ｆの値を持つ平面を構成してそれを平面乙とし、
前記胞甲の北東方向の格子点における前記平面乙の値が胞甲の中心の値より大きい場合は、前記平面乙の該北東方向の格子点での値が、該格子点を共有する胞甲以外の３胞の中心でのｆの値を超えないようにＤｆ東西及びＤｆ南北を１以下の最大の定数により定数倍することにより前記平面乙を変更し、
胞甲の北東方向の格子点における平面乙の値が胞甲の中心の値より小さい場合は、前記平面乙の該北東方向の格子点での値が、該格子点を共有する胞甲以外の３胞の中心でのｆの値を下回らないように前記Ｄｆ東西及び前記Ｄｆ南北を１以下の最大の定数により定数倍することにより前記平面乙を変更し、
得られた平面乙に対して、該北東方向の格子点に対して行った操作を、北西方向の格子点、南西方向の格子点及び南東方向の格子点に対してそれぞれ行うことにより、該平面乙の傾きを変更させ、得られた平面を該胞甲内での前記補間関数Ｆとして定める事を特徴する補間方法。
【請求項４】
３次元実領域上の構成された３次元構造格子と、該３次元構造格子の各胞中心において値が定まる事によって定義される該３次元構造格子上の関数ｆとが用意されたとき、該関数ｆと該格子の幾何を用いることにより、前記３次元実領域上の関数Ｆを定義する補間方法において、
着目する胞を胞甲とし、
胞甲の各辺を東西、南北あるいは天地にそれぞれ合わせた際に、東西南北天地方向に隣合う６胞に対して、前記胞甲の中心点と各６胞の中心点との前記関数ｆの片側差分値の２倍の値をそれぞれＤｆ東、Ｄｆ西、Ｄｆ南、Ｄｆ北、Ｄｆ天、Ｄｆ地とし、前記Ｄｆ東と前記Ｄｆ西が異符号の場合はゼロ、同符号の場合はその絶対値の小さいものをＤｆ東西とし、前記Ｄｆ南と前記Ｄｆ北が異符号の場合はゼロ、同符号の場合はその絶対値の小さいものをＤｆ南北とし、前記Ｄｆ天と前記Ｄｆ地が異符号の場合はゼロ、同符号の場合はその絶対値の小さいものをＤｆ天地とし、前記Ｄｆ東西と前記Ｄｆ南北と前記Ｄｆ天地の傾きを有し前記胞甲の中心における関数ｆの値を通る補間関数の候補を構成してそれを平面乙とし、前記胞甲の北東天方向の格子点における平面乙の値が胞甲の中心の値より大きい場合は、平面乙の該北東天方向の格子点での値が、該格子点を共有する胞甲以外の７胞の中心でのｆ値を超えないように、前記Ｄｆ東西、前記Ｄｆ南北及び前記Ｄｆ天地を１以下の最大の定数により定数倍することにより平面乙を変更し、前記胞甲の北東天方向の格子点における平面乙の値が胞甲の中心の値より小さい場合は、平面乙の該北東天方向の格子点での値が、該格子点を共有する胞甲以外の７胞の中心でのｆ値を下回らないように前記Ｄｆ東西、前記Ｄｆ南北及び前記Ｄｆ天地を１以下の最大の定数により定数倍することにより平面乙を変更し、得られた補間関数の候補平面乙に対して、該北東天方向の格子点に対して行った操作を、北西天方向の格子点、南西天方向の格子点、南東天方向の格子点、北東地方向の格子点、北西地方向の格子点、南西地方向の格子点、及び南東地方向の格子点に対してそれぞれ行うことにより、平面乙の傾きを変更させ、得られた平面を該胞甲内での前記補間関数Ｆとして定める事を特徴とする補間方法。
【請求項５】
移流型微分方程式の数値解法における補間方法として用いることを特徴とする請求項４記載の補間方法。
【請求項６】
１次元実領域上に構成された１次元構造格子と、該１次元構造格子の各胞中心において値が定まる事によって定義される該１次元構造格子上の関数ｆとが用意されたとき、該関数ｆと該格子の幾何を用いることにより、前記１次元実領域上の関数Ｆを定義する補間方法を実行する実行装置において、
着目する胞に対して、該格子上の点を前進差分と後進差分の符号が異なるときには傾きをゼロに設定し、該格子上の点を前進差分と後進差分の符号が同じ場合は、その絶対値の小さい方の２倍を傾きに設定する手段と、
中心での値として前記着目する胞の中心の値ｆを持ち、該傾きを持つ１次関数によって、該胞で定まる１次元実領域の部分領域上での該関数Ｆの値を定める手段とを備えたことを特徴とする実行装置。
【請求項７】
２次元実領域上の構成された２次元構造格子と、該２次元構造格子の各胞中心において値が定まる事によって定義される該２次元構造格子上の関数ｆとが用意されたとき、該関数ｆと該格子の幾何を用いることにより、前記２次元実領域上の関数Ｆを定義する補間方法を実行する実行装置において、
着目する胞を胞甲とし、胞甲の各辺を東西あるいは南北にそれぞれ合わせた際に東西南北方向に隣合う４胞に対して、前記胞甲の中心点と各４胞の中心点との前記関数ｆの片側差分値の２倍の値をそれぞれＤｆ東、Ｄｆ西、Ｄｆ南、Ｄｆ北とし、前記Ｄｆ東と前記Ｄｆ西が異符号の場合はゼロ、同符号の場合はその絶対値の小さいものをＤｆ東西とし、前記Ｄｆ南と前記Ｄｆ北が異符号の場合はゼロ、同符号の場合はその絶対値の小さいものをＤｆ南北とし、前記Ｄｆ東西と前記Ｄｆ南北の傾きを有し前記胞甲の中心において胞甲の中心での前記関数ｆの値を持つ平面を構成してそれを平面乙とし、
前記胞甲の北東方向の格子点における前記平面乙の値が胞甲の中心の値より大きい場合は、前記平面乙の該北東方向の格子点での値が、該格子点を共有する胞甲以外の３胞の中心でのｆの値を超えないようにＤｆ東西及びＤｆ南北を１以下の最大の定数により定数倍することにより前記平面乙を変更する手段と、
胞甲の北東方向の格子点における平面乙の値が胞甲の中心の値より小さい場合は、前記平面乙の該北東方向の格子点での値が、該格子点を共有する胞甲以外の３胞の中心での前記関数ｆの値を下回らないように前記Ｄｆ東西及び前記Ｄｆ南北を１以下の最大の定数により定数倍することにより前記平面乙を変更する手段と、
得られた平面乙に対して、該北東方向の格子点に対して行った操作を、北西方向の格子点、南西方向の格子点及び南東方向の格子点に対してそれぞれ行うことにより、該平面乙の傾きを変更させ、得られた平面を該胞甲内での前記補間関数Ｆとして定める手段とを備えた事を特徴する実行装置。
【請求項８】
３次元実領域上の構成された３次元構造格子と、該３次元構造格子の各胞中心において値が定まる事によって定義される該３次元構造格子上の関数ｆとが用意されたとき、該関数ｆと該格子の幾何を用いることにより、前記３次元実領域上の関数Ｆを定義する補間方法を実行する実行装置において、
着目する胞を胞甲とし、
胞甲の各辺を東西、南北あるいは天地にそれぞれ合わせた際に、東西南北天地方向に隣合う６胞に対して、前記胞甲の中心点と各６胞の中心点との前記関数ｆの片側差分値の２倍の値をそれぞれＤｆ東、Ｄｆ西、Ｄｆ南、Ｄｆ北、Ｄｆ天、Ｄｆ地とし、前記Ｄｆ東と前記Ｄｆ西が異符号の場合はゼロ、同符号の場合はその絶対値の小さいものをＤｆ東西とし、前記Ｄｆ南と前記Ｄｆ北が異符号の場合はゼロ、同符号の場合はその絶対値の小さいものをＤｆ南北とし、前記Ｄｆ天と前記Ｄｆ地が異符号の場合はゼロ、同符号の場合はその絶対値の小さいものをＤｆ天地とし、前記Ｄｆ東西と前記Ｄｆ南北と前記Ｄｆ天地の傾きを有し前記胞甲の中心における関数ｆの値を通る補間関数の候補を構成してそれを平面乙とし、前記胞甲の北東天方向の格子点における平面乙の値が胞甲の中心の値より大きい場合は、平面乙の該北東天方向の格子点での値が、該格子点を共有する胞甲以外の７胞の中心でのｆ値を超えないように、前記Ｄｆ東西、前記Ｄｆ南北及び前記Ｄｆ天地を１以下の最大の定数により定数倍することにより平面乙を変更する手段と、
前記胞甲の北東天方向の格子点における平面乙の値が胞甲の中心の値より小さい場合は、平面乙の該北東天方向の格子点での値が、該格子点を共有する胞甲以外の７胞の中心でのｆ値を下回らないように前記Ｄｆ東西、前記Ｄｆ南北及び前記Ｄｆ天地を１以下の最大の定数により定数倍することにより平面乙を変更する手段と、
得られた補間関数の候補平面乙に対して、該北東天方向の格子点に対して行った操作を、北西天方向の格子点、南西天方向の格子点、南東天方向の格子点、北東地方向の格子点、北西地方向の格子点、南西地方向の格子点、及び南東地方向の格子点に対してそれぞれ行うことにより、平面乙の傾きを変更させ、得られた平面を該胞甲内での前記補間関数Ｆとして定める手段とを備えた事を特徴とする実行装置。
【請求項９】
１次元実領域上に構成された１次元構造格子と、該１次元構造格子の各胞中心において値が定まる事によって定義される該１次元構造格子上の関数ｆとが用意されたとき、該関数ｆと該格子の幾何を用いることにより、前記１次元実領域上の関数Ｆを定義する補間方法を実行するためのコンピュータで実行可能な制御プログラムであって、
着目する胞に対して、該格子上の点を前進差分と後進差分の符号が異なるときには傾きをゼロとし、該格子上の点を前進差分と後進差分の符号が同じ場合は、その絶対値の小さい方の２倍を傾きに選び、中心での値として前記着目する胞の中心の値ｆを持ち、該傾きを持つ１次関数によって、該胞で定まる１次元実領域の部分領域上での該関数Ｆの値を定めるような内容を有することを特徴とする制御プログラム。
【請求項１０】
２次元実領域上の構成された２次元構造格子と、該２次元構造格子の各胞中心において値が定まる事によって定義される該２次元構造格子上の関数ｆとが用意されたとき、該関数ｆと該格子の幾何を用いることにより、前記２次元実領域上の関数Ｆを定義する補間方法を実行するためのコンピュータで実行可能な制御プログラムであって、
着目する胞を胞甲とし、
胞甲の各辺を東西あるいは南北にそれぞれ合わせた際に東西南北方向に隣合う４胞に対して、前記胞甲の中心点と各４胞の中心点との前記関数ｆの片側差分値の２倍の値をそれぞれＤｆ東、Ｄｆ西、Ｄｆ南、Ｄｆ北とし、前記Ｄｆ東と前記Ｄｆ西が異符号の場合はゼロ、同符号の場合はその絶対値の小さいものをＤｆ東西とし、前記Ｄｆ南と前記Ｄｆ北が異符号の場合はゼロ、同符号の場合はその絶対値の小さいものをＤｆ南北とし、前記Ｄｆ東西と前記Ｄｆ南北の傾きを有し前記胞甲の中心において胞甲の中心での前記関数ｆの値を持つ平面を構成してそれを平面乙とし、
前記胞甲の北東方向の格子点における前記平面乙の値が胞甲の中心の値より大きい場合は、前記平面乙の該北東方向の格子点での値が、該格子点を共有する胞甲以外の３胞の中心でのｆの値を超えないようにＤｆ東西及びＤｆ南北を１以下の最大の定数により定数倍することにより前記平面乙を変更し、
胞甲の北東方向の格子点における平面乙の値が胞甲の中心の値より小さい場合は、前記平面乙の該北東方向の格子点での値が、該格子点を共有する胞甲以外の３胞の中心でのｆの値を下回らないように前記Ｄｆ東西及び前記Ｄｆ南北を１以下の最大の定数により定数倍することにより前記平面乙を変更し、
得られた平面乙に対して、該北東方向の格子点に対して行った操作を、北西方向の格子点、南西方向の格子点及び南東方向の格子点に対してそれぞれ行うことにより、該平面乙の傾きを変更させ、得られた平面を該胞甲内での前記補間関数Ｆとして定めるような内容を有する事を特徴する制御プログラム。
【請求項１１】
３次元実領域上の構成された３次元構造格子と、該３次元構造格子の各胞中心において値が定まる事によって定義される該３次元構造格子上の関数ｆとが用意されたとき、該関数ｆと該格子の幾何を用いることにより、前記３次元実領域上の関数Ｆを定義する補間方法において、
着目する胞を胞甲とし、
胞甲の各辺を東西、南北あるいは天地にそれぞれ合わせた際に、東西南北天地方向に隣合う６胞に対して、前記胞甲の中心点と各６胞の中心点との前記関数ｆの片側差分値の２倍の値をそれぞれＤｆ東、Ｄｆ西、Ｄｆ南、Ｄｆ北、Ｄｆ天、Ｄｆ地とし、前記Ｄｆ東と前記Ｄｆ西が異符号の場合はゼロ、同符号の場合はその絶対値の小さいものをＤｆ東西とし、前記Ｄｆ南と前記Ｄｆ北が異符号の場合はゼロ、同符号の場合はその絶対値の小さいものをＤｆ南北とし、前記Ｄｆ天と前記Ｄｆ地が異符号の場合はゼロ、同符号の場合はその絶対値の小さいものをＤｆ天地とし、前記Ｄｆ東西と前記Ｄｆ南北と前記Ｄｆ天地の傾きを有し前記胞甲の中心における関数ｆの値を通る補間関数の候補を構成してそれを平面乙とし、前記胞甲の北東天方向の格子点における平面乙の値が胞甲の中心の値より大きい場合は、平面乙の該北東天方向の格子点での値が、該格子点を共有する胞甲以外の７胞の中心でのｆ値を超えないように、前記Ｄｆ東西、前記Ｄｆ南北及び前記Ｄｆ天地を１以下の最大の定数により定数倍することにより平面乙を変更し、前記胞甲の北東天方向の格子点における平面乙の値が胞甲の中心の値より小さい場合は、平面乙の該北東天方向の格子点での値が、該格子点を共有する胞甲以外の７胞の中心でのｆ値を下回らないように前記Ｄｆ東西、前記Ｄｆ南北及び前記Ｄｆ天地を１以下の最大の定数により定数倍することにより平面乙を変更し、得られた補間関数の候補平面乙に対して、該北東天方向の格子点に対して行った操作を、北西天方向の格子点、南西天方向の格子点、南東天方向の格子点、北東地方向の格子点、北西地方向の格子点、南西地方向の格子点、及び南東地方向の格子点に対してそれぞれ行うことにより、平面乙の傾きを変更させ、得られた平面を該胞甲内での前記補間関数Ｆとして定めるような内容を有する事を特徴とする制御プログラム。
【発明の詳細な説明】
【０００１】
【発明の属する技術分野】
流体の運動を計算機上で実現する計算流体力学などで使用される補間方法等に関する。
【０００２】
【従来の技術】
流体の運動を計算機上で実現する計算流体力学は、近年大きく発展してきた。
しかしながら、流体中を物体が運動する様子を表現する、あるいは、２流体あるいは多層流体の運動を表現するということは困難を伴い、様々な試みがなされている。
【０００３】
静止した物体に対する流体の計算は、物体の形状に応じて格子を生成するということができるため、有限要素法等の非構造格子等を利用する方法が知られている。物体を動かす場合も、その動きに応じて格子を変化させるという方法もあるが、大きな変形を伴う場合、一般に困難を生じる。
【０００４】
他方、スタッガード格子として知られる双対格子と格子との双対性を利用した構造格子上での差分方程式による流体方程式の近似的実現は、その安定性から信頼性が高い。このことから、物質の占有率を格子上の関数として表現し、その関数の移動を行うという考え方がある。その１つとして例えば、非圧縮性流体における流体体等法（以下ＶＯＦ法と呼ぶ）を利用する方法がある。この場合は、占有率関数の移動問題が重要な技術となる。
【０００５】
また、圧縮性流体の場合は、様々な物の移動が重要であり、特に衝撃波の伝播等は、上記の物質、あるいは界面の運動と非常に関連している。
【０００６】
これらの移動を表現する方程式（移流方程式と呼ばれる）の初期値問題のことをリーマンの空気中の波動方程式の論文「Ｕｅｂｅｒ　ｄｉｅ　Ｆｏｒｔｐｆｌａｎｚｕｎｇ　ｅｂｅｎｅｒ　Ｌｕｆｔｗｅｌｌｅｎ　ｖｏｎ　ｅｎｄｌｉｃｈｅｒ　Ｓｃｈｗｉｎｇｕｎｇｓｗｅｉｔｅ　（Ａｂｈａｎｄｌｕｎｇｅｎ　ｄｅｒ　Ｋｏｎｉｇｌｉｃｈｅｎ　Ｇｅｓｅｌｌｓｃｈａｆｔ　ｄｅｒ　Ｗｉｓｓｅｎｓｃｈａｆｔｅｎ　ｚｕ　Ｇｏｔｔｉｎｇｅｎ，　８　（１８６０），　４３−６５）」に因んで、リーマン問題と呼ぶ。
【０００７】
このリーマン問題の数値解法として、近年全変動逓減法（Ｔｏｔａｌ　Ｖａｌｕａｔｉｏｎ　Ｄｉｍｉｎｉｓｈｉｎｇ法、以下ＴＶＤ法と記す）と呼ばれる高次精度上流化スキーム、制限関数を導入する本質的非振動法（Ｅｓｓｅｎｔｉａｌｌｙ　Ｎｏｎ−Ｏｓｃｉｌｌａｔｏｒｙ法、以下ＥＮＯ法と記す）、保存単調上流法（Ｍｏｎｅｔｏｎｅ　Ｕｐｓｔｒｅａｍ　Ｓｃｈｅｍｅ　ｆｏｒ　Ｃｏｎｓｅｒｖａｔｉｏｎ　Ｌａｗｓ、以下ＭＵＳＣＬ法と記す）、正値保存法（Ｐｏｓｉｔｉｖｅｌｙ　Ｐｒｅｓｅｒｖｉｎｇ，以下ＰＰ法と記す）等々の概念や方法が存在し、それらは独立であったり、従属であったりしている。
【０００８】
本発明で取り扱うものの一つは、
∂　Ｆ（ｘ，ｔ）　／　∂ｔ　＝∂Ｃ（ｘ，ｔ）　Ｆ（ｘ，ｔ）　／　∂ｘ…（１）
とする１次元の偏微分方程式であり、この方程式を差分化して計算を行う問題を考える。Ｆ（Ｘ，Ｔ）として、Ｘ，Ｔは以下離散的空間格子点及び離散時間を指すものとする。
【０００９】
式（１）を解く問題に対して、上述したように、計算流体力学においては様々な計算手法が提案されている。ＴＶＤ法とは、ＴＶＤ条件
ＴＶ（Ｆ［・，ｎ］））　＝　Σｊ　｜Ｆ［ｊ＋１，ｎ］　−　Ｆ［ｊ，ｎ］｜…（２）
が時間発展に対して増加しないという条件、
ＴＶ（Ｆ［・，ｎ＋１］）　≦　ＴＶ（Ｆ［・，ｎ］）…（３）
を満たす方法のことである。
【００１０】
広い意味での保存単調上流法（Ｍｏｎｅｔｏｎｅ　Ｕｐｓｔｒｅａｍ　Ｓｃｈｅｍｅ　ｆｏｒ　Ｃｏｎｓｅｒｖａｔｉｏｎ　Ｌａｗｓ　、ＭＵＳＣＬ法と記す）は、データの再構築と近似リーマン解法との２つから構成される計算手法である。
【００１１】
ＭＵＳＣＬ法とＴＶＤ法を統合することにより、関数に不連続な飛びが存在する場合に発生する振動を除去する制限関数付き上流法が、ＥＮＯ法である。
【００１２】
また、データの再構築方法としては、区分的定関数で補間する方法（これを１次精度と記す）、関数を区分的線形関数で補間する方法（これを２次精度と記す）、及び区分的２次曲線で補間する方法（これを３次精度と記す）が多く使われている。
【００１３】
本発明は、２次精度に関するものである。
【００１４】
従来の方法における２次精度のＥＮＯ法について説明を行う。上記式（１）を数値的に解く際に、区分的１次関数
Ｈ［Ｉ，Ｊ］（ｘ）　＝　Ｆ（Ｉ，Ｊ）　＋　Ｇ（Ｉ，Ｊ）ｘ…（４）
が存在しているとして差分近似解を求める方法を述べる。もしも関数Ｃ（Ｉ＋１／２，Ｊ）が正であれば、
Ｑ（Ｉ＋１／２，Ｊ）　：＝　Ｃ（Ｉ＋１／２，Ｊ）Ｈ［Ｉ，Ｊ］（Δｘ／２　−　Ｃ（Ｉ，Ｊ）　Δｔ／２）…（５）
とし、そうでなければ、
Ｑ（Ｉ＋１／２，Ｊ）　：＝　Ｃ（Ｉ＋１／２，Ｊ）Ｈ［Ｉ＋１，Ｊ］（−Δｘ／２　−　Ｃ（Ｉ＋１，Ｊ）　Δｔ／２）…（６）
と定義し、
Ｆ（Ｉ，Ｊ＋１）　＝　Ｆ（Ｉ，Ｊ）　＋　Δｔ　（Ｑ（Ｉ＋１／２，Ｊ）　−　Ｑ（Ｉ−１／２，Ｊ）　）／Δｘ…（７）
とすることにより、式（１）を近似的に解くことができる。
【００１５】
他方、式（４）の区分的１次関数において勾配を表すスロープ関数Ｇ（Ｉ，Ｊ）は以下のように定める。
【００１６】
Ｇ（Ｉ，Ｊ）＝ΔＦ（Ｉ，Ｊ）／Δｘ…（８）
ここで、ΔＦ（Ｉ，Ｊ）の求め方は以下に示すように、様々な定義の仕方がある。ここでは特徴的な３つの方法について述べる。
【００１７】
まず、
2004062749.tif 000003
を導入する。１次元ＴＶＤ制限関数として、Ｍｉｎｍｏｄと呼ばれるものにおいては、ΔＦ（Ｉ，Ｊ）は、
2004062749.tif 000004
で与えられるものであり、平均制限関数と呼ばれるものでは、
2004062749.tif 000005
であり、Ｓｕｐｅｒｂｅｅと呼ばれるものでは、
2004062749.tif 000006
であるとしている。但し、ｓｇｎ（）は符号関数である。つまり、その引数が正のときには１、負の時には−１を返す関数である。
【００１８】
これらのことは、例えば「Ａ．　Ｓｕｒｅｓｈ，　“Ｐｏｓｉｔｉｖｉｔｙ−Ｐｅｒｓｅｒｖｉｎｇ　Ｓｃｈｅｍｅｓ　ｉｎ　Ｍｕｌｔｉｄｉｍｅｎｓｉｏｎｓ”ＳＩＡＭ　Ｊ．　Ｓｃｉ．　Ｃｏｍｐｕｔ．　２２，　１１８４−１１９８　（２０００）」に記載されている。
【００１９】
これらの補間方法で、矩形形状の関数を上記式（７）に従って移流計算した結果が図７〜図９に示されている。周期的境界条件で、ステップ数１００回の後の形状と、初期形状をグラフに示した。図７がＭｉｎｍｏｄ、図８が平均制限関数、図９がＳｕｐｅｒｂｅｅと呼ばれるものである。これらの図から明らかなように、矩形形状が激しく崩れ鈍っていることが判る。
【００２０】
次に、上述のような１次元に対するアルゴリズムを高次元に拡張したものについて説明する。ここでは、上記文献に記載されている方法について述べる。
【００２１】
二次元格子上に与えられたある時間Ｔでの関数をＦ（Ｉ，Ｊ，Ｔ）とする。Ｉ，Ｊは整数であり、２次元の格子の各格子を指している。
【００２２】
この時、次の値を定義する。
2004062749.tif 000007
また
2004062749.tif 000008
とする。
【００２３】
2004062749.tif 000009
この時、
2004062749.tif 000010
とした後に、（Ｉ，Ｊ）格子内の補間関数Ｈ（Ｉ，Ｊ，Ｔ）（ｘ，ｙ）を
2004062749.tif 000011
とするものである。
【００２４】
これを利用すると、図１０（ａ）に示す初期値に対して、図１０（ｂ）に示すような移流の計算ができることになる。
【００２５】
【発明が解決しようとする課題】
しかしながら上記従来の補間方法では、次のような問題点があった。
【００２６】
２相の非圧縮性流体における流体体積法（ＶＯＦ法）において、形状を如何に保持させたまま時間発展させるかという問題は、非常に大きな関心事であるにも関わらず、従来の方法では、計算ステップを経る毎に形状が崩れてなだらかになってゆく傾向があった。
【００２７】
このような点から、より形状を保持できるようなような計算アルゴリズムが必須であり、そのために、一次元に対するアルゴリズムについては、上記式（４）における補間関数をより良いものに取り替えることが必要であり、高次元に対するアルゴリズムについては、上記式（１７）における補間関数をより良いものに取り替えることが必要であった。
【００２８】
本発明は上記従来の問題点に鑑み、２相の非圧縮性流体における流体体積法（ＶＯＦ法）において、形状をシャープに保持させたまま時間発展させることができる補間方法等を提供することを目的とする。
【００２９】
【課題を解決するための手段】
上記目的を達成するために、本発明の第一の補間方法では、１次元実領域上に構成された１次元構造格子と、該１次元構造格子の各胞中心において値が定まる事によって定義される該１次元構造格子上の関数ｆとが用意されたとき、該関数ｆと該格子の幾何を用いることにより、前記１次元実領域上の関数Ｆを定義する補間方法において、着目する胞に対して、該格子上の点を前進差分と後進差分の符号が異なるときには傾きをゼロとし、該格子上の点を前進差分と後進差分の符号が同じ場合は、その絶対値の小さい方の２倍を傾きに選び、中心での値として前記着目する胞の中心の値ｆを持ち、該傾きを持つ１次関数によって、該胞で定まる１次元実領域の部分領域上での該関数Ｆの値を定めることを特徴とする。
【００３０】
本発明の第二の補間方法では、２次元実領域上の構成された２次元構造格子と、該２次元構造格子の各胞中心において値が定まる事によって定義される該２次元構造格子上の関数ｆとが用意されたとき、該関数ｆと該格子の幾何を用いることにより、前記２次元実領域上の関数Ｆを定義する補間方法において、着目する胞を胞甲とし、胞甲の各辺を東西あるいは南北にそれぞれ合わせた際に東西南北方向に隣合う４胞に対して、前記胞甲の中心点と各４胞の中心点との前記関数ｆの片側差分値の２倍の値をそれぞれＤｆ東、Ｄｆ西、Ｄｆ南、Ｄｆ北とし、前記Ｄｆ東と前記Ｄｆ西が異符号の場合はゼロ、同符号の場合はその絶対値の小さいものをＤｆ東西とし、前記Ｄｆ南と前記Ｄｆ北が異符号の場合はゼロ、同符号の場合はその絶対値の小さいものをＤｆ南北とし、前記Ｄｆ東西と前記Ｄｆ南北の傾きを有し前記胞甲の中心において胞甲の中心での前記関数ｆの値を持つ平面を構成してそれを平面乙とし、前記胞甲の北東方向の格子点における前記平面乙の値が胞甲の中心の値より大きい場合は、前記平面乙の該北東方向の格子点での値が、該格子点を共有する胞甲以外の３胞の中心でのｆの値を超えないようにＤｆ東西及びＤｆ南北を１以下の最大の定数により定数倍することにより前記平面乙を変更し、胞甲の北東方向の格子点における平面乙の値が胞甲の中心の値より小さい場合は、前記平面乙の該北東方向の格子点での値が、該格子点を共有する胞甲以外の３胞の中心でのｆの値を下回らないように前記Ｄｆ東西及び前記Ｄｆ南北を１以下の最大の定数により定数倍することにより前記平面乙を変更し、得られた平面乙に対して、該北東方向の格子点に対して行った操作を、北西方向の格子点、南西方向の格子点及び南東方向の格子点に対してそれぞれ行うことにより、該平面乙の傾きを変更させ、得られた平面を該胞甲内での前記補間関数Ｆとして定める事を特徴とする。
【００３１】
本発明の第三の補間方法では、３次元実領域上の構成された３次元構造格子と、該３次元構造格子の各胞中心において値が定まる事によって定義される該３次元構造格子上の関数ｆとが用意されたとき、該関数ｆと該格子の幾何を用いることにより、前記３次元実領域上の関数Ｆを定義する補間方法において、着目する胞を胞甲とし、胞甲の各辺を東西、南北あるいは天地にそれぞれ合わせた際に、東西南北天地方向に隣合う６胞に対して、前記胞甲の中心点と各６胞の中心点との前記関数ｆの片側差分値の２倍の値をそれぞれＤｆ東、Ｄｆ西、Ｄｆ南、Ｄｆ北、Ｄｆ天、Ｄｆ地とし、前記Ｄｆ東と前記Ｄｆ西が異符号の場合はゼロ、同符号の場合はその絶対値の小さいものをＤｆ東西とし、前記Ｄｆ南と前記Ｄｆ北が異符号の場合はゼロ、同符号の場合はその絶対値の小さいものをＤｆ南北とし、前記Ｄｆ天と前記Ｄｆ地が異符号の場合はゼロ、同符号の場合はその絶対値の小さいものをＤｆ天地とし、前記Ｄｆ東西と前記Ｄｆ南北と前記Ｄｆ天地の傾きを有し前記胞甲の中心における関数ｆの値を通る補間関数の候補を構成してそれを平面乙とし、前記胞甲の北東天方向の格子点における平面乙の値が胞甲の中心の値より大きい場合は、平面乙の該北東天方向の格子点での値が、該格子点を共有する胞甲以外の７胞の中心でのｆ値を超えないように、前記Ｄｆ東西、前記Ｄｆ南北及び前記Ｄｆ天地を１以下の最大の定数により定数倍することにより平面乙を変更し、前記胞甲の北東天方向の格子点における平面乙の値が胞甲の中心の値より小さい場合は、平面乙の該北東天方向の格子点での値が、該格子点を共有する胞甲以外の７胞の中心でのｆ値を下回らないように前記Ｄｆ東西、前記Ｄｆ南北及び前記Ｄｆ天地を１以下の最大の定数により定数倍することにより平面乙を変更し、得られた補間関数の候補平面乙に対して、該北東天方向の格子点に対して行った操作を、北西天方向の格子点、南西天方向の格子点、南東天方向の格子点、北東地方向の格子点、北西地方向の格子点、南西地方向の格子点、及び南東地方向の格子点に対してそれぞれ行うことにより、平面乙の傾きを変更させ、得られた平面を該胞甲内での前記補間関数Ｆとして定める事を特徴とする。
【００３２】
本発明の第一の実行装置では、１次元実領域上に構成された１次元構造格子と、該１次元構造格子の各胞中心において値が定まる事によって定義される該１次元構造格子上の関数ｆとが用意されたとき、該関数ｆと該格子の幾何を用いることにより、前記１次元実領域上の関数Ｆを定義する補間方法を実行する実行装置において、着目する胞に対して、該格子上の点を前進差分と後進差分の符号が異なるときには傾きをゼロに設定し、該格子上の点を前進差分と後進差分の符号が同じ場合は、その絶対値の小さい方の２倍を傾きに設定する手段と、中心での値として前記着目する胞の中心の値ｆを持ち、該傾きを持つ１次関数によって、該胞で定まる１次元実領域の部分領域上での該関数Ｆの値を定める手段とを備えたことを特徴とする。
【００３３】
本発明の第二の実行装置では、２次元実領域上の構成された２次元構造格子と、該２次元構造格子の各胞中心において値が定まる事によって定義される該２次元構造格子上の関数ｆとが用意されたとき、該関数ｆと該格子の幾何を用いることにより、前記２次元実領域上の関数Ｆを定義する補間方法を実行する実行装置において、着目する胞を胞甲とし、胞甲の各辺を東西あるいは南北にそれぞれ合わせた際に東西南北方向に隣合う４胞に対して、前記胞甲の中心点と各４胞の中心点との前記関数ｆの片側差分値の２倍の値をそれぞれＤｆ東、Ｄｆ西、Ｄｆ南、Ｄｆ北とし、前記Ｄｆ東と前記Ｄｆ西が異符号の場合はゼロ、同符号の場合はその絶対値の小さいものをＤｆ東西とし、前記Ｄｆ南と前記Ｄｆ北が異符号の場合はゼロ、同符号の場合はその絶対値の小さいものをＤｆ南北とし、前記Ｄｆ東西と前記Ｄｆ南北の傾きを有し前記胞甲の中心において胞甲の中心での前記関数ｆの値を持つ平面を構成してそれを平面乙とし、前記胞甲の北東方向の格子点における前記平面乙の値が胞甲の中心の値より大きい場合は、前記平面乙の該北東方向の格子点での値が、該格子点を共有する胞甲以外の３胞の中心でのｆの値を超えないようにＤｆ東西及びＤｆ南北を１以下の最大の定数により定数倍することにより前記平面乙を変更する手段と、胞甲の北東方向の格子点における平面乙の値が胞甲の中心の値より小さい場合は、前記平面乙の該北東方向の格子点での値が、該格子点を共有する胞甲以外の３胞の中心での前記関数ｆの値を下回らないように前記Ｄｆ東西及び前記Ｄｆ南北を１以下の最大の定数により定数倍することにより前記平面乙を変更する手段と、得られた平面乙に対して、該北東方向の格子点に対して行った操作を、北西方向の格子点、南西方向の格子点及び南東方向の格子点に対してそれぞれ行うことにより、該平面乙の傾きを変更させ、得られた平面を該胞甲内での前記補間関数Ｆとして定める手段とを備えた事を特徴とする。
【００３４】
本発明の第三の実行装置では、３次元実領域上の構成された３次元構造格子と、該３次元構造格子の各胞中心において値が定まる事によって定義される該３次元構造格子上の関数ｆとが用意されたとき、該関数ｆと該格子の幾何を用いることにより、前記３次元実領域上の関数Ｆを定義する補間方法を実行する実行装置において、着目する胞を胞甲とし、胞甲の各辺を東西、南北あるいは天地にそれぞれ合わせた際に、東西南北天地方向に隣合う６胞に対して、前記胞甲の中心点と各６胞の中心点との前記関数ｆの片側差分値の２倍の値をそれぞれＤｆ東、Ｄｆ西、Ｄｆ南、Ｄｆ北、Ｄｆ天、Ｄｆ地とし、前記Ｄｆ東と前記Ｄｆ西が異符号の場合はゼロ、同符号の場合はその絶対値の小さいものをＤｆ東西とし、前記Ｄｆ南と前記Ｄｆ北が異符号の場合はゼロ、同符号の場合はその絶対値の小さいものをＤｆ南北とし、前記Ｄｆ天と前記Ｄｆ地が異符号の場合はゼロ、同符号の場合はその絶対値の小さいものをＤｆ天地とし、前記Ｄｆ東西と前記Ｄｆ南北と前記Ｄｆ天地の傾きを有し前記胞甲の中心における関数ｆの値を通る補間関数の候補を構成してそれを平面乙とし、前記胞甲の北東天方向の格子点における平面乙の値が胞甲の中心の値より大きい場合は、平面乙の該北東天方向の格子点での値が、該格子点を共有する胞甲以外の７胞の中心でのｆ値を超えないように、前記Ｄｆ東西、前記Ｄｆ南北及び前記Ｄｆ天地を１以下の最大の定数により定数倍することにより平面乙を変更する手段と、前記胞甲の北東天方向の格子点における平面乙の値が胞甲の中心の値より小さい場合は、平面乙の該北東天方向の格子点での値が、該格子点を共有する胞甲以外の７胞の中心でのｆ値を下回らないように前記Ｄｆ東西、前記Ｄｆ南北及び前記Ｄｆ天地を１以下の最大の定数により定数倍することにより平面乙を変更する手段と、得られた補間関数の候補平面乙に対して、該北東天方向の格子点に対して行った操作を、北西天方向の格子点、南西天方向の格子点、南東天方向の格子点、北東地方向の格子点、北西地方向の格子点、南西地方向の格子点、及び南東地方向の格子点に対してそれぞれ行うことにより、平面乙の傾きを変更させ、得られた平面を該胞甲内での前記補間関数Ｆとして定める手段とを備えた事を特徴とする。
【００３５】
本発明の第一の制御プログラムでは、１次元実領域上に構成された１次元構造格子と、該１次元構造格子の各胞中心において値が定まる事によって定義される該１次元構造格子上の関数ｆとが用意されたとき、該関数ｆと該格子の幾何を用いることにより、前記１次元実領域上の関数Ｆを定義する補間方法を実行するためのコンピュータで実行可能な制御プログラムであって、着目する胞に対して、該格子上の点を前進差分と後進差分の符号が異なるときには傾きをゼロとし、該格子上の点を前進差分と後進差分の符号が同じ場合は、その絶対値の小さい方の２倍を傾きに選び、中心での値として前記着目する胞の中心の値ｆを持ち、該傾きを持つ１次関数によって、該胞で定まる１次元実領域の部分領域上での該関数Ｆの値を定めるような内容を有することを特徴とする。
【００３６】
本発明の第二の制御プログラムでは、２次元実領域上の構成された２次元構造格子と、該２次元構造格子の各胞中心において値が定まる事によって定義される該２次元構造格子上の関数ｆとが用意されたとき、該関数ｆと該格子の幾何を用いることにより、前記２次元実領域上の関数Ｆを定義する補間方法を実行するためのコンピュータで実行可能な制御プログラムであって、着目する格子を胞甲とし、胞甲の各辺を東西あるいは南北にそれぞれ合わせた際に東西南北方向に隣合う４胞に対して、前記胞甲の中心点と各４胞の中心点との前記関数ｆの片側差分値の２倍の値をそれぞれＤｆ東、Ｄｆ西、Ｄｆ南、Ｄｆ北とし、前記Ｄｆ東と前記Ｄｆ西が異符号の場合はゼロ、同符号の場合はその絶対値の小さいものをＤｆ東西とし、前記Ｄｆ南と前記Ｄｆ北が異符号の場合はゼロ、同符号の場合はその絶対値の小さいものをＤｆ南北とし、前記Ｄｆ東西と前記Ｄｆ南北の傾きを有し前記胞甲の中心において胞甲の中心での前記関数ｆの値を持つ平面を構成してそれを平面乙とし、前記胞甲の北東方向の格子点における前記平面乙の値が胞甲の中心の値より大きい場合は、前記平面乙の該北東方向の格子点での値が、該格子点を共有する胞甲以外の３胞の中心でのｆの値を超えないようにＤｆ東西及びＤｆ南北を１以下の最大の定数により定数倍することにより前記平面乙を変更し、胞甲の北東方向の格子点における平面乙の値が胞甲の中心の値より小さい場合は、前記平面乙の該北東方向の格子点での値が、該格子点を共有する胞甲以外の３胞の中心でのｆの値を下回らないように前記Ｄｆ東西及び前記Ｄｆ南北を１以下の最大の定数により定数倍することにより前記平面乙を変更し、得られた平面乙に対して、該北東方向の格子点に対して行った操作を、北西方向の格子点、南西方向の格子点及び南東方向の格子点に対してそれぞれ行うことにより、該平面乙の傾きを変更させ、得られた平面を該胞甲内での前記補間関数Ｆとして定めるような内容を有する事を特徴する。
【００３７】
本発明の第三の制御プログラムでは、３次元実領域上の構成された３次元構造格子と、該３次元構造格子の各胞中心において値が定まる事によって定義される該３次元構造格子上の関数ｆとが用意されたとき、該関数ｆと該格子の幾何を用いることにより、前記３次元実領域上の関数Ｆを定義する補間方法において、着目する格子を胞甲とし、胞甲の各辺を東西、南北あるいは天地にそれぞれ合わせた際に、東西南北天地方向に隣合う６胞に対して、前記胞甲の中心点と各６胞の中心点との前記関数ｆの片側差分値の２倍の値をそれぞれＤｆ東、Ｄｆ西、Ｄｆ南、Ｄｆ北、Ｄｆ天、Ｄｆ地とし、前記Ｄｆ東と前記Ｄｆ西が異符号の場合はゼロ、同符号の場合はその絶対値の小さいものをＤｆ東西とし、前記Ｄｆ南と前記Ｄｆ北が異符号の場合はゼロ、同符号の場合はその絶対値の小さいものをＤｆ南北とし、前記Ｄｆ天と前記Ｄｆ地が異符号の場合はゼロ、同符号の場合はその絶対値の小さいものをＤｆ天地とし、前記Ｄｆ東西と前記Ｄｆ南北と前記Ｄｆ天地の傾きを有し前記胞甲の中心における関数ｆの値を通る補間関数の候補を構成してそれを平面乙とし、前記胞甲の北東天方向の格子点における平面乙の値が胞甲の中心の値より大きい場合は、平面乙の該北東天方向の格子点での値が、該格子点を共有する胞甲以外の７胞の中心でのｆ値を超えないように、前記Ｄｆ東西、前記Ｄｆ南北及び前記Ｄｆ天地を１以下の最大の定数により定数倍することにより平面乙を変更し、前記胞甲の北東天方向の格子点における平面乙の値が胞甲の中心の値より小さい場合は、平面乙の該北東天方向の格子点での値が、該格子点を共有する胞甲以外の７胞の中心でのｆ値を下回らないように前記Ｄｆ東西、前記Ｄｆ南北及び前記Ｄｆ天地を１以下の最大の定数により定数倍することにより平面乙を変更し、得られた補間関数の候補平面乙に対して、該北東天方向の格子点に対して行った操作を、北西天方向の格子点、南西天方向の格子点、南東天方向の格子点、北東地方向の格子点、北西地方向の格子点、南西地方向の格子点、及び南東地方向の格子点に対してそれぞれ行うことにより、平面乙の傾きを変更させ、得られた平面を該胞甲内での前記補間関数Ｆとして定めるような内容を有する事を特徴とする。
【００３８】
【発明の実施の形態】
以下、本発明の実施の形態を図面に基づいて説明する。
【００３９】
［第１実施形態］
＜第１実施形態の補間方法＞
始めに本実施形態の補間方法について説明する。本実施形態は、データの再構築方法として、関数を区分的線形関数で補間する方法である。つまり、移流方程式の数値解法としては２次精度のものに対応するものである。
【００４０】
本実施形態では、１次元の問題として、
∂　Ｆ（ｘ，ｔ）　／　∂ｔ　＝∂Ｃ（ｘ，ｔ）　Ｆ（ｘ，ｔ）　／　∂ｘ…（１８）
とする偏微分方程式を差分化して計算を行う問題を考える。Ｆ（ｘ，ｔ）を離散化したものとしてＦ（Ｘ，Ｔ）を考える。Ｆの引数Ｘ，Ｔは、以下離散的な空間格子及び離散時間を意味するものとする。
【００４１】
まず、本実施形態の補間関数が、区分的１次関数
Ｈ［Ｉ，Ｊ］（ｘ）　＝　Ｆ（Ｉ，Ｊ）　＋　Ｇ（Ｉ，Ｊ）ｘ…（１９）
として、既に得られているとして、上記の偏微分方程式の差分近似解を求める方法を述べる。本発明の補間関数Ｈ［Ｉ，Ｊ］（ｘ）の係数Ｇ（Ｉ，Ｊ）の求め方に関しては後で述べる。そこで述べるように、Ｇ（Ｉ，Ｊ）の決定方法が本発明の補間アルゴリズムの特徴である。Ｇ（Ｉ，Ｊ）の決定に際し、格子の幾何学的性質である隣接情報を利用することによって決定を行っていることは重要である。
【００４２】
要するに、１次元実領域上の構成された１次元構造格子と、該１次元構造格子の各胞中心において値が定まる事によって定義される該１次元構造格子上の関数Ｆ（Ｘ，Ｔ）とが用意されたとき、該関数Ｆ（Ｘ，Ｔ）と該格子の幾何を利用することによって該関数Ｆ（Ｘ，Ｔ）を、該１次元実領域上のＨ［Ｉ，Ｊ］（ｘ）として定義する事になる。すなわち、離散的な点でのみ定義されているＦという写像が、Ｈではその定義域がｘ軸上全体の関数に広がる。
【００４３】
もしも関数Ｃ（Ｉ＋１／２，Ｊ）が正であれば、
Ｑ（Ｉ＋１／２，Ｊ）　：＝　Ｃ（Ｉ＋１／２，Ｊ）Ｈ［Ｉ，Ｊ］（Δｘ／２　−　Ｃ（Ｉ，Ｊ）　Δｔ／２）…（２０）
とし、そうでなければ、
Ｑ（Ｉ＋１／２，Ｊ）　：＝　Ｃ（Ｉ＋１／２，Ｊ）Ｈ［Ｉ＋１，Ｊ］（−Δｘ／２　−　Ｃ（Ｉ，Ｊ）　Δｔ／２）…（２１）
と定義し、
Ｆ（Ｉ，Ｊ＋１）　＝　Ｆ（Ｉ，Ｊ）　＋　Δｔ　（Ｑ（Ｉ＋１／２，Ｊ）　−　Ｑ（Ｉ−１／２，Ｊ）　）／Δｘ…（２２）
とする。
【００４４】
他方、区分的１次関数において勾配を表すスロープ関数Ｇ（Ｉ，Ｊ）を以下のように定める。まず、
2004062749.tif 000012
を導入する。本実施形態の１次元制限関数として、
2004062749.tif 000013
を導入する。但し、ｓｇｎ（）は符号関数である。つまり、その引数が正のときには１、負の時には−１を返す関数である。従来例のｍｉｎｍｏｄ制限関数より２だけ大きいことが重要である。このΔＦ　（Ｉ，Ｊ）を利用して
2004062749.tif 000014
と定義すればよい。
【００４５】
Ｓｐ（Ｉ，Ｊ）／Δｘは、上記式（２３）の定義から前進差分、Ｓｍ（Ｉ，Ｊ）／Δｘは上記式（２３）の定義から後進差分を意味する。それらが異符号のときは、（ｓｇｎ（Ｓｐ（Ｉ，Ｊ）　＋　ｓｇｎ（Ｓｍ（Ｉ，Ｊ））がゼロとなり、全体がゼロとなる。そうでないときは、Ｇ（Ｉ，Ｊ）は、（ｓｇｎ（Ｓｐ（Ｉ，Ｊ）　＋　ｓｇｎ（Ｓｍ（Ｉ，Ｊ））が２または−２になるから、後進差分あるいは前進差分の絶対値の小さいものの２倍と一致する。
【００４６】
このように本実施形態の補間方法は、１次元構造格子上の着目する胞に対して、該胞でのＦの前進差分と後進差分の符号が異なるときには傾きをゼロ、該格子上の点を前進差分と後進差分の符号が同じ場合は、その絶対値の小さい方の２倍を傾きに選び、中心での値を該着目する格子の中心での値Ｆ（Ｘ，Ｔ）を与え、１次関数によって、該胞で定まる１次元実領域の部分領域上での関数Ｈ（Ｉ，Ｊ）の値を定めるものである。
【００４７】
＜本実施形態の補間方法を利用した移流計算の流れ＞
図１は、上記した第１実施形態の補間方法を利用した移流計算の流れ図である。
【００４８】
まず、ステップＳ１で初期設定をして、例えば、格子の生成、時間間隔を決定した後に関数の初期形状Ｆ（・，１）を初期化する。
【００４９】
次のステップＳ２では、上述のＱ（・，Ｊ）の計算を行う（上記式（２０）及び式（２１）））。この時、上記式（２０）及び式（２１）に示すように、本実施形態の補間関数Ｈを利用する。そして、ステップＳ３でＱ（・，Ｊ）のデータを基に、Ｆ（・，Ｊ＋１）の計算を行う（上記式（２２））。
【００５０】
その後はステップＳ４で時刻を更新する。すなわち、今解こうとしている問題は上記式（１８）の時間発展方程式であるが、通常の時間発展方程式の数値解法の慣例に従い、時刻Ｔ＝ΔｔＪでのＦの値Ｆ（・，Ｊ）を、Ｆ（・，Ｊ−１）から求めた後に、時間をＴ＝Δｔ＊（Ｊ−１）からＴ＝Δｔ＊Ｊに変更を行う。
【００５１】
そして、ステップＳ５で最終ステップか否かの判断を行い、最終でなければ、ステップＳ２に戻り、最終ステップになるまで計算を行う。
【００５２】
本実施形態の補間方法を利用すると、関数のスロープがゼロのものとシャープなものの２つの領域に分離され、２相の非圧縮性流体の形状表現関数としては非常に良い結果を得ることができる。その結果を表示したのが図２である。従来例の計算結果である図７〜図９に対応した同じ条件で計算を行っており、初期値と１０００ステップ後の形状を図示している。図２で示す本実施形態の移流計算結果は、図７〜図９に示したものに比べて非常にシャープな矩形形状を保っていることが分かる。
【００５３】
＜本実施形態のハードウェア＞
図３は、本実施形態の特徴を成す補間方法を利用した移流計算を実行する計算機の概観図であり、図４はその要部の電気的構成を示すブロック図である。
【００５４】
図３及び図４に示すように、この計算機は、複数のＣＰＵ１１６やメモリ１１８等からなる演算処理部、ハードディスク１１７やフロッピー（登録商標）ディスク１１０からなる記憶装置部、キーボード１１５やマウス１１４等からなる入力部、そしてディスプレイ１１２等からなる出力表示部などで構成されているシステムである。各演算部等は、本体１１１に収まっている。
【００５５】
この計算機は、メモリ１１８上に図１に示したフローチャートに沿って設計しコード化されたプログラムをロードし、必要な計算領域を確保し、適当な方法で入力された領域形状や格子点座標および物理量に基づき所定の演算処理を行い、その結果得られた領域上での物理量の値をハードディスク１１７等に書き込んで記憶させたり、ディスプレイ１１２上に表示させる。
【００５６】
なお、本実施形態の特徴を成す上記補間方法を利用すると、不連続点を含みながらも、格子点上に与えられた値から実軸上の関数を構成することができる。このことより、本発明の補間方法を上記実施形態とは異なる状況、例えば信号処理のデータ抽出や、１次元画像処理等々で使用しても、同様にシャープな形状を表現するものとなる。
【００５７】
［第２実施形態］
第２実施形態では、上記第１実施形態で説明した一次元の格子に対する補間アルゴリズムを二次元に拡張したものである。
【００５８】
本実施形態で取り扱うものは、
2004062749.tif 000015
とする２次元空間での偏微分方程式であり、これをを差分化して計算を行う問題を考える。Ｆ（ｘ，ｙ，ｔ）を離散化したものとしてＦ（Ｘ，Ｙ，Ｔ）を考える。Ｆの引数Ｘ，Ｙ，Ｔは、以下離散的な空間格子及び離散時間を意味するものとする。
【００５９】
まず、本実施形態の補間関数が、区分的１次関数
Ｈ［Ｉ，Ｊ，Ｔ］（ｘ，ｙ）　＝　Ｆ（Ｉ，Ｊ，Ｔ）　＋　Ｇｘ（Ｉ，Ｊ，Ｔ）ｘ　＋　Ｇｙ（Ｉ，Ｊ，Ｔ）ｙ…（２７）
として、既に得られているとして、上記の偏微分方程式の差分近似解を求める方法を述べる。本実施形態の補間関数Ｈ［Ｉ，Ｊ，Ｔ］（ｘ，ｙ）の係数Ｇｘ（Ｉ，Ｊ，Ｔ）及びＧｙ（Ｉ，Ｊ，Ｔ）の求め方に関しては後で述べる。そこで述べるように、Ｇｘ（Ｉ，Ｊ，Ｔ）及びＧｙ（Ｉ，Ｊ，Ｔ）の決定方法が本実施形態の補間アルゴリズムの特徴である。Ｇｘ（Ｉ，Ｊ，Ｔ）及びＧｙ（Ｉ，Ｊ，Ｔ）の決定に際し、格子の幾何学的性質である隣接情報を利用することによって決定を行っていることは重要である。
【００６０】
要するに、２次元実領域上の構成された２次元構造格子と、該２次元構造格子の各胞中心において値が定まる事によって定義されるある時刻Ｔでの該２次元構造格子上の関数Ｆ（Ｘ，Ｙ，Ｔ）とが用意されたとき、該関数Ｆ（Ｘ，Ｙ，Ｔ）と該格子の幾何を利用することによって該関数Ｆ（Ｘ，Ｙ，Ｔ）を、該２次元実領域上のＨ［Ｉ，Ｊ，Ｔ］（ｘ，ｙ）として定義する事になる。
【００６１】
もしも　関数Ｃｘ（Ｉ＋１／２，Ｊ，Ｔ）が正であれば、
2004062749.tif 000016
とし、そうでなければ、
2004062749.tif 000017
と定義し、
また、もしも　関数Ｃｙ（Ｉ，Ｊ＋１／２，Ｔ）が正であれば、
2004062749.tif 000018
とし、そうでなければ、
2004062749.tif 000019
と定義し、
2004062749.tif 000020
とする。
【００６２】
他方、区分的１次関数において勾配を表すスロープ関数Ｇ（Ｉ，Ｊ，Ｔ）を以下のように定める。
【００６３】
まず、
2004062749.tif 000021
を導入する。本実施形態の１次元制限関数として、
2004062749.tif 000022
を導入する。
【００６４】
次に、
2004062749.tif 000023
と、
2004062749.tif 000024
を定義する。
【００６５】
更に
ｈｐｐ（Ｉ，Ｊ，Ｔ）　＝　　　Ｓｘ（Ｉ，　Ｊ，　Ｔ）／２．　＋　　Ｓｘ（Ｉ，　Ｊ，　Ｔ）／２．
ｈｍｐ（Ｉ，Ｊ，Ｔ）　＝　−　Ｓｘ（Ｉ，　Ｊ，　Ｔ）／２．　＋　　Ｓｘ（Ｉ，　Ｊ，　Ｔ）／２．
ｈｐｍ（Ｉ，Ｊ，Ｔ）　＝　　　Ｓｘ（Ｉ，　Ｊ，　Ｔ）／２．　−　　Ｓｘ（Ｉ，　Ｊ，　Ｔ）／２．
ｈｍｍ（Ｉ，Ｊ，Ｔ）　＝　−　Ｓｘ（Ｉ，　Ｊ，　Ｔ）／２．　−　　Ｓｘ（Ｉ，　Ｊ，　Ｔ）／２．…（３７）
とする。
【００６６】
更に、「ｐｈｉｐｐ＝ｐｈｉｍｐ＝ｐｈｉｐｍ＝ｐｈｉｍｍ＝１としておき、もし（ｈｐｐ　＞　Ｖｐｐｍａｘ　）ならば、　ｐｈｉｐｐ　＝　Ｖｐｐｍａｘ／ｈｐｐと置き直す。もし、（ｈｐｐ　＜　Ｖｐｐｍｉｎ　）ならば、　ｐｈｉｐｐ　＝　Ｖｐｐｍｉｎ／ｈｐｐと置き直す。
【００６７】
同様に、もし（ｈｐｍ　＞　Ｖｐｍｍａｘ　）ならば、　ｐｈｉｐｍ　＝　Ｖｐｍｍａｘ／ｈｐｍと置き直す。
またもし、（ｈｐｍ　＜　Ｖｐｍｍｉｎ　）ならば、　ｐｈｉｐｍ　＝　Ｖｐｍｍｉｎ／ｈｐｍと置き直す。もし（ｈｍｐ　＞　Ｖｍｐｍａｘ　）ならば、　ｐｈｉｍｐ　＝　Ｖｍｐｍａｘ／ｈｍｐと置き直す。
またもし、（ｈｍｐ　＜　Ｖｍｐｍｉｎ　）ならば、　ｐｈｉｍｐ　＝　Ｖｍｐｍｉｎ／ｈｍｐと置き直す。もし（ｈｍｍ　＞　Ｖｍｍｍａｘ　）ならば、　ｐｈｉｍｍ　＝　Ｖｍｍｍａｘ／ｈｍｍと置き直す。またもし、（ｈｍｍ　＜　Ｖｍｍｍｉｎ　）ならば、　ｐｈｉｍｍ　＝　Ｖｍｍｍｉｎ／ｈｍｍと置き直す。」という計算を行う。
【００６８】
このようにして得られたｐｈｉｐｐ，ｐｈｉｍｐ，ｐｈｉｐｍ，ｐｈｉｍｍに対して、
Ｇｘ（Ｉ，Ｊ，　Ｔ）　＝　Ｍｉｎ（ｐｈｉｐｐ，　ｐｈｉｐｍ，　ｐｈｉｍｐ，　ｐｈｉｍｍ　）　Ｓｘ（Ｉ，Ｊ，Ｔ）／Δｘ
Ｇｙ（Ｉ，Ｊ，　Ｔ）　＝　Ｍｉｎ（ｐｈｉｐｐ，　ｐｈｉｐｍ，　ｐｈｉｍｐ，　ｐｈｉｍｍ　）　Ｓｙ（Ｉ，Ｊ，Ｔ）／Δｙ…（３８）
とした後に、（Ｉ，Ｊ）格子内の補間関数Ｈ（Ｉ，Ｊ，Ｔ）（ｘ，ｙ）を
Ｈ（Ｉ，Ｊ，Ｔ）（ｘ，ｙ）　＝　Ｆ（Ｉ，Ｊ，　Ｔ）　＋　Ｇｘ（Ｉ，Ｊ，　Ｔ）　ｘ　＋　Ｇｙ（Ｉ，Ｊ，　Ｔ）ｙ
…（３９）
とするものである。
【００６９】
このように本実施形態の補間方法を以下にまとめると、着目する格子を胞甲とし、胞甲の各辺を東西あるいは南北にそれぞれ合わせた際に東西南北方向に隣合う４胞に対して、胞甲の中心点と各４胞の中心点との前記関数ｆの片側差分値の２倍の値をそれぞれＤｆ東、Ｄｆ西、Ｄｆ南、Ｄｆ北とし、Ｄｆ東とＤｆ西が異符号の場合はゼロ、同符号の場合はその絶対値の小さいものをＤｆ東西とし、Ｄｆ南とＤｆ北が異符号の場合はゼロ、同符号の場合はその絶対値の小さいものをＤｆ南北とし、Ｄｆ東西とＤｆ南北の傾きを持ち、胞甲の中心におけるｆの値を通る平面を構成し、それを乙と呼ぶと、胞甲の北東方向の格子点における乙の値が胞甲の中心の値より大きい場合は、乙の該北東方向の格子点での値が、該格子点を共有する胞甲以外の３胞の中心でのｆ値を超えないようにＤｆ東西及びＤｆ南北を１以下の最大の定数により定数倍することにより乙を変更する。また、胞甲の北東方向の格子点における乙の値が胞甲の中心の値より小さい場合は、乙の該北東方向の格子点での値が、該格子点を共有する胞甲以外の３胞の中心でのｆ値を下回らないようにＤｆ東西及びＤｆ南北を１以下の最大の定数により定数倍することにより乙を変更する。このようにして得られた平面乙に対して、該北東方向の格子点に対して行った操作を、北西方向の格子点、南西方向の格子点及び南東方向の格子点に対してそれぞれ行うことにより、乙の傾きを変更させ、得られた平面を該胞甲内での前記補間関数Ｈとして定める。このアルゴリズムと上記のアルゴリズムは等価である。
【００７０】
＜第２実施形態の補間方法を利用した移流計算の流れ＞
図５は、上記した第２実施形態の補間方法を利用した移流計算の流れ図である。
【００７１】
まず、ステップＳ１１で初期設定をして、例えば、格子の生成、時間間隔を決定した後に関数の初期形状Ｆ（・，・，１）を初期化する。
【００７２】
次のステップＳ１２では、上述のＱ（・，・，Ｊ）の計算を行う（上記式（２８）、式（２９）、式（３０）、式（３１））。この時、上記式（２８）及び式（２９）に示すように、本発明の補間関数Ｈを利用する。そして、ステップＳ１３でＱ（・，・，Ｊ）のデータを基に、Ｆ（・，・，Ｊ＋１）の計算を行う（上記式（２２））。
【００７３】
その後はステップＳ１４で時刻を更新する。すなわち、今解こうとしている問題は上記式（１８）の時間発展方程式であるが、通常の時間発展方程式の数値解法の慣例に従い、時刻Ｔ＝ΔｔＪでのＦの値Ｆ（・，・，Ｊ）を、Ｆ（・，・，Ｊ−１）から求めた後に、時間をＴ＝Δｔ＊（Ｊ−１）からＴ＝Δｔ＊Ｊに変更を行う。
【００７４】
そして、ステップＳ１５で最終ステップか否かの判断を行い、最終でなければ、ステップＳ１２に戻り、最終ステップになるまで計算を行う。
【００７５】
本実施形態の補間方法を利用すると、関数のスロープがゼロのものとシャープなものの２つの領域に分離され、２相の非圧縮性流体の形状表現関数としては非常に良い結果を得ることができる。その結果を表示したのが図６（ａ），（ｂ）である。従来例の計算結果である図１０（ａ），（ｂ）に対応した同じ条件で計算を行っており、初期値と１０００ステップ後の形状を図示している。図６（ｂ）で示す本実施形態の移流計算結果は、図１０（ｂ）に示したものに比べて非常にシャープな矩形形状を保っていることが分かる。
【００７６】
＜本実施形態のハードウェア＞
本実施形態の特徴を成す補間方法を利用した移流計算を実行する計算機は、図３及び図４に示したものと同様である。
【００７７】
［第３実施形態］
第３実施形態では、上記第２実施形態で説明した二次元の格子に対する補間アルゴリズムを三次元に拡張したものである。本実施形態で取り扱うものは、
∂　Ｆ（ｘ，ｙ，ｚ，ｔ）　／　∂ｔ　＝∂　Ｃｘ（ｘ，ｙ，ｚ，ｔ）　Ｆ（ｘ，ｙ，ｚ，ｔ）　／　∂ｘ
＋　∂　Ｃｙ（ｘ，ｙ，ｚ，ｔ）　Ｆ（ｘ，ｙ，ｚ，ｔ）　／　∂ｙ
＋　∂　Ｃｚ（ｘ，ｙ，ｚ，ｔ）　Ｆ（ｘ，ｙ，ｚ，ｔ）　／　∂ｚ…（４０）
とする３次元空間の偏微分方程式であり、これを差分化して計算を行う問題を考える。Ｆ（ｘ，ｙ，ｚ，ｔ）を離散化したものとしてＦ（Ｘ，Ｙ，Ｚ，Ｔ）を考える。Ｆの引数Ｘ，Ｙ，Ｚ，Ｔは、以下離散的な空間格子及び離散時間を意味するものとする。
【００７８】
まず、本実施形態の補間関数が、区分的１次関数
2004062749.tif 000025
として、既に得られているとして、上記の偏微分方程式の差分近似解を求める方法を述べる。本実施形態の補間関数Ｈ［Ｉ，Ｊ，Ｋ，Ｔ］（ｘ，ｙ，ｚ）の係数Ｇｘ（Ｉ，Ｊ，Ｋ，Ｔ）、Ｇｙ（Ｉ，Ｊ，Ｋ，Ｔ）及びＧｚ（Ｉ，Ｊ，Ｋ，Ｔ）の求め方に関しては後で述べる。そこで述べるように、Ｇｘ（Ｉ，Ｊ，Ｋ，Ｔ）、Ｇｙ（Ｉ，Ｊ，Ｋ，Ｔ）及びＧｚ（Ｉ，Ｊ，Ｋ，Ｔ）の決定方法が本実施形態の補間アルゴリズムの特徴である。Ｇｘ（Ｉ，Ｊ，Ｔ）及びＧｙ（Ｉ，Ｊ，Ｔ）の決定に際し、格子の幾何学的性質である隣接情報を利用することによって決定を行っていることは重要である。
【００７９】
要するに、３次元実領域上の構成された３次元構造格子と、該３次元構造格子の各胞中心において値が定まる事によって定義されるある時刻Ｔでの該３次元構造格子上の関数Ｆ（Ｘ，Ｙ，Ｚ，Ｔ）とが用意されたとき、該関数Ｆ（Ｘ，Ｙ，Ｚ，Ｔ）と該格子の幾何を利用することによって該関数Ｆ（Ｘ，Ｙ，Ｚ，Ｔ）を拡張・補間して該３次元実領域上の関数Ｈ［Ｉ，Ｊ，Ｋ，Ｔ］（ｘ，ｙ，ｚ）を定義するものである。
【００８０】
もしも関数Ｃｘ（Ｉ＋１／２，Ｊ，Ｋ，Ｔ）が正であれば、
2004062749.tif 000026
とし、そうでなければ、
2004062749.tif 000027
と定義し、
また、もしも関数Ｃｙ（Ｉ，Ｊ＋１／２，Ｋ，Ｔ）が正であれば、
2004062749.tif 000028
とし、そうでなければ、
2004062749.tif 000029
と定義し、
また、もしも関数Ｃｚ（Ｉ，Ｊ，Ｋ＋１／２，Ｔ）が正であれば、
2004062749.tif 000030
とし、そうでなければ、
2004062749.tif 000031
と定義し、
2004062749.tif 000032
とする。
【００８１】
他方、区分的１次関数において勾配を表すスロープ関数Ｇ（Ｉ，Ｊ，Ｋ，Ｔ）を以下のように定める。
【００８２】
まず、
2004062749.tif 000033
を導入する。本実施形態の１次元制限関数として、
2004062749.tif 000034
を導入する。
【００８３】
次に、
2004062749.tif 000035
2004062749.tif 000036
2004062749.tif 000037
を定義する。
【００８４】
更に、
ｈｐｐｐ（Ｉ，Ｊ，Ｋ，Ｔ）　＝　　Ｓｘ（Ｉ，Ｊ，Ｋ，　Ｔ）／２．　＋　　Ｓｘ（Ｉ，Ｊ，Ｋ，　Ｔ）／２　＋　Ｓｚ（Ｉ，Ｊ，Ｋ，Ｔ）／２．
ｈｍｐｐ（Ｉ，Ｊ，Ｋ，Ｔ）　＝　−　Ｓｘ（Ｉ，Ｊ，Ｋ，　Ｔ）／２．　＋　　Ｓｘ（Ｉ，Ｊ，Ｋ，　Ｔ）／２　＋　Ｓｚ（Ｉ，Ｊ，Ｋ，　Ｔ）／２．
ｈｐｍｐ（Ｉ，Ｊ，Ｋ，Ｔ）　＝　　　Ｓｘ（Ｉ，Ｊ，Ｋ，　Ｔ）／２．　−　　Ｓｘ（Ｉ，Ｊ，Ｋ，　Ｔ）／２　＋　Ｓｚ（Ｉ，Ｊ，Ｋ，　Ｔ）／２．
ｈｍｍｐ（Ｉ，Ｊ，Ｋ，Ｔ）　＝　−　Ｓｘ（Ｉ，Ｊ，Ｋ，　Ｔ）／２．　−　　Ｓｘ（Ｉ，Ｊ，Ｋ，　Ｔ）／２　＋　Ｓｚ（Ｉ，Ｊ，Ｋ，　Ｔ）／２．
ｈｐｐｍ（Ｉ，Ｊ，Ｋ，Ｔ）　＝　　　Ｓｘ（Ｉ，Ｊ，Ｋ，　Ｔ）／２．　＋　　Ｓｘ（Ｉ，Ｊ，Ｋ，　Ｔ）／２　−　Ｓｚ（Ｉ，Ｊ，Ｋ，　Ｔ）／２．
ｈｍｐｍ（Ｉ，Ｊ，Ｋ，Ｔ）　＝　−　Ｓｘ（Ｉ，Ｊ，Ｋ，　Ｔ）／２．　＋　　Ｓｘ（Ｉ，Ｊ，Ｋ，　Ｔ）／２　−　Ｓｚ（Ｉ，Ｊ，Ｋ，　Ｔ）／２．
ｈｐｍｍ（Ｉ，Ｊ，Ｋ，Ｔ）　＝　　　Ｓｘ（Ｉ，Ｊ，Ｋ，　Ｔ）／２．　−　　Ｓｘ（Ｉ，Ｊ，Ｋ，　Ｔ）／２　−　Ｓｚ（Ｉ，Ｊ，Ｋ，　Ｔ）／２．
ｈｍｍｍ（Ｉ，Ｊ，Ｋ，Ｔ）　＝　−　Ｓｘ（Ｉ，Ｊ，Ｋ，　Ｔ）／２．　−　　Ｓｘ（Ｉ，Ｊ，Ｋ，　Ｔ）／２　−　Ｓｚ（Ｉ，Ｊ，Ｋ，　Ｔ）／２．…（５３）
とする。
【００８５】
更に、「ｐｈｉｐｐｐ＝ｐｈｉｍｐｐ＝ｐｈｉｐｍｐ＝ｐｈｉｍｍｐ＝　ｐｈｉｐｐｍ＝ｐｈｉｍｐｍ＝ｐｈｉｐｍｍ＝ｐｈｉｍｍｍ＝１としておき、もし（ｈｐｐｐ　＞　Ｖｐｐｐｍａｘ　）ならば、　ｐｈｉｐｐｐ　＝　Ｖｐｐｐｍａｘ／ｈｐｐｐと置き直す。もし（ｈｐｐｐ　＜　Ｖｐｐｐｍｉｎ　）ならば、　ｐｈｉｐｐｐ　＝　Ｖｐｐｐｍｉｎ／ｈｐｐｐと置き直す。同様に、もし（ｈｐｍｐ　＞　Ｖｐｍｐｍａｘ　）ならば、　ｐｈｉｐｍｐ　＝　Ｖｐｍｐｍａｘ／ｈｐｍｐと置き直す。
またもし、（ｈｐｍ　＜　Ｖｐｍｍｉｎ　）ならば、　ｐｈｉｐｍ　＝　Ｖｐｍｐｍｉｎ／ｈｐｍと置き直す。もし（ｈｍｐｐ　＞　Ｖｍｐｐｍａｘ　）ならば、　ｐｈｉｍｐｐ　＝　Ｖｍｐｐｍａｘ／ｈｍｐｐと置き直す。またもし、（ｈｍｐｐ　＜　Ｖｍｐｐｍｉｎ　）ならば、　ｐｈｉｍｐｐ　＝　Ｖｍｐｐｍｉｎ／ｈｍｐｐと置き直す。もし（ｈｍｍｐ
＞　Ｖｍｍｐｍａｘ　）ならば、　ｐｈｉｍｍｐ　＝　Ｖｍｍｐｍａｘ／ｈｍｍｐと置き直す。
【００８６】
またもし、（ｈｍｍｐ　＜　Ｖｍｍｐｍｉｎ　）ならば、　ｐｈｉｍｍｐ　＝　Ｖｍｍｐｍｉｎ／ｈｍｍｐと置き直す。もし（ｈｐｐｍ　＞　Ｖｐｐｍｍａｘ　）ならば、　ｐｈｉｐｐｍ　＝　Ｖｐｐｍｍａｘ／ｈｐｐｍと置き直す。またもし（ｈｐｐｍ　＜　Ｖｐｐｍｍｉｎ　）ならば、　ｐｈｉｐｐｍ　＝　Ｖｐｐｍｍｉｎ／ｈｐｐｍと置き直す。もし（ｈｐｍｍ　＞　Ｖｐｍｍｍａｘ　）ならば、　ｐｈｉｐｍｍ　＝　Ｖｐｍｍｍａｘ／ｈｐｍｍと置き直す。またもし、（ｈｐｍｍ　＜　Ｖｐｍｍｍｉｎ　）ならば、　ｐｈｉｐｍｍ　＝　Ｖｐｍｍｍｉｎ／ｈｐｍｍと置き直す。もし（ｈｍｐｍ　＞　Ｖｍｐｍｍａｘ　）ならば、　ｐｈｉｍｐｐ　＝　Ｖｍｐｐｍａｘ／ｈｍｐｐと置き直す。
【００８７】
またもし、（ｈｍｐｍ　＜　Ｖｍｐｍｍｉｎ　）ならば、　ｐｈｉｍｐｍ　＝　Ｖｍｐｍｍｉｎ／ｈｍｐｍと置き直す。もし（ｈｍｍｍ　＞　Ｖｍｍｍｍａｘ　）ならば、　ｐｈｉｍｍｍ　＝　Ｖｍｍｍｍａｘ／ｈｍｍｍと置き直す。またもし、（ｈｍｍｍ　＜　Ｖｍｍｍｍｉｎ　）ならば、　ｐｈｉｍｍｍ　＝　Ｖｍｍｍｍｉｎ／ｈｍｍｍと置き直す。」という計算を行う。
【００８８】
このようにして得られたｐｈｉｐｐｐ，ｐｈｉｍｐｐ，ｐｈｉｐｍｐ，ｐｈｉｍｍｐ，　ｐｈｉｐｐｍ，ｐｈｉｍｐｍ，ｐｈｉｐｍｍ，ｐｈｉｍｍｍに対して、
Ｇｘ（Ｉ，Ｊ，Ｋ，Ｔ）＝Ｍｉｎ（ｐｈｉｐｐｐ，ｐｈｉｍｐｐ，ｐｈｉｐｍｐ，ｐｈｉｍｍｐ，ｐｈｉｐｐｍ，ｐｈｉｍｐｍ，ｐｈｉｐｍｍ，ｐｈｉｍｍｍ）Ｓｘ（Ｉ，Ｊ，Ｋ，Ｔ）／Δｘ
Ｇｙ（Ｉ，Ｊ，Ｋ，Ｔ）＝Ｍｉｎ（ｐｈｉｐｐｐ，ｐｈｉｍｐｐ，ｐｈｉｐｍｐ，ｐｈｉｍｍｐ，ｐｈｉｐｐｍ，ｐｈｉｍｐｍ，ｐｈｉｐｍｍ，ｐｈｉｍｍｍ）Ｓｙ（Ｉ，Ｊ，Ｋ，Ｔ）／Δｙ
Ｇｚ（Ｉ，Ｊ，Ｋ，Ｔ）＝Ｍｉｎ（ｐｈｉｐｐｐ，ｐｈｉｍｐｐ，ｐｈｉｐｍｐ，ｐｈｉｍｍｐ，ｐｈｉｐｐｍ，ｐｈｉｍｐｍ，ｐｈｉｐｍｍ，ｐｈｉｍｍｍ）Ｓｚ（Ｉ，Ｊ，Ｋ，Ｔ）／Δｚ
とした後に、（Ｉ，Ｊ，Ｋ）格子内の補間関数Ｈ（Ｉ，Ｊ，Ｋ，Ｔ）（ｘ，ｙ，ｚ）を
2004062749.tif 000038
とするものである。
【００８９】
このように本実施形態の補間方法をまとめると、着目する格子を胞甲とし、胞甲の各辺を東西あるいは南北あるいは天地にそれぞれ合わせた際に東西南北天地方向に隣合う６胞に対して、胞甲の中心点と各６胞の中心点との前記関数ｆの片側差分値の２倍の値をそれぞれ　Ｄｆ東、Ｄｆ西、Ｄｆ南、Ｄｆ北、Ｄｆ天、Ｄｆ地とし、Ｄｆ東とＤｆ西が異符号の場合はゼロ、同符号の場合はその絶対値の小さいものをＤｆ東西とし、Ｄｆ南とＤｆ北が異符号の場合はゼロ、同符号の場合はその絶対値の小さいものをＤｆ南北とし、Ｄｆ天とＤｆ地が異符号の場合はゼロ、同符号の場合はその絶対値の小さいものをＤｆ天地とし、Ｄｆ東西とＤｆ南北とＤｆ天地の傾きを持ち、胞甲の中心におけるｆの値を通る補間関数の候補を構成し、それを乙と呼ぶと、胞甲の北東天方向の格子点における乙の値が胞甲の中心の値より大きい場合は、乙の該北東天方向の格子点での値が、該格子点を共有する胞甲以外の７胞の中心でのｆ値を超えないように、Ｄｆ東西及びＤｆ南北及びＤｆ天地を１以下の最大の定数により定数倍することにより乙を変更する。また、胞甲の北東天方向の格子点における乙の値が胞甲の中心の値より小さい場合は、乙の該北東天方向の格子点での値が、該格子点を共有する胞甲以外の７胞の中心でのｆ値を下回らないように、Ｄｆ東西及びＤｆ南北及びＤｆ天地を１以下の最大の定数により定数倍することにより乙を変更する。
【００９０】
このようにして得られた補間関数の候補乙に対して、該北東天方向の格子点に対して行った操作を、北西天方向の格子点、南西天方向の格子点、南東天方向の格子点、北東地方向の格子点、北西地方向の格子点、南西地方向の格子点、及び南東地方向の格子点に対してそれぞれ行うことにより、乙の傾きを変更させ、得られた平面を該胞甲内での前記補間関数Ｆとして定める事により補間関数を得る。このアルゴリズムと上記に示したアルゴリズムと等価である。
【００９１】
上記した本実施形態の補間方法を利用した移流計算の流れ図は、図５に示したものと同様である。また、本第３実施形態の補間方法を図３及び図４で説明したと同様にコンピュータ上で実現することができる。本実施形態においても、前述した第２実施形態における２次元のものと同様に形状をシャープに保ったまま移流計算が行うことができる。
【００９２】
なお、上記第２及び第３実施形態の補間方法を利用すると、不連続点を含みながらも格子点上に与えられた値から２次元あるいは３次元空間上の関数を構成することができる。このことより、本発明の補間方法を上記実施形態とは異なる状況、例えば画像データの再現、２次元あるいは３次元のデジタル対象のデータの詳細化等々で使用しても、同様にシャープな形状を表現するものとなる。
【００９３】
本発明は、上述した実施形態の装置に限定されず、複数の機器から構成されるシステムに適用しても、１つの機器から成る装置に適用してもよい。前述した実施形態の機能を実現するソフトウェアのプログラムコードを記憶した記憶媒体をシステムあるいは装置に供給し、そのシステムあるいは装置のコンピュータ（またはＣＰＵやＭＰＵ）が記憶媒体に格納されたプログラムコードを読み出し実行することによっても、完成されることは言うまでもない。
【００９４】
この場合、記憶媒体から読み出されたプログラムコード自体が前述した実施形態の機能を実現することになり、そのプログラムコードを記憶した記憶媒体は本発明を構成することになる。プログラムコードを供給するための記憶媒体としては、例えば、フロッピー（登録商標）ディスク、ハードディスク、光ディスク、光磁気ディスク、ＣＤ−ＲＯＭ、ＣＤ−Ｒ、磁気テープ、不揮発性のメモリカード、ＲＯＭを用いることができる。また、コンピュータが読み出したプログラムコードを実行することにより、前述した実施形態の機能が実現されるだけではなく、そのプログラムコードの指示に基づき、コンピュータ上で稼動しているＯＳなどが実際の処理の一部または全部を行い、その処理によって前述した実施形態の機能が実現される場合も含まれることは言うまでもない。
【００９５】
さらに、記憶媒体から読み出されたプログラムコードが、コンピュータに挿入された機能拡張ボードやコンピュータに接続された機能拡張ユニットに備わるメモリに書き込まれた後、次のプログラムコードの指示に基づき、その拡張機能を拡張ボードや拡張ユニットに備わるＣＰＵなどが処理を行って実際の処理の一部または全部を行い、その処理によって前述した実施形態の機能が実現される場合も含まれることは言うまでもない。
【００９６】
【発明の効果】
以上詳細に説明したように本発明によれば、２相の非圧縮性流体における流体体積法（ＶＯＦ法）において、形状をシャープに保持させたまま時間発展させることが可能な補間アルゴリズムを実現することが可能であり、２相の非圧縮性流体の形状表現関数としては非常に良い結果を得ることができる。
【図面の簡単な説明】
【図１】第１実施形態の補間方法を利用した移流計算の流れ図である。
【図２】第１実施形態の補間方法を利用した移流計算の結果を示す図である。
【図３】第１実施形態の特徴を成す補間方法を利用した移流計算を実行する計算機の概観図である。
【図４】本発明の補間方法を利用した移流計算用の計算機の要部構成を示すブロック図である。
【図５】第２実施形態の補間方法を利用した移流計算の流れ図である。
【図６】第２及び第３実施形態の補間方法を利用した移流計算の結果を示す図である。
【図７】従来の補間方法を利用した移流計算の結果を示す図である。
【図８】従来の補間方法を利用した移流計算の結果を示す図である。
【図９】従来の補間方法を利用した移流計算の結果を示す図である。
【図１０】従来の補間方法（高次元）を利用した移流計算の結果を示す図である。
【符号の説明】
１１０　フロッピー（登録商標）ディスク
１１１　計算機の本体
１１２　ディスプレイ
１１３　フロッピー（登録商標）ディスクドライブ、
１１４　マウス
１１５　キーボード
１１６　ＣＰＵ、
１１７　ハードディスク
１１８　メモリ
【図１】
2004062749.tif 000039
【図２】
2004062749.tif 000040
【図３】
2004062749.tif 000041
【図４】
2004062749.tif 000042
【図５】
2004062749.tif 000043
【図６】
2004062749.tif 000044
【図７】
2004062749.tif 000045
【図８】
2004062749.tif 000046
【図９】
2004062749.tif 000047
【図１０】
2004062749.tif 000048