JP7299613B2

JP7299613B2 - システム、方法およびプログラム

Info

Publication number: JP7299613B2
Application number: JP2019153033A
Authority: JP
Inventors: 正裕岩永
Original assignee: Individual
Current assignee: Individual
Priority date: 2019-08-23
Filing date: 2019-08-23
Publication date: 2023-06-28
Anticipated expiration: 2039-08-23
Also published as: JP2021032694A

Description

特許法第３０条第２項適用平成３０年９月３日、一般社団法人日本機械学会２０１８年度年次大会のＤＶＤ講演論文集事前公開サイトにて発表。該当アドレスｈｔｔｐｓ：／／ｗｗｗ．ｊｓｍｅ．ｏｒ．ｊｐ／ｃｏｎｆｅｒｅｎｃｅ／ｎｅｎｊｉ２０１８／平成３０年９月１１日、一般社団法人日本機械学会２０１８年度年次大会にて発表。

本発明は、ナビエ－ストークス方程式の解析において適切に解を収束させるシステム、方法およびプログラムに関する。

流体の運動を記述するナビエ－ストークス（Navier-Stokes）方程式は、複雑な偏微分方程式で記述され、解析には多大な労力を要する。そこで、数値シミュレーションによって流体の挙動を近似的に予測するＣＦＤ（Computational Fluid Dynamics、数値流体力学）技術が開発されている。特に、非圧縮性流体の流れを差分法で解析する場合において、連続の式を厳密に成立させるために、種々の手法が開発されてきた。

例えば、特開平３－２２９１５６号公報（特許文献１）では、ナビエ－ストークス方程式と連続の式に対して、同じ長さのブロック長に保った連立離散化方式とそれを計数とする連立一次方程式を不完全三角分解付き共役勾配系の解法で数値解を求める技術が開示されている。

しかしながら、特許文献１などの従来技術では、解が発散したり、収束に時間が掛かったりし、実用の観点から充分ではなかった。そこで、ナビエ－ストークス方程式の解を安定的に収束させるさらなる技術が求められていた。

本発明は、上記従来技術における課題に鑑みてなされたものであり、ナビエ－ストークス方程式の解析において適切に解を収束させるシステム、方法およびプログラムを提供することを目的とする。

すなわち、本発明によれば、
ナビエ－ストークス方程式を解析するシステムであって、
解析のモデルの圧力分布および速度分布の初期条件を定義する定義手段と、
ある時刻の速度分布、圧力分布および付加項を用いて、前記ある時刻の次の時刻の速度の湧き出しが０に収束するように、圧力分布を算出する第１の算出手段と、
前記ある時刻の速度分布と前記第１の算出手段によって算出された圧力分布に基づいて、前記ある時刻の次の時刻の速度分布を算出する第２の算出手段と
を含み、
前記第１の算出手段および前記第２の算出手段は、所定回数に達するまで算出を繰り返すことを特徴とする、システムが提供される。

本発明によれば、ナビエ－ストークス方程式の解析において適切に解を収束させるシステム、方法およびプログラムが提供できる。

本実施形態の情報処理装置を構成する機能ブロック図。本実施形態における解法の手順を示すフローチャート。 κの値に応じた∇ｕの収束の様子を示すグラフ。 κ＝０の場合の∇ｕの変化を示すグラフ。 κ＝３５０の場合の∇ｕの変化を示すグラフ。 κの値による∇ｕの変化を示すグラフ。レイノルズ数を変化させて計算した球体の抗力係数を示すグラフ。解析結果と風洞実験結果とを比較した図。

以下、本発明を、実施形態をもって説明するが、本発明は後述する実施形態に限定されるものではない。また、明細書中に表示されている数式は、ベクトルの符号などといった、いわゆる環境依存文字を含むことから、一部の符号などを省略して記載している箇所がある点に留意されたい（イメージで表示されている数式を除く）。

ナビエ－ストークスの運動方程式は、下記の式（１）～（３）のように記述される。

なお、上記式は、代表長さＬおよび代表速度Ｕで、下記式（４）のように無次元化している。なお、式（４）中のｐは圧力分布を、ρは流体の密度を、Ｒｅはレイノルズ数を、μは流体の粘性係数をそれぞれ示している。

上記式（１）～（３）について、下記式（５）の演算を行うと、下記式（６）が得られる。

本実施形態では、ポテンシャル流れを初期条件として、陽解法によって次のステップの速度場を求める。ポテンシャル流れは、連続の式を満たし、滑り無しの条件以外の境界条件を満足している。

上記式（６）において、∇ｕは、速度の湧き出しを示している。速度の湧き出しの時間変化がゼロ、すなわち次のステップで連続の式を維持するためには、下記式（７）を満たさなければならない。

そして、上記式（７）を満たすためには、圧力分布は、下記式（８）を満足する必要がある。

ここで、∇ｕがゼロでない場合であっても、速やかに∇ｕがゼロに収束するために、上記式（８）に付加項を追加し、下記式（９）のようにする。すなわち、式（９）のように、速度の湧き出しである∇ｕに正の定数κを乗じた付加項を右辺に加えることで、各場所の∇ｕの値を小さく保ち、速度の湧き出しの総量をゼロに収束させることができる。なお、付加項は、式（９）に示す∇ｕに定数κを乗じたものに限定されず、後述するように、∇ｕの種々の関数を付加項とすることができる。

ここで、式（６）に上記式（９）を代入すると、下記式（１０）を得る。

上記式（１０）に示すように、付加項によって、速度の湧き出しが存在する場合、すなわち∇ｕが発生する場合であっても、すぐに消滅し、近似的に連続の式を満足する解を得ることができる。

上記のような解法の実行は、コンピュータなどの情報処理装置や、種々の計算システムによって行われる。ここでは、情報処理装置が実行する場合を例に説明する。図１は、本実施形態の情報処理装置１００を構成する機能ブロック図である。

本実施形態の情報処理装置１００は、初期条件定義部１０１、圧力分布算出部１０２、速度分布算出部１０３を含む。

初期条件定義部１０１は、解析の対象となるモデルの圧力分布および速度分布の初期条件をポテンシャル流れによって定義する手段である。また、初期条件定義部１０１は、境界条件などといった種々の条件を定義することができる。

圧力分布算出部１０２は、ある時刻の速度分布、圧力分布および付加項を用いて、次の時刻の速度の湧き出しが０に収束するように、圧力分布を算出する手段である。圧力分布算出部１０２は、解を適切に収束させるために、式（９）のように付加項を加算したうえで、圧力分布を算出する。

速度分布算出部１０３は、ある時刻の速度分布と圧力分布算出部１０２が算出した圧力分布に基づいて、次の時刻の速度分布を算出する手段である。

圧力分布算出部１０２および速度分布算出部１０３における詳細な処理については、後述する。

なお、上述したソフトウェアブロックは、ＣＰＵが本実施形態のプログラムを実行することで、各ハードウェアを機能させることにより、実現される機能手段に相当する。また、各実施形態に示した機能手段は、全部がソフトウェア的に実現されても良いし、その一部または全部を同等の機能を提供するハードウェアとして実装することもできる。

さらに、上述した各機能手段は、必ずしも全てが図１に示すような構成で含まれていなくてもよい。例えば、他の好ましい実施形態では、各機能手段は、複数の情報処理装置１００の協働によって実現されてもよい。

次に、本実施形態の具体的な解法について説明する。図２は、本実施形態における解法の手順を示すフローチャートである。なお、図２に示す処理は、図１に示した情報処理装置１００や、同様の構成を備える計算システムが実行することができる。まず、Ｓ１０００から処理を開始する。

Ｓ１００１では、ポテンシャル流れの速度分布および圧力分布の初期値を定義する。解法の最初の速度分布および圧力分布は、滑り無しの条件以外の境界条件を満足するポテンシャル流れから出発する。ポテンシャル流れでは全領域で連続の式が満足されている。

Ｓ１００２では、ある時刻ｔの速度分布、圧力分布と付加項を用いて、次の時刻ｔ＋ｄｔの速度の湧き出しが０に収束するように、圧力分布を算出する。付加項を追加することにより、例え∇ｕが発生しても、数ステップ後（例えばｔ＋５～６ｄｔ）には∇ｕの値を０に収束させることができる。具体的には、Ｓ１００２では、圧力分布に付加項を追加し、下記式（１１）のようにする。なお、式（１１）は、式（９）の右辺に相当する。

ここで、上記式（１１）中の付加項ｆ（∇ｕ）は、∇ｕの関数であり、発生した∇ｕを自然に数ステップのうちに消滅させる効果をもつ関数である。すなわち、ｆ（∇ｕ）を付加項として追加することで、容易に収束解が得られる。付加項ｆ（∇ｕ）の例としては、式（９）に示したκ∇ｕが代表的であるが、これ以外の関数であってもよい。したがって、例えば、下記式（１２）のような関数を付加項として追加してもよい。なお、下記式（１２）中のλ、αは、実定数である。すなわち、本解法は厳密に∇ｕ＝０を満たすように解くのではなく、多少の∇ｕの存在を許しながら、付加項の効果により発生した∇ｕを消滅させることにより、近似的に連続の式∇ｕ＝０を満足させるものである。

ここで、隣のセルとの距離を下記式（１３）のようにおくと、下記式（１４）を得る。なお、式（１３）中のｄｌは、立方体メッシュの一辺の長さを示す。

次に、Ｓ１００３では、下記式（１５）の演算を繰り返し、圧力pの収束値を求める。

その後、Ｓ１００４では、Ｓ１００３で求めた圧力ｐの収束値を適用して、次のステップの速度分布を求める。次のステップの速度分布は、下記式（１６）より、次のステップの速度をｕ’とし、時間ステップをｄｔとすると、下記式（１７）から求めることができる。

上記式（１７）を以て求めた速度分布は、式（９）を満足する圧力分布を用いたものであることから、連続の式（∇ｕ’＝０）をほぼ満足する。仮に∇ｕ’が値を持った場合であっても、付加項の作用により、演算を数ステップ繰り返すことで、ゼロに減衰する。すなわち、Ｔ＝ｔで∇ｕが発生しても、数ステップ後（例えば、５～６ステップ後）には、ゼロに減衰する。そこで、Ｓ１００５では、速度分布を算出する演算を所定回数実行したか否かによって処理を分岐する。所定回数の演算をした場合には（ＹＥＳ）、Ｓ１００６にて終了する。一方で、所定回数の演算をしていない場合には（ＮＯ）、Ｓ１００２に戻り、さらに次のステップの速度分布を求め、所定回数に達するまで、上記の処理を繰り返す。

次に、本実施形態の解法を適用して解析した結果について説明する。式（１８）は式（１０）の∇ｕをｙに置き換えたものであり、この微分方程式の特性を調べた。図３は、κの値に応じたｙの収束の様子を示すグラフである。図３は、下記式（１８）をｔ＝０においてｙ＝１、κを１００，２００，５００，９５０として、ｄｔ＝０．００１ステップで解析した結果を示している。

図３から、ｙが値を有していても、迅速にゼロに収束していることが示された。式（１０）も、付加項κ∇ｕの効果により、∇ｕが迅速にゼロに収束する。

次に、κの値と∇ｕとの関係を解析した結果について説明する。レイノルズ数Ｒｅを６００として、球体の周りの流れ模様を、κの値を変えて計算し、∇ｕの時間変化を調べた。計算は、代表長さを球体の直径とし、代表速度を入口速度として無次元化し、時間の単位は、球体の直径分の距離を入口速度で以て通過する時間としている。時間の刻みは、０．００１、すなわち球体の直径の１／１０００の距離を入口速度で以て通過する時間間隔として計算した。計算領域は、断面が６×６、流れ方向１２の直方体で、セルは、一辺が０．０５の立方体とし、セルの総数は、３４５６０００である。差分はすべて最も単純で簡単な中心差分法を用いている。

図４は、κ＝０の場合の∇ｕの時間変化を示すグラフである。また、図５は、κ＝３５０の場合の∇ｕの時間変化を示すグラフである。図４および図５の横軸は、無次元時間を示している。また、図４（ａ）および図５（ａ）の縦軸は、計算領域全体の∇ｕの平均値を示し、図４（ｂ）および図５（ｂ）の縦軸は、∇ｕの絶対値の平均値を示している。

κ＝０の場合には、図４に示すように、解は、Ｔ＝０．２７で発散した。一方で、κ＝３５０の場合には、図５に示すように、付加項κ∇ｕの効果により、∇ｕの値はゼロに収束し、また、∇ｕの絶対値も無次元時間が１０を超えると、０．００１５と小さい値となった。

図６は、κの値による∇ｕの変化を示すグラフである。なお、図６における解析の条件は、図４や図５のものと同様である。図６（ａ）は、Ｒｅ＝６００で球体の周りの流れを解析したときの流れ場全体の∇ｕの空間平均値をＴ＝２０～２５で時間平均した値を示している。また、図６（ｂ）は、Ｒｅ＝６００で球体の周りの流れを解析したときの流れ場全体の∇ｕの絶対値の空間平均値をＴ＝２０～２５で時間平均した値を示している。図６（ａ）、（ｂ）より、κの値が大きくなることで、∇ｕやその絶対値がゼロに近づくことが示された。

次に、球体の抗力係数の解析について説明する。図７は、レイノルズ数を変化させて計算した球体の抗力係数を示すグラフである。図７の横軸はレイノルズ数Ｒｅを、縦軸は抗力係数ＣＤを示している。また、図７では、κ＝３５０として計算を行っている。図７の丸で示されるプロットは、断面４×４、流れ方向９の直方体をコントロールボリュームとして計算した抗力係数である。なお、各セルの時間１／１０００あたりの速度変化による非定常項を加味して計算している。また、図７の四角形で示されるプロットは、流れに垂直な球体の表面を構成するセルに加わる応力と、流れに平行な球体の表面を構成するセルに加わる応力との和から求めた抗力係数である。図７の実線は、従来の実験結果を示している。図７より、各プロットとも従来の実験値によく一致していることが示された。

次に、解析結果と風洞実験との比較について説明する。図８は、解析結果と風洞実験結果とを比較した図である。図８（ａ）は、Ｒｅが６００の球体の周りの時刻Ｔ＝５０における流れ模様を示す解析結果である。また、図８（ｂ）は、鉛直円形可視化風洞を用いて流速０．４ｍ／ｓにおいて、直径２４．６ｍｍの球体の周りの流れを煙で以て可視化した図である。実験のＲｅ数は計算条件とほぼ同じ６４０である。図８（ａ）、（ｂ）を比較すると、両者の流れの様子はよく一致しており、本実施形態の解析法によって実際の流れをシミュレーションできることが示された。

以上、説明した本発明の実施形態によれば、ナビエ－ストークス方程式の解析において適切に解を収束させるシステム、方法およびプログラムを提供することができる。

上述した本発明の実施形態の各機能は、Ｃ、Ｃ＋＋、Ｃ＃、Ｊａｖａ（登録商標）等で記述された装置実行可能なプログラムにより実現でき、本実施形態のプログラムは、ハードディスク装置、ＣＤ－ＲＯＭ、ＭＯ、ＤＶＤ、フレキシブルディスク、ＥＥＰＲＯＭ（登録商標）、ＥＰＲＯＭ等の装置可読な記録媒体に格納して頒布することができ、また他装置が可能な形式でネットワークを介して伝送することができる。

以上、本発明について実施形態をもって説明してきたが、本発明は上述した実施形態に限定されるものではなく、当業者が推考しうる実施態様の範囲内において、本発明の作用・効果を奏する限り、本発明の範囲に含まれるものである。

１００…情報処理装置、
１０１…初期条件定義部、
１０２…圧力分布算出部、
１０３…速度分布算出部

特開平３－２２９１５６号公報

Claims

ナビエ－ストークス方程式を解析するシステムであって、
解析のモデルの圧力分布および速度分布の初期条件をポテンシャル流れによって定義する定義手段と、
ある時刻の速度分布、圧力分布および付加項を用いて、前記ある時刻の次の時刻の速度の湧き出しが０に収束するように、圧力分布を算出する第１の算出手段と、
前記ある時刻の速度分布と前記第１の算出手段によって算出された圧力分布に基づいて、前記ある時刻の次の時刻の速度分布を算出する第２の算出手段と
を含み、
前記第１の算出手段および前記第２の算出手段は、所定回数に達するまで算出を繰り返すことを特徴とする、システム。
前記付加項は、速度の湧き出しである∇ｕの関数である、請求項１に記載のシステム。
ナビエ－ストークス方程式を解析する方法であって、
解析のモデルの圧力分布および速度分布の初期条件を定義するステップと、
ある時刻の速度分布、圧力分布および付加項を用いて、前記ある時刻の次の時刻の速度の湧き出しが０に収束するように、圧力分布を算出する第１の算出ステップと、
前記ある時刻の速度分布と前記第１の算出ステップにおいて算出された圧力分布に基づいて、前記ある時刻の次の時刻の速度分布を算出する第２の算出ステップと
を含み、
前記第１の算出ステップおよび前記第２の算出ステップを、所定回数に達するまで繰り返すことを特徴とする、方法。
前記付加項は、速度の湧き出しである∇ｕの関数である、請求項３に記載の方法。
情報処理装置に、請求項３または４に記載の方法の各ステップを実行させるためのプログラム。