JP3978534B2 - Method for setting moving boundary moving on fixed grid and computer program for realizing the same - Google Patents

Description

【００１２】
また、レベルセット関数ｆ_lsの時間発展は、次式（２）の移流方程式に従って計算される。なお、次式（２）において、ｕは流体の速度を表している。
【００１３】
【数２】

【００１４】
ここで、図９に示すような１次元の流体問題を例にとると、上式（１）に示すレベルセット関数ｆ_lsは、符号２０に示すようなものとなり、液体と気体との境界２１，２２がレベルセット関数ｆ_lsのゼロ等高面（ゼロレベルセット）として求められる。
【００１５】
このため、上式（１）に示すレベルセット関数ｆ_lsを上式（２）の移流方程式に従って時間発展させ、ある時刻でのレベルセット関数ｆ_lsのゼロレベルセットを求めることにより、境界を追跡することができる。
【００１６】
なお、上式（１）および図９から分かるように、レベルセット関数ｆ_lsは傾きが１の関数であり、ある点２３におけるｆ_lsの値ａは、点２３の境界２１からの距離ｂと等しい。すなわち、レベルセット関数ｆ_lsは、境界からの距離を表す距離関数としての性質を有している。
【００１７】
このようなレベルセット法では、境界を直接追跡することなく、レベルセット関数ｆ_lsを用いて間接的に追跡するので、境界のトポロジーが変化するような場合（例えば、図１０(a)(b)(c)に示すように、２つの水滴（液体）が衝突するような場合）でも、特別な処理を行うことなく自然に対応することができる。また、直交固定格子上で曲面からなる境界を精度良く捕らえるには、境界をある程度スムースにする必要があるが、レベルセット法では、境界のスムーシングの幅を容易にコントロールすることができるので、曲面を精度良く捕らえることができる。さらに、表面張力の計算等に必要な曲率κの計算には、境界に対する単位法線ベクトルの発散を求める必要があり、境界の周辺での単位法線ベクトルが必要になるが、レベルセット法では、レベルセット関数ｆ_lsを用いることにより、次式（３）に従って曲率κを精度良く計算することができる。
【００１８】
【数３】

【００１９】
なお、このようにして求められた曲率κに基づいて、ＣＳＦ（Continum Surface Force）モデル（文献４（J.U. Brackbill, D.B. Kothe and C. Zemach, J. Comput. Phys. 100, 335 (1992)）参照）等を用いて、容易に表面張力を計算することができる。
【００２０】
【発明が解決しようとする課題】
上述したように、レベルセット法は、移動境界を設定（捕獲）する手法として優れたものであるが、移動境界の設定に必要とされるレベルセット関数ｆ_lsには一般に、上式（２）の移流方程式に従って時間発展させるにつれてその関数の性質（特に、距離関数としての性質）が損なわれていくという問題がある。
【００２１】
すなわち、ある時刻でのレベルセット関数ｆ_lsが図１１(a)に示すように距離関数の性質を満たしていても、上式（２）の移流方程式により時間発展させた後のレベルセット関数ｆ_lsは、図１１(b)に示すように一般に距離関数の性質を満たさなくなる。
【００２２】
このため、従来のレベルセット法では、時間発展させた後のゼロレベルセットが境界を正しく表現しているものと仮定し、このゼロレベルセットを基準としてその他の部分に対して距離関数の性質を回復させるための処理、すなわち再初期化処理（re-initialization）を行っている。
【００２３】
しかしながら、このような再初期化処理に伴って、ゼロレベルセットも若干ずれることとなるので、これにより境界の位置に関してエラーが混入する（図１１(c)）。このエラーは１計算ステップ分で見れば僅かであるが、このエラーが混入したレベルセット関数（図１１(d)）を用いて次回以降の計算ステップを繰り返すこととなるので、このようなエラーが毎計算ステップごとに蓄積してしまう。このため、時間発展させるにつれて境界の位置が正しい位置から大きくずれてしまう。このようなエラーの蓄積は、流体問題の種類によっては深刻な問題となり、特に複雑な移動境界を有する流体問題の場合には、このエラーは無視できないものとなる。
【００２４】
本発明はこのような点を考慮してなされたものであり、レベルセット法において問題となる再初期化処理時のエラーの蓄積を効果的に防止して、固定格子上で移動境界を精度良く捕らえることができる、移動境界の設定方法およびそれを実現するコンピュータプログラムを提供することを目的とする。
【００２５】
【課題を解決するための手段】
本発明は、その第１の解決手段として、移動境界を有する流体現象を数値的に解析する流体シミュレーションで用いられ、固定格子上を移動する移動境界をレベルセット法を用いて設定する、移動境界の設定方法において、移流方程式に従って移動する移動境界の設定に必要とされるレベルセット関数から独立した関数として、移動境界を表すゼロレベルセットを追跡することを目的とするゼロレベルセット追跡関数であって、そのゼロ等高面がゼロレベルセットであり、レベルセット関数と同一の形を有するゼロレベルセット追跡関数を準備するステップと、前記ゼロレベルセット追跡関数を前記移流方程式に従って、当該移流方程式の計算を安定的に行うために要求される安定性の条件から決定される時間間隔分だけ時間発展させる時間発展ステップと、前記時間発展ステップにより時間発展された後のゼロレベルセット追跡関数に基づいて、移動境界を表すゼロレベルセットを求めるとともに、この求められたゼロレベルセットに基づいて、移動境界の設定に必要とされるレベルセット関数を構築する境界設定ステップとを含み、前記時間発展ステップおよび前記境界設定ステップを、所定の終了時間に達するまで繰り返すことを特徴とする、移動境界の設定方法を提供する。
【００２６】
なお、上述した第１の解決手段においては、前記境界設定ステップにより構築されたレベルセット関数に基づいて、移動境界の内外にある物質の識別に用いられる密度関数を求めるステップをさらに含むことが好ましい。また、前記時間発展ステップにおいて、前記移流方程式をＣＩＰ法により解くことが好ましい。さらに、前記固定格子は直交固定格子であることが好ましい。
【００２７】
本発明は、その第２の解決手段として、移動境界を有する流体現象を数値的に解析する流体シミュレーションで用いられ、固定格子上を移動する移動境界をレベルセット法を用いて設定する、移動境界の設定方法を実現するコンピュータプログラムにおいて、コンピュータに対して実行させる手順として、移流方程式に従って移動する移動境界の設定に必要とされるレベルセット関数から独立した関数として、移動境界を表すゼロレベルセットを追跡することを目的とするゼロレベルセット追跡関数であって、そのゼロ等高面がゼロレベルセットであり、レベルセット関数と同一の形を有するゼロレベルセット追跡関数を準備する手順と、前記ゼロレベルセット追跡関数を前記移流方程式に従って、当該移流方程式の計算を安定的に行うために要求される安定性の条件から決定される時間間隔分だけ時間発展させる時間発展手順と、前記時間発展手順により時間発展された後のゼロレベルセット追跡関数に基づいて、移動境界を表すゼロレベルセットを求めるとともに、この求められたゼロレベルセットに基づいて、移動境界の設定に必要とされるレベルセット関数を構築する境界設定手順とを含み、前記時間発展手順および前記境界設定手順を、所定の終了時間に達するまで繰り返しコンピュータに対して実行させることを特徴とする、移動境界の設定方法を実現するコンピュータプログラムを提供する。
【００２８】
本発明によれば、レベルセット法で用いられるレベルセット関数ｆ_lsから独立した関数としてゼロレベルセット追跡関数を準備し、このゼロレベルセット追跡関数を移流方程式に従って時間発展させ、ゼロレベルセット追跡関数に基づいて、移動境界を表すゼロレベルセットを求めるとともに、この求められたゼロレベルセットに基づいて、移動境界の設定に必要とされるレベルセット関数を構築しているので、時間発展の計算によりレベルセット関数にエラーが蓄積することがなく、従来のレベルセット法で問題となっていた再初期化処理時のエラーの蓄積の問題を効果的に防止することができる。
【００２９】
【発明の実施の形態】
以下、図面を参照して本発明の実施の形態について説明する。ここで、図１は本発明による移動境界の設定方法の一実施の形態を説明するためのフローチャート、図２(a)(b)(c)(d)は本実施の形態の概要を１次元の移動境界問題を例に挙げて説明するための概念図である。なお、本実施の形態に係る移動境界の設定方法は、直交固定格子上を移動する移動境界をレベルセット法を用いて設定するものであり、移動境界を有する流体現象等の解析に用いられる。なお、固定格子としては、直交固定格子以外にも、三角格子等の任意の固定格子を用いることができる。
【００３０】
図１に示すように、本実施の形態に係る移動境界の設定方法においては、まず、レベルセット法で用いられるレベルセット関数ｆ_lsから独立した関数として、移動境界を表すゼロレベルセットを追跡することのみを目的とするゼロレベルセット追跡関数ｆ_zeroを準備する（ステップ１０１）。なお、広い意味でいえば、ゼロレベルセット追跡関数ｆ_zeroもレベルセット関数であり、そのゼロ等高面はゼロレベルセットであるが、本明細書では、距離関数の性質を持つ関数のみをレベルセット関数と呼び、このようなレベルセット関数のゼロ等高面のみをゼロレベルセットと呼ぶことにする。
【００３１】
次に、ゼロレベルセット追跡関数ｆ_zeroを次式（４）の移流方程式に従って、安定性の条件から決定される時間間隔分だけ時間発展させる（ステップ１０２）。具体的には例えば、図２(a)(b)に示すように、ゼロレベルセット追跡関数ｆ_zeroを境界の移動に合わせて右方向に移動させる。なおこのとき、ゼロレベルセット追跡関数ｆ_zeroは時間発展に伴ってその形状が崩れる。
【００３２】
【数４】

【００３３】
なお、上式（４）の移流方程式の解法には、高精度でかつロバストな方法であるＣＩＰ法（Cubic Interpolated Propagation Method）（文献５（T. Yabe and T. Aoki, Comput. Phys. Commun. 66, 219 (1991)）および文献６（T. Nakamura and T. Yabe, Comput. Phys. Commun. 120, 122 (1999)）参照）を用いることが好ましく、これにより、ゼロレベルセット追跡関数ｆ_zeroの形状の崩れを抑えることができる。なお、移流方程式の解法には、ＣＩＰ法以外にも、ＰＰＭ（Piecewise Parabolic Method）法等の任意の方法を用いることができる。
【００３４】
その後、時間発展された後のゼロレベルセット追跡関数ｆ_zeroに基づいて、移動境界を表すゼロレベルセットを求めるとともに、この求められたゼロレベルセットに基づいて、移動境界の設定に必要とされるレベルセット関数ｆ_lsを構築する（ステップ１０３）。具体的には例えば、図２(c)に示すように、ゼロレベルセット追跡関数ｆ_zeroのゼロ等高面に基づいて、移動境界を表すゼロレベルセットを求め、この求められたゼロレベルセットに基づいて、距離関数の性質を持つレベルセット関数ｆ_lsを構築する。なお、レベルセット関数ｆ_lsは、境界のスムーシングや曲率の計算等のために用いられる。
【００３５】
以下、ステップ１０３における処理（ゼロレベルセット追跡関数ｆ_zeroからレベルセット関数ｆ_lsを構築する処理）の詳細について説明する。
【００３６】
ステップ１０３においては、まず、ゼロレベルセット追跡関数ｆ_zeroのゼロ等高面を、各格子点間を補間することにより求める。
【００３７】
次に、このようにして求められたゼロ等高面を、移動境界を表すゼロレベルセットとし、境界から１セル以内の格子点のレベルセット関数ｆ_lsを、ファーストマーチング法（fast marching method）（上記文献２および文献７（D. Adalsteinsson and J.A. Sethian, J. Comput. Phys. 148, 2 (1999)）参照）に従って、次式（５）のEikonal方程式を計算することにより求める。
【００３８】
【数５】

【００３９】
次に、このようにして求められたレベルセット関数ｆ_lsの値を固定し、Sussmann等の反復計算による方法（文献８（M. Sussman, P. Smereka and S. Osher, J. Comput. Phys. 114, 146 (1994)）参照）に従って、次式（６）が収束するまで収束計算を行うことにより、１セル以遠の部分のレベルセット関数ｆ_lsを構築する。なお、Sussman等の方法では一般に、反復計算の間にゼロレベルセットがずれるという問題があるが、ここでは、ゼロレベルセットの周辺の値がファーストマーチング法により事前に計算され、実際の収束計算時には固定されているため、その問題は起こらない。
【００４０】
【数６】

【００４１】
ここで、τは反復計算のための時間であり、実時間と同じである必要はない。また、Ｓは平滑化符号関数（smoothed sign function）であり、次式（７）により表される。
【００４２】
【数７】

【００４３】
なお、ゼロレベルセットの周辺のレベルセット関数ｆ_lsを固定しているため、計算の安定性が問題とならない場合には、Ｓ（ｆ_ls）＝１とすることも可能である。
【００４４】
ここで、上式（６）による収束計算における初期値としては、前の計算ステップのｆ_lsの値や、１、−１等を用いるのではなく、上式（２）の移流方程式により計算した値を用いるようにしてもよく、これにより、収束に要する時間を大幅に短縮することができる。このような収束計算には、上述したＣＩＰ法を用いることができる他、１次の風上差分法を用いることができる。これは、１次の風上差分法のような比較的低精度の方法を用いても、その方法により混入されたエラーが上式（６）による収束計算により補正されるからである。なお、計算時間に拘らない場合には、上式（６）による収束計算ではなく、上述したファーストマーチング法を用いることもできる。
【００４５】
その後、このようにして構築されたレベルセット関数ｆ_lsに基づいて、移動境界の内外にある物質の識別に用いられる密度関数（color (density) function）ｆ_colorを求める（ステップ１０４）。なお、密度関数ｆ_colorは、次式（８）（９）に従って、レベルセット関数ｆ_lsから平滑化ヘビサイド関数（smoothed Heaviside function）Ｈ_αを用いて計算することができる。ここで、２αは境界のスムーシングの幅である。
【００４６】
【数８】

【００４７】
そして、所定の終了時間に達するまでステップ１０２〜１０４における処理を繰り返す（ステップ１０５）。なお、ステップ１０２に戻って処理を続ける場合、時間発展の計算の基礎となる関数はゼロレベルセット追跡関数ｆ_zeroである。すなわち、ステップ１０３において、ゼロレベルセット追跡関数ｆ_zeroに基づいてレベルセット関数ｆ_lsを構築するときには、従来のレベルセット法と同様にエラーが混入することとなるが（図２(c)参照）、ステップ１０２における時間発展の計算には、エラーが混入していないゼロレベルセット追跡関数ｆ_zeroが用いられる（図２(d)参照）。
【００４８】
このように本実施の形態によれば、レベルセット法で用いられるレベルセット関数ｆ_lsから独立した関数としてゼロレベルセット追跡関数ｆ_zeroを準備し、このゼロレベルセット追跡関数ｆ_zeroを移流方程式に従って時間発展させ、ゼロレベルセット追跡関数ｆ_zeroに基づいて、移動境界を表すゼロレベルセットを求めるとともに、この求められたゼロレベルセットに基づいて、移動境界の設定に必要とされるレベルセット関数ｆ_lsを構築しているので、時間発展の計算によりレベルセット関数ｆ_lsにエラーが蓄積することがなく、従来のレベルセット法で問題となっていた再初期化処理時のエラーの蓄積の問題を効果的に防止することができる。
【００４９】
なお、上述した実施の形態に係る移動境界の設定方法は、図３に示すようなコンピュータシステム１０上で稼働するコンピュータプログラムとして実現することができる。ここで、コンピュータシステム１０は、バス１１と、バス１１に接続されたプロセッサ１２、メモリ１３およびハードディスク１４と、バス１１に接続された周辺機器（キーボードやマウス等の入力装置１５、ディスプレイやプリンタ等の出力装置１６、ＦＤドライブ１７およびＣＤ−ＲＯＭドライブ１８）とを備えている。そして、上述したようなコンピュータプログラムは、フレキシブルディスク１９ａおよびＣＤ−ＲＯＭ１９ｂ等のようなコンピュータ読み取り可能な記録媒体に格納され、メモリ１３やハードディスク１４等に読み込まれる。そして、コンピュータシステム１０において、プロセッサ１２から逐次読み出されて実行されることにより上述したような手順を実現することができる。
【００５０】
【実施例】
次に、上述した実施の形態の具体的実施例について述べる。
【００５１】
まず、上述した実施の形態をZalesak問題（文献９（S.T. Zalesak, J. Comput. Phys. 31, 335 (1979)）参照）に適用した場合の実施例について述べる。
【００５２】
図４に示すように、本実施例の移動境界問題は、一部に切り欠きを有する円形状の二次元の剛体３１を矢印３２に沿って反時計方向に回転させるものである。
【００５３】
ここで、剛体３１の初期条件は、次式（１０）とした。
【００５４】
【数９】

【００５５】
また、計算のための格子は１００×１００の直交固定格子とし、剛体３１の回転速度は、次式（１１）のような速度場の下で、１回転の計算を６２８回の計算ステップで行うようなものとした。（すなわち、時間間隔Δｔを１／６２８とした。）
【００５６】
【数１０】

【００５７】
さらに、レベルセット関数ｆ_lsは図５(a)に示すものとした。なお、このようなレベルセット関数ｆ_lsに対応する、初期状態でのゼロレベルセット（境界）および密度関数ｆ_colorはそれぞれ、図５(b)(c)に示すようなものとなる。
【００５８】
このような条件の下で、レベルセット関数ｆ_lsと同一の形を有するゼロレベルセット追跡関数ｆ_zeroを準備し、このゼロレベルセット追跡関数ｆ_zeroを上式（１１）に従ってＣＩＰ法により時間発展させ、ゼロレベルセット追跡関数ｆ_zeroに基づいて、移動境界を表すゼロレベルセットを求めるとともに、この求められたゼロレベルセットに基づいて、移動境界の設定に必要とされるレベルセット関数ｆ_lsを構築し、移動境界を設定した。
【００５９】
図６(a)(b)および図６(c)(d)はそれぞれ、剛体３１の１回転後および１０回転後の境界の様子を示す図である。なお、図６(a)(c)はレベルセット関数ｆ_lsに基づいて求められたゼロレベルセット（境界）を厳密解（初期条件）とともに示す図であり、図６(b)(d)はレベルセット関数ｆ_lsに基づいて求められた密度関数ｆ_colorを示す図である。
【００６０】
図６(a)(b)(c)(d)に示すように、本実施例では、剛体３１の１回転後はもちろん１０回転後であっても境界を精度良く捕らえていることが分かる。
【００６１】
ここで、比較例として、上述した従来のレベルセット法を用いて剛体３１を１回転させた場合の境界の様子を図７(a)(b)に示す。図７(a)はゼロレベルセット（境界）を厳密解（初期条件）とともに示す図であり、図７(b)は密度関数ｆ_colorを示す図である。なお、この比較例におけるレベルセット法では、レベルセット関数ｆ_lsを時間発展させる際の移流方程式の解法としてＣＩＰ法を用いた。
【００６２】
図７(a)(b)に示すように、この比較例では、１回転後にもかかわらず、移動境界が初期状態からかなりずれていることが分かる。なお、このずれは時間とともに拡大した。
【００６３】
次に、上述した実施の形態を流体問題に適用した場合の実施例について述べる。
【００６４】
図８(a)(b)に示すように、本実施例の流体問題は、長方形の剛体３３を水面３４に落下させるものである。なお、計算のための格子は４０×４０の直交固定格子とし、時間間隔（Δｔ）は安定性の条件を満たすように決定した。また、流体シミュレーションの手法としては、圧力のポアソン方程式をベースにした混相流の計算手法であるＣＵＰ法（文献１０（T. Yabe and P.Y. Wang, J. Phys. Soc. Jpn. 60, 2105 (1991)）参照）を用いた。
【００６５】
その結果、図８(b)に示すように、本実施例においては、比較的少ない格子数にもかかわらず、剛体３３の形状を良く捕らえていることが分かる。これに対し、図８(c)に示す比較例（従来のレベルセット法を用いたもの）では、再初期化処理時のエラーにより、剛体３３の角が丸くなる傾向があることが分かる。
【００６６】
【発明の効果】
以上説明したように本発明によれば、レベルセット法において問題となる再初期化処理時のエラーの蓄積を効果的に防止して、固定格子上で移動境界を精度良く捕らえることができる。従って、本発明は、河川の激流や、海岸の大規模波動、工業装置における流体現象等に対応する構造物の設計等に極めて有用である。
【図面の簡単な説明】
【図１】本発明による移動境界の設定方法の一実施の形態を説明するためのフローチャート。
【図２】図１に示す実施の形態の概要を１次元の移動境界問題を例に挙げて説明するための概念図。
【図３】図１に示す実施の形態を実現するコンピュータプログラムが実行されるコンピュータシステムの一例を示す図。
【図４】本発明の一実施例としてのZalesak問題（二次元の剛体の回転問題）の概要を説明するための模式図。
【図５】図４に示す実施例におけるレベルセット関数、ゼロレベルセットおよび密度関数の一例を示す図。
【図６】図４に示す実施例における剛体の回転後（１回転後および１０回転後）の境界の様子を示す図。
【図７】従来のレベルセット法を用いた場合の剛体の回転後（１回転後）の境界の様子を示す図。
【図８】本発明の一実施例としての流体問題（長方形の剛体を水面に落下させる問題）の計算結果を示す図。
【図９】１次元の流体問題で用いられるレベルセット関数の一例を示す図。
【図１０】境界のトポロジーが変化するような１次元の流体問題の一例を示す図。
【図１１】従来のレベルセット法の概要を説明するための概念図。
【符号の説明】
ｆ_zero ゼロレベルセット追跡関数
ｆ_ls レベルセット関数
１０ コンピュータシステム
２０ レベルセット関数
２１，２２ 境界[0001]
BACKGROUND OF THE INVENTION
The present invention relates to a method for numerically analyzing a fluid phenomenon or the like having a moving boundary, and in particular, a moving boundary setting method for setting a moving boundary moving on a fixed grid using a level set method, and the realization thereof. Relates to a computer program.
[0002]
[Prior art]
Liquid phenomena with complex moving boundaries (gas-liquid, liquid-liquid, gas-liquid solid, etc.) appear everywhere around us (eg rivers, kitchens, bodies, etc.). Conversely, as an actual fluid phenomenon, it can be said that there are few fluid phenomena that do not have complicated moving boundaries.
[0003]
Research on such fluid phenomena is considered important in various fields such as physics, engineering, and medicine. In addition to analysis by experimental methods, analysis by numerical methods is also required. Yes.
[0004]
For example, the design of structures (dams, etc.) for attenuating torrents that cause river erosion and bank breakage is an appropriate fluid simulation method considering the capture of moving boundaries with large deformations and complex moving boundaries. Therefore, there is no choice but to rely on large-scale experiments and empirical rules alone, and the development of effective methods is desired. Similar problems are common to the design and arrangement of jetty and tetrapots on the coast, and the design of pipe shapes to prevent vibration breakage of cooling pipes in industrial equipment. Moreover, if an effective method is obtained, it can be applied to increase the efficiency of inkjet, coating, cleaning, and the like.
[0005]
However, when analyzing fluid phenomena with complex moving boundaries by numerical methods, the handling of boundary conditions is complicated, and the scale of phenomena and the scale of complex moving boundaries must be handled simultaneously. From the above, the analysis of phenomena involves various difficulties.
[0006]
For fluid problems with complex moving boundaries, it is effective to use a fixed grid (Euler grid). This is because, according to the fixed grid, even a complicated boundary shape can be easily expressed, and it is not necessary to generate a grid each time the boundary moves. In addition, since the boundary only moves on a fixed lattice point, the lattice is not twisted and the calculation is not broken, and it is possible to easily cope with a large boundary deformation or a topology change. Because.
[0007]
However, when a fixed grid is used, problems such as numerical diffusion at the boundary and the influence of meshes can naturally occur. Therefore, various devices are required.
[0008]
At present, the level set method is considered to be effective as a method for dealing with such problems (Reference 1 (S. Osher and JA Sethian, J. Comput. Phys. 79, 12 (1988)) and Reference 2 (JA Sethian, Level Set Methods and Fast Marching Methods, Cambridge University Press (1999)). In addition, the level set method has been applied to various problems so far, and has been successful (see the above-mentioned literature 2 and literature 3 (K. Yokoi and F. Xiao, Phys. Rev. E 61 R1016 (2000)). ).
[0009]
The level set method is one of the methods for setting (capturing) the moving boundary. The (N−1) -dimensional boundary is defined as a zero-contour (value of 0) of an N-dimensional function. It is expressed using a certain isoline or isosurface). This N-dimensional function is called a level set function, and its zero contour is called a zero-level set. In this method, the boundary itself is not tracked directly, but the boundary determined as the zero level set is indirectly tracked through a level set function.
[0010]
Here, a function f _ls that satisfies the condition of the following equation (1) is used as the level set function.
[0011]
[Expression 1]

[0012]
Further, the time evolution of the level set function f _ls is calculated according to the advection equation of the following equation (2). In the following equation (2), u represents the velocity of the fluid.
[0013]
[Expression 2]

[0014]
Here, taking the one-dimensional fluid problem as shown in FIG. 9 as an example, the level set function f _ls shown in the above equation (1) is as shown by reference numeral 20, and the boundary 21 between the liquid and the gas 21 , 22 is obtained as a zero _contour surface (zero level set) of the level set function _fls .
[0015]
For this reason, the level set function f _ls shown in the above equation (1) is time-developed according to the advection equation of the above equation (2), and the zero level set of the level set function f _ls at a certain time is obtained to track the boundary. can do.
[0016]
As can be seen from the above equation (1) and FIG. 9, the level set function f _ls is a function having a slope of 1, and the value a of f _ls at a certain point 23 is the distance b from the boundary 21 of the point 23. equal. That is, the level set function f _ls has a property as a distance function representing a distance from the boundary.
[0017]
In such a level set method, since the boundary is not tracked directly but is tracked indirectly using the level set function f _ls , the boundary topology changes (for example, FIG. 10 (a) (b) (As shown in (c), even when two water droplets (liquid) collide with each other), it can be handled naturally without any special treatment. In addition, in order to accurately capture a boundary consisting of curved surfaces on an orthogonal fixed grid, it is necessary to smooth the boundaries to some extent, but in the level set method, the width of the smoothing of the boundaries can be easily controlled. Can be accurately captured. Furthermore, in calculating the curvature κ required for the calculation of surface tension, etc., it is necessary to obtain the divergence of the unit normal vector with respect to the boundary, and the unit normal vector around the boundary is required. By using the level set function f _ls , the curvature κ can be accurately calculated according to the following equation (3).
[0018]
[Equation 3]

[0019]
The CSF (Continum Surface Force) model (Reference 4 (JU Brackbill, DB Kothe and C. Zemach, J. Comput. Phys. 100, 335 (1992)) is referred to based on the curvature κ thus obtained. ) Etc., the surface tension can be easily calculated.
[0020]
[Problems to be solved by the invention]
As described above, the level set method is excellent as a method for setting (capturing) a moving boundary. Generally, the level set function f _ls required for setting a moving boundary is generally expressed by the above equation (2). As the time evolves according to the advection equation, the property of the function (particularly, the property as a distance function) is impaired.
[0021]
That is, even if the level set function f _ls at a certain time satisfies the property of the distance function as shown in FIG. 11 (a), the level set function f after time evolution by the advection equation of the above equation (2). _As shown in FIG. 11B, _ls generally does not satisfy the properties of the distance function.
[0022]
For this reason, in the conventional level set method, it is assumed that the zero level set after time evolution correctly represents the boundary, and the properties of the distance function are expressed with respect to other parts based on this zero level set. Processing for recovery, that is, re-initialization is performed.
[0023]
However, with such re-initialization processing, the zero level set is also slightly shifted, which causes an error regarding the boundary position (FIG. 11 (c)). This error is slight in one calculation step, but the next calculation step is repeated using the level set function (FIG. 11 (d)) in which this error is mixed. It accumulates at every calculation step. For this reason, the position of the boundary greatly deviates from the correct position as time progresses. Accumulation of such errors becomes a serious problem depending on the type of fluid problem, and this error cannot be ignored particularly in the case of a fluid problem having a complicated moving boundary.
[0024]
The present invention has been made in consideration of such points, and effectively prevents accumulation of errors during re-initialization processing, which is a problem in the level set method, so that the moving boundary can be accurately defined on the fixed grid. An object of the present invention is to provide a moving boundary setting method that can be captured and a computer program that realizes the method.
[0025]
[Means for Solving the Problems]
As a first solution, the present invention is used in a fluid simulation for numerically analyzing a fluid phenomenon having a moving boundary, and sets a moving boundary that moves on a fixed grid using a level set method. This is a zero level set tracking function for tracking the zero level set representing the moving boundary as a function independent of the level set function required for setting the moving boundary moving according to the advection equation. Preparing a zero level set tracking function whose zero contour surface is a zero level set and having the same shape as the level set function, and applying the zero level set tracking function to the advection equation according to the advection equation. Time to develop the time by the time interval determined from the stability condition required for stable calculation Based on the expansion step and the zero level set tracking function after time evolution by the time evolution step, a zero level set representing the movement boundary is obtained, and the movement boundary is set based on the obtained zero level set. And a boundary setting step for constructing a level set function required for the method, and the time evolution step and the boundary setting step are repeated until a predetermined end time is reached. To do.
[0026]
In the first solving means described above, it is preferable that the method further includes a step of obtaining a density function used for identifying a substance inside and outside the moving boundary based on the level set function constructed by the boundary setting step. . In the time development step, the advection equation is preferably solved by a CIP method. Further, the fixed grating is preferably an orthogonal fixed grating.
[0027]
As a second solution, the present invention is used in a fluid simulation for numerically analyzing a fluid phenomenon having a moving boundary, and sets a moving boundary that moves on a fixed grid using a level set method. In the computer program that realizes the setting method, the zero level set representing the moving boundary is set as a function independent of the level set function required for setting the moving boundary moving according to the advection equation. A zero level set tracking function intended to be tracked, the zero contour surface of which is a zero level set, and preparing a zero level set tracking function having the same shape as the level set function; To stably calculate the advection equation according to the advection equation with the level set tracking function A zero level set representing a moving boundary based on a time evolution procedure that evolves time by a time interval determined from a required stability condition, and a zero level set tracking function after time evolution by the time evolution procedure. And a boundary setting procedure for constructing a level set function required for setting the moving boundary based on the determined zero level set, and the time evolution procedure and the boundary setting procedure are set to a predetermined value. Provided is a computer program for realizing a moving boundary setting method, wherein the computer repeatedly executes the computer until an end time is reached.
[0028]
According to the present invention, a zero level set tracking function is prepared as a function independent of the level set function f _ls used in the level set method, the zero level set tracking function is time-developed according to the advection equation, and the zero level set tracking function is provided. Based on the above, the zero level set representing the moving boundary is obtained, and the level set function required for setting the moving boundary is constructed based on the obtained zero level set. Errors do not accumulate in the level set function, and it is possible to effectively prevent the accumulation of errors during the reinitialization process, which has been a problem with the conventional level set method.
[0029]
DETAILED DESCRIPTION OF THE INVENTION
Embodiments of the present invention will be described below with reference to the drawings. Here, FIG. 1 is a flowchart for explaining an embodiment of a moving boundary setting method according to the present invention, and FIGS. 2 (a), (b), (c) and (d) are one-dimensional outlines of the present embodiment. It is a conceptual diagram for exemplifying and explaining the moving boundary problem. The moving boundary setting method according to the present embodiment sets the moving boundary moving on the orthogonal fixed grid using the level set method, and is used for analysis of a fluid phenomenon having the moving boundary. In addition to the orthogonal fixed grating, an arbitrary fixed grating such as a triangular grating can be used as the fixed grating.
[0030]
As shown in FIG. 1, in the moving boundary setting method according to the present embodiment, first, a zero level set representing a moving boundary is tracked as a function independent of the level set function f _ls used in the level set method. A zero level set tracking function f _zero is prepared for the purpose only (step 101). In a broad sense, the zero level set tracking function f _zero is also a level set function, and the zero contour surface is a zero level set. This is called a set function, and only the zero contour surface of such a level set function is called a zero level set.
[0031]
Next, the zero level set tracking function f _zero is evolved by the time interval determined from the stability condition according to the advection equation of the following equation (4) (step 102). Specifically, for example, as shown in FIGS. 2A and 2B, the zero level set tracking function f _zero is moved in the right direction in accordance with the movement of the boundary. At this time, the shape of the zero level set tracking function f _zero collapses with time development.
[0032]
[Expression 4]

[0033]
It should be noted that the CIP method (Cubic Interpolated Propagation Method) (Reference 5 (T. Yabe and T. Aoki, Comput. Phys. Commun.)) Is a highly accurate and robust method for solving the advection equation (4). 66, 219 (1991)) and literature 6 (T. Nakamura and T. Yabe , Comput. Phys. Commun. 120, 122 (1999)) is preferably used reference), Thus, the zero level set tracking function f _zero The collapse of the shape can be suppressed. In addition to the CIP method, an arbitrary method such as a PPM (Piecewise Parabolic Method) method can be used for solving the advection equation.
[0034]
Thereafter, a zero level set representing a moving boundary is obtained based on the zero level set tracking function f _zero after time evolution, and the moving boundary is set based on the obtained zero level set. A level set function f _ls is constructed (step 103). Specifically, for example, as shown in FIG. 2 (c), a zero level set representing a moving boundary is obtained based on the zero level surface of the zero level set tracking function f _zero , and the obtained zero level set is obtained. Based on this, a level set function f _ls having the property of a distance function is constructed. The level set function f _ls is used for boundary smoothing, curvature calculation, and the like.
[0035]
Hereinafter, details of the processing in step 103 (processing for constructing the level set function f _ls from the zero level set tracking function f _zero ) will be described.
[0036]
In step 103, first, the zero level surface of the zero level set tracking function f _zero is obtained by interpolating between the lattice points.
[0037]
Next, the zero contour surface obtained in this way is set as a zero level set representing a moving boundary, and a level set function f _ls of lattice points within one cell from the boundary is defined as a fast marching method (fast marching method ( According to Reference 2 and Reference 7 (see D. Adalsteinsson and JA Sethian, J. Comput. Phys. 148, 2 (1999)), the Eikonal equation of the following equation (5) is calculated.
[0038]
[Equation 5]

[0039]
Next, the value of the level set function f _ls obtained in this way is fixed, and a method based on iterative calculation by Sussmann et al. (Reference 8 (M. Sussman, P. Smereka and S. Osher, J. Comput. Phys. 114, 146 (1994))), the level set function f _ls of the part beyond one cell is constructed by performing convergence calculation until the following expression (6) converges. In general, the method of Sussman et al. Has a problem that the zero level set is shifted during the iterative calculation, but here, the values around the zero level set are calculated in advance by the first marching method, and at the time of actual convergence calculation The problem does not occur because it is fixed.
[0040]
[Formula 6]

[0041]
Here, τ is the time for the iterative calculation and need not be the same as the real time. S is a smoothed sign function and is expressed by the following equation (7).
[0042]
[Expression 7]

[0043]
Since the level set function f _ls around the zero level set is fixed, it is possible to set S (f _ls ) = 1 when the stability of calculation does not matter.
[0044]
Here, as the initial value in the convergence calculation by the above equation (6), the value of f _{ls in} the previous calculation step, 1, -1, etc. are not used, but the calculation is performed by the advection equation of the above equation (2). A value may be used, which can greatly reduce the time required for convergence. For such convergence calculation, the above-described CIP method can be used, and a first-order upwind difference method can be used. This is because even if a relatively low accuracy method such as the first-order upwind difference method is used, an error mixed by the method is corrected by the convergence calculation according to the above equation (6). If the calculation time is not concerned, the first marching method described above can be used instead of the convergence calculation according to the above equation (6).
[0045]
Thereafter, based on the level set function f _ls constructed in this way, a density function (color (function) f _color ) used to identify the substance inside and outside the moving boundary is obtained (step 104). The density function f _color can be calculated from the level set function f _ls by using a smoothed Heaviside function H _α according to the following equations (8) and (9). Here, 2α is the width of the boundary smoothing.
[0046]
[Equation 8]

[0047]
Then, the processes in steps 102 to 104 are repeated until a predetermined end time is reached (step 105). When returning to step 102 and continuing the processing, the function that is the basis of the time evolution calculation is the zero level set tracking function f _zero . That is, in step 103, when the level set function f _ls is constructed based on the zero level set tracking function f _zero , an error is mixed as in the conventional level set method (see FIG. 2 (c)). The time evolution calculation in step 102 uses a zero level set tracking function f _{zero in} which no error is mixed (see FIG. 2D).
[0048]
Thus, according to the present embodiment, the zero level set tracking function f _zero is prepared as a function independent of the level set function f _ls used in the level set method, and the zero level set tracking function f _zero is determined according to the advection equation. Based on the zero level set tracking function f _zero , a zero level set representing a moving boundary is obtained, and a level set function f required for setting the moving boundary is obtained based on the obtained zero level set. _{Since ls} is constructed, errors are not accumulated in the level set function f _ls due to the calculation of time evolution, and the problem of error accumulation during the reinitialization process, which has been a problem with the conventional level set method, is solved. It can be effectively prevented.
[0049]
Note that the moving boundary setting method according to the above-described embodiment can be realized as a computer program that runs on a computer system 10 as shown in FIG. Here, the computer system 10 includes a bus 11, a processor 12, a memory 13 and a hard disk 14 connected to the bus 11, and peripheral devices (an input device 15 such as a keyboard and a mouse, a display, a printer, etc.) connected to the bus 11. Output device 16, FD drive 17 and CD-ROM drive 18). The computer program as described above is stored in a computer-readable recording medium such as the flexible disk 19a and the CD-ROM 19b, and is read into the memory 13, the hard disk 14, or the like. In the computer system 10, the above-described procedure can be realized by being sequentially read from the processor 12 and executed.
[0050]
【Example】
Next, specific examples of the above-described embodiment will be described.
[0051]
First, an example in which the above-described embodiment is applied to the Zalesak problem (see Reference 9 (ST Zalesak, J. Comput. Phys. 31, 335 (1979))) will be described.
[0052]
As shown in FIG. 4, the moving boundary problem of this embodiment is to rotate a circular two-dimensional rigid body 31 having a notch in a part in a counterclockwise direction along an arrow 32.
[0053]
Here, the initial condition of the rigid body 31 was set to the following formula (10).
[0054]
[Equation 9]

[0055]
In addition, the lattice for calculation is a 100 × 100 orthogonal fixed lattice, and the rotation speed of the rigid body 31 is calculated in one rotation in 628 calculation steps under the velocity field as in the following equation (11). It was like that. (That is, the time interval Δt is set to 1/628.)
[0056]
[Expression 10]

[0057]
Further, the level set function f _ls is assumed to be as shown in FIG. Note that corresponds to such a level set function f _ls, zero level set (boundary) and density function f _color in the initial state, respectively, is as shown in FIG. 5 (b) (c).
[0058]
Under such conditions, a zero level set tracking function f _zero having the same form as the level set function f _ls is prepared, and this zero level set tracking function f _zero is time-developed by the CIP method according to the above equation (11). Based on the zero level set tracking function f _zero , a zero level set representing a moving boundary is obtained, and a level set function f _ls required for setting the moving boundary is obtained based on the obtained zero level set. Constructed and set moving boundaries.
[0059]
6 (a), 6 (b) and 6 (c), (d) are diagrams showing the state of the boundary after one rotation and ten rotations of the rigid body 31, respectively. FIGS. _6A and 6C are diagrams showing the zero level set (boundary) obtained based on the level set function f _ls together with an exact solution (initial condition). FIGS. It is a figure which shows the density function _fcolor calculated | required based on the level set function _fls .
[0060]
As shown in FIGS. 6 (a), (b), (c), and (d), in this embodiment, it can be seen that the boundary is accurately captured even after 10 revolutions as well as after 10 revolutions of the rigid body 31.
[0061]
Here, as a comparative example, FIGS. 7A and 7B show a boundary state when the rigid body 31 is rotated once by using the above-described conventional level set method. FIG. 7A shows a zero level set (boundary) together with an exact solution (initial condition), and FIG. 7B shows a density function f _color . Note that, in the level set method in this comparative example, the CIP method was used as a solution method of the advection equation when the level set function _fls is developed over time.
[0062]
As shown in FIGS. 7 (a) and 7 (b), it can be seen that in this comparative example, the movement boundary is considerably deviated from the initial state, even after one rotation. Note that this deviation increased with time.
[0063]
Next, an example in which the above-described embodiment is applied to a fluid problem will be described.
[0064]
As shown in FIGS. 8A and 8B, the fluid problem of this embodiment is that the rectangular rigid body 33 is dropped on the water surface 34. The lattice for calculation was a 40 × 40 orthogonal fixed lattice, and the time interval (Δt) was determined so as to satisfy the stability condition. As a fluid simulation method, a CUP method (Reference 10 (T. Yabe and PY Wang, J. Phys. Soc. Jpn. 60, 2105 (1991)) is a calculation method of multiphase flow based on the Poisson equation of pressure. )))) Was used.
[0065]
As a result, as shown in FIG. 8B, it can be seen that in this embodiment, the shape of the rigid body 33 is well captured despite the relatively small number of lattices. In contrast, in the comparative example shown in FIG. 8C (using the conventional level set method), it can be seen that the corners of the rigid body 33 tend to be rounded due to an error during the reinitialization process.
[0066]
【The invention's effect】
As described above, according to the present invention, it is possible to effectively prevent accumulation of errors during re-initialization processing, which is a problem in the level set method, and to accurately capture the moving boundary on the fixed grid. Therefore, the present invention is extremely useful for designing a structure corresponding to a rapid stream of a river, a large-scale wave on a coast, a fluid phenomenon in an industrial apparatus, and the like.
[Brief description of the drawings]
FIG. 1 is a flowchart for explaining an embodiment of a moving boundary setting method according to the present invention.
FIG. 2 is a conceptual diagram for explaining the outline of the embodiment shown in FIG. 1 taking a one-dimensional moving boundary problem as an example;
FIG. 3 is a diagram showing an example of a computer system that executes a computer program that implements the embodiment shown in FIG. 1;
FIG. 4 is a schematic diagram for explaining the outline of the Zalesak problem (two-dimensional rigid body rotation problem) as an embodiment of the present invention.
FIG. 5 is a diagram showing an example of a level set function, a zero level set, and a density function in the embodiment shown in FIG. 4;
6 is a diagram showing a boundary state after rotation (after one rotation and after ten rotations) of the rigid body in the embodiment shown in FIG. 4;
FIG. 7 is a view showing a boundary after a rigid body is rotated (after one rotation) when a conventional level set method is used.
FIG. 8 is a diagram showing a calculation result of a fluid problem (problem of dropping a rectangular rigid body onto the water surface) as an embodiment of the present invention.
FIG. 9 is a diagram illustrating an example of a level set function used in a one-dimensional fluid problem.
FIG. 10 is a diagram illustrating an example of a one-dimensional fluid problem in which the boundary topology changes.
FIG. 11 is a conceptual diagram for explaining an overview of a conventional level set method.
[Explanation of symbols]
f _zero level set tracking function f _ls level set function 10 computer system 20 level set functions 21, 22 boundary

Claims (5)

In a moving boundary setting method, which is used in a fluid simulation that numerically analyzes a fluid phenomenon having a moving boundary and sets a moving boundary that moves on a fixed grid using a level set method,
A zero level set tracking function intended to track the zero level set representing the moving boundary as a function independent of the level set function required for setting the moving boundary moving according to the advection equation, such as zero Providing a zero level set tracking function whose elevation is a zero level set and having the same shape as the level set function;
A time evolution step of evolving the zero level set tracking function according to the advection equation by a time interval determined from a stability condition required to stably calculate the advection equation;
Based on the zero level set tracking function after time evolution by the time evolution step, a zero level set representing a movement boundary is obtained, and the movement boundary is set based on the obtained zero level set. A boundary setting step for constructing a level set function
A method for setting a moving boundary, wherein the time development step and the boundary setting step are repeated until a predetermined end time is reached.   The method according to claim 1, further comprising: determining a density function used for identifying substances inside and outside the moving boundary based on the level set function constructed by the boundary setting step.   3. The method according to claim 1, wherein the advection equation is solved by a CIP method in the time evolution step.   The method according to claim 1, wherein the fixed grating is an orthogonal fixed grating. In a computer program for realizing a moving boundary setting method, which is used in a fluid simulation for numerically analyzing a fluid phenomenon having a moving boundary and sets a moving boundary moving on a fixed grid using a level set method.
As a procedure to be executed on the computer,
A zero level set tracking function intended to track the zero level set representing the moving boundary as a function independent of the level set function required for setting the moving boundary moving according to the advection equation, such as zero A procedure for preparing a zero level set tracking function whose elevation is a zero level set and having the same shape as the level set function;
A time evolution procedure for evolving the zero level set tracking function according to the advection equation by a time interval determined from a stability condition required to stably calculate the advection equation;
Based on the zero level set tracking function after time evolution by the time evolution procedure, a zero level set representing a movement boundary is obtained, and the movement boundary is set based on the obtained zero level set. A boundary setting procedure for constructing a level set function
A computer program for realizing a moving boundary setting method, wherein the computer repeatedly executes the time evolution procedure and the boundary setting procedure until a predetermined end time is reached.

Images