(57)【要約】【課題】  ダンピング関数を規定する物理量を精度よく計算する。【解決手段】  本実施の形態におけるダンピング関数を、コルモゴロフの速度スケールを用いて定義するが、乱流エネルギの全ての波数成分に対する散逸率の代わりに格子フィルタサイズ以下の波数成分に対するSGS乱流エネルギの散逸率εSGSをそのまま用いるのではなく、乱流をモデル化する際の所定波数（格子フィルタサイズ）及び壁面からの距離により、格子フィルタサイズ以下の波数成分に対するSGS乱流エネルギの散逸率εSGSを乱流エネルギの全ての波数成分に対する散逸率に変換する（ステップＡ６）。【選択図】        図３
2005166015.tif 000002
【特許請求の範囲】
【請求項１】
各々異なる複数種類の乱流を規定する波数範囲内の所定波数以上の乱流成分をモデル化して解析する際に用いられるダンピング関数を規定する物理量を、該物理量を規定するパラメータを用いてコンピュータに計算させる物理量計算プログラムであって、
前記パラメータの値として、前記モデル化された領域内において定まる前記パラメータの値に基づいて前記波数範囲の全域において定まる前記パラメータの値に予測された値を用いて、前記物理量を計算することを特徴とする物理量計算プログラム。
【請求項２】
前記予測された値は、前記モデル化された領域内において定まる前記パラメータの値を、前記所定波数及び前記壁面からの距離を用いて、前記波数範囲の全域において定まる前記パラメータの値に変換した値であることを特徴とする請求項１に記載の物理量計算プログラム。
【請求項３】
前記物理量は、前記乱流変動が消散する波数に関連の強い物理量であることを特徴とする請求項１又は請求項２に記載の物理量計算プログラム。
【請求項４】
前記物理量は、コルモゴロフの速度スケールであることを特徴とする請求項１乃至請求項３の何れか１項に記載の物理量計算プログラム。
【請求項５】
前記物量量は、乱流エネルギーの平方根であることを特徴とする請求項１又は請求項２に記載の物理量計算プログラム。
【請求項６】
コンピュータに、請求項１乃至請求項５の何れか１項に記載の物理量計算プログラムにより計算された前記物理量に基づいてダンピング関数を計算させるダンピング関数計算プログラム。
【請求項７】
コンピュータに、請求項１乃至請求項５の何れか１項に記載の物理量計算プログラムにより計算された前記物理量に基づいて乱流粘性を計算させる乱流粘性計算プログラム。
【請求項８】
コンピュータに、請求項１乃至請求項５の何れか１項に記載の物理量計算プログラムにより計算された前記物理量に基づいて前記流体を解析させる流体解析プログラム。
【発明の詳細な説明】
【技術分野】
【０００１】
本発明は、物理量計算プログラム、ダンピング関数計算プログラム、乱流粘性計算プログラム、及び流体解析プログラムにかかり、特に、ダンピング関数を規定する物理量を計算するための物理量計算プログラム、該物理量計算プログラムに関連するダンピング関数計算プログラム、乱流粘性計算プログラム、及び流体解析プログラムに関する。
【背景技術】
【０００２】
乱流の高精度数値解析手法の一つであるLarge Eddy Simulation (LES) では（非特許文献１等）、支配方程式に格子フィルタを施すことにより、格子フィルタサイズより大きなスケールの乱流成分（ＧＳ成分）は直接計算し、格子フィルタサイズより小さなスケールの乱流成分はサブグリッドスケール(SGS)モデルによってモデル化して解く。従って、LESの計算精度は、用いたSGSモデルの精度に大きく依存する。代表的なSGSモデルの多くは渦粘性近似に基づいており、乱流粘性を与えることによってモデル化をする。このとき、物体壁面近傍では壁面の効果を考慮するためのダンピング関数を乗じることが必要とされる。即ち、壁面付近の流体も含めて流体全体を、壁面から遠いところの十分に乱れた流れ場を想定した乱流のモデル化により解析すると、壁面付近の流体については、乱流粘性が大きくなりすぎる。そこで、この壁面近傍での減衰効果を表現するために、ダンピング関数を用いている。よって、このダンピング関数の与え方が計算精度に大きな影響を与えることがよく知られている。
【０００３】
ところで、乱流を解析する手法としては、ＬＥＳの他に、K- εモデルを用いる手法がよく知られる。K- εモデルであっても上記ダンピング関数を必要とする。K- εモデルのダンピング関数として、コルモゴロフの速度スケールをパラメータとして用いると適切なダンピング関数を構成できることが報告されている。
【０００４】
なお、乱流変動は流体のもつ分子粘性によってある波数において消散する性質をもつが、コルモゴロフの速度スケールは、この波数に関連の強い速度スケールであり、乱流エネルギの散逸率εを用いて表される。
【０００５】
ここで、上記コルモゴロフの速度スケールを用いた上記ダンピング関数を、K- εモデルではなく、LESにおいて用いることも考えられる。
【０００６】
しかし、ＬＥＳでは、上記のように、格子フィルタサイズより小さなスケールの乱流成分のみをモデル化して解いており、乱流エネルギの散逸率εを解いていないので、そのまま用いることはできない。
【０００７】
そこで、ＬＥＳにおいて、上記コルモゴロフの速度スケールを用いた上記ダンピング関数を用いる場合に、乱流エネルギの全体の散逸率εに代えて、格子フィルタサイズ以下の波数成分に対するSGS乱流エネルギの散逸率を用いることが提案されている。
【非特許文献１】日本機会学会論文集６４巻623号、Ｂ（1998）、ｐｐ.2208-2215  「ＬＥＳデータに基づくこう配拡散型熱流束モデルの高精度化の検討」
【発明の開示】
【発明が解決しようとする課題】
【０００８】
しかしながら、格子フィルタサイズ以下の波数成分に対するSGS乱流エネルギの散逸率は、乱流エネルギの全ての波数成分に対する散逸率εより小さくかつ格子サイズに依存する。従って、上記格子フィルタサイズ以下の波数成分に対するSGS乱流エネルギの散逸率を用いて上記ダンピング関数を規定すると、適切なダンピング効果が得られないとともに計算に用いる格子サイズに大きく依存したダンピング関数となってしまい、計算精度の低下を招く。
【０００９】
本発明は、上記事実に鑑み成されたもので、ダンピング関数を規定する物理量を精度よく計算することのできる物理量計算プログラム、該物理量計算プログラムに関連するダンピング関数計算プログラム、乱流粘性計算プログラム、及び流体解析プログラムを提供することを目的とする。
【課題を解決するための手段】
【００１０】
上記目的を達成するために請求項１に記載の発明は、各々異なる複数種類の乱流を規定する波数範囲内の所定波数以上の高波数の乱流成分をモデル化して解析する際に用いられるダンピング関数を規定する物理量を、該物理量を規定するパラメータを用いてコンピュータに計算させる物理量計算プログラムであって、前記パラメータの値として、前記モデル化された領域内において定まる前記パラメータの値に基づいて前記波数範囲の全域において定まる前記パラメータの値に予測された値を用いて、前記物理量を計算することを特徴とする。
【００１１】
即ち、本発明の物理量計算プログラムは、各々異なる複数種類の乱流を規定する波数範囲内の所定波数以上の乱流成分をモデル化して解析する際に用いられるダンピング関数を規定する物理量を、該物理量を規定するパラメータを用いてコンピュータに計算させるプログラムである。
【００１２】
上記のように、ダンピング関数は、各々異なる複数種類の乱流を規定する波数範囲内の所定波数以上の乱流成分をモデル化して解析する際に用いられ、より詳細には、物体壁面付近を流れる流体への壁面の効果を定義するものである。
【００１３】
本発明では、このダンピング関数を規定する物理量を、該物理量を規定するパラメータを用いてコンピュータに計算させる。この場合、上記のようにモデル化された領域内において定まるパラメータの値を用いることも考えられるが、このパラメータの値は、当該モデル化された領域内において定まるものであるので、ダンピング関数を規定する物理量を精度よく計算することができない。
【００１４】
そこで、本発明は、上記パラメータの値として、上記モデル化された領域内において定まるパラメータの値に基づいて上記波数範囲の全域において定まるパラメータの値に予測された値を用いて、ダンピング関数を規定する物理量を計算するようにしている。
【００１５】
このように、モデル化された領域内において定まる前記パラメータの値に基づいて前記波数範囲の全域において定まる前記パラメータの値に予測された値を用いてダンピング関数を規定する物理量を計算するので、ダンピング関数を規定する物理量を精度よく計算することができる。
【００１６】
ところで、上記予測された値は、例えば、上記モデル化された領域内において定まる前記パラメータの値を、前記所定波数及び前記壁面からの距離を用いて、波数範囲の全域において定まるパラメータの値に変換した値を用いるようにしてもよい。
【００１７】
また、上記ダンピング関数を規定する物理量としては、乱流変動が消散する波数に関連の強い物理量、即ち、流体のもつ分子粘性によって定まりかつ乱流変動が消散する波数に関連の強い物理量である、例えば、コルモゴロフの速度スケールとしてもよく、また、乱流エネルギーの平方根としていもよい。
【００１８】
以上説明したように、ダンピング関数を規定する物理量を精度よく計算することができる物理量計算プログラムを得ることができるので、このプログラムにより計算される該物理量を用いて他の物理量を計算するプログラムが得られる。
【００１９】
即ち、請求項６記載の発明にかかるダンピング関数計算プログラムは、コンピュータに、請求項１乃至請求項５の何れか１項に記載の物理量計算プログラムにより計算された前記物理量に基づいてダンピング関数を計算させるプログラムである。
【００２０】
また、請求項７記載の発明にかかる乱流粘性計算プログラムは、コンピュータに、請求項１乃至請求項５の何れか１項に記載の物理量計算プログラムにより計算された前記物理量に基づいて乱流粘性を計算させるプログラムである。なお、乱流粘性計算プログラムとしては、コンピュータに、請求項６に記載のダンピング関数計算プログラムにより計算されたダンピング関数に基づいて乱流粘性を計算させるプログラムでもよい。
【００２１】
更に、請求項８に記載の流体解析プログラムは、コンピュータに、請求項１乃至請求項５の何れか１項に記載の物理量計算プログラムにより計算された前記物理量に基づいて前記流体を解析させるプログラムである。なお、流体解析プログラムは、コンピュータに、請求項６記載のダンピング関数計算プログラムにより計算された前記ダンピング関数に基づいて前記流体を解析させるプログラムでもよい、更には、コンピュータに、乱流粘性計算プログラムにより計算された乱流粘性に基づいて前記流体を解析させるプログラムでもよい。
【００２２】
また、次の装置も提案される。
【００２３】
即ち、各々異なる複数種類の乱流を規定する波数範囲内の所定波数以上の乱流成分をモデル化して解析する際に用いられるダンピング関数を規定する物理量を、該物理量を規定するパラメータを用いて計算する物理量計算装置であって、前記パラメータの値として、前記モデル化された領域内において定まる前記パラメータの値に基づいて前記波数範囲の全域において定まる前記パラメータの値に予測された値を用いて、前記物理量を計算することを特徴とする物理量計算装置である。
【００２４】
なお、前記予測された値を、前記モデル化された領域内において定まる前記パラメータの値を、前記所定波数及び前記壁面からの距離を用いて、前記波数範囲の全域において定まる前記パラメータの値に変換した値としてもよい。
【００２５】
また、前記物理量を、前記乱流変動が消散する波数に関連の強い物理量としてもよく、コルモゴロフの速度スケールとしてもよい。
【００２６】
更に、前記物量量を、乱流エネルギーの平方根としてもよい。
【００２７】
また、上記物理量に基づいてダンピング関数を計算するダンピング関数計算装置や、上記物理量に基づいて乱流粘性を計算する乱流粘性計算装置や、上記物理量に基づいて前記流体を解析する流体解析装置が提案される。
【発明の効果】
【００２８】
以上説明したように本発明によれば、モデル化された領域内において定まる前記パラメータの値に基づいて前記波数範囲の全域において定まる前記パラメータの値に予測された値を用いてダンピング関数を規定する物理量を計算するので、ダンピング関数を規定する物理量を精度よく計算することができる、という効果がある。
【発明を実施するための最良の形態】
【００２９】
以下、図面を参照して本発明の実施の形態の一例を詳細に説明する。
【００３０】
最初に、本実施の形態のＬＥＳにおける流体解析の原理を説明する。
【００３１】
本実施の形態におけるダンピング関数を次のように定義する。
【００３２】
【数１】
2005166015.tif 000003
【００３３】
ここで、（２）式のy*は、ダンピング関数を規定する物理量であり、（３）式のｕεは、該物理量を規定するパラメータであり、コルモゴロフの速度スケールである。
【００３４】
また、yは壁面からの最短距離、ν（ヌー）は動粘性係数、ΔはLESにおける格子フィルター幅（= 計算格子幅）である。またＡ、Ｂおよびはモデル定数である。εSGSはＳＧＳの乱流エネルギkSGSの散逸率である。
【００３５】
εSGSの与え方には、εSGSの輸送方程式を解いて与える方法と、kSGSなどから予測する方法がある。後者の場合、例えば次のように与えられる。
【００３６】
【数２】
2005166015.tif 000004
【００３７】
上記式(4)に従ってεSGSを予測する場合にはkSGSを見積もる必要があるが、これにも幾つかの方法がある。代表的なものを以下に示す。
【００３８】
第１に、乱流エネルギkSGSの輸送方程式を解いて与える方法がある。
【００３９】
ここで、代表的な乱流エネルギkSGSの輸送方程式は、
【００４０】
【数３】
2005166015.tif 000005
【００４１】
第２に、計算で得られる瞬時の流れ場に、さらに空間フィルタ
【００４２】
【数４】
2005166015.tif 000006
【００４３】
を施すことによって予測する方法がある。
【００４４】
【数５】
2005166015.tif 000007
【００４５】
上記（３）式において、(y/Δ)成分がある理由を説明する。
【００４６】
上記のように、乱流を解析する手法としては、ＬＥＳの他に、K- εモデルを用いる手法がある。K- εモデルにおいてもコルモゴロフの速度スケールを用いるダンピング関数を用いている。この場合、K- εモデルにおいて用いられるコルモゴロフの速度スケールは、方程式を解いて得られる乱流エネルギ全体に対する散逸率εを用いてｕε＝（νε）1/4と与えられる。
【００４７】
ここで、上記コルモゴロフの速度スケールを用いた上記ダンピング関数を、K- εモデルではなく、LESにおいて用いることも考えられるが、上記K- εモデルにおいて用いられるコルモゴロフの速度スケールでは、図７（Ａ）に示すように、乱流エネルギの全ての波数成分に対する散逸率を用いているので、図７（Ｂ）に示すように、格子フィルタサイズより小さなスケールの乱流成分をモデル化するＬＥＳでは、乱流エネルギの全ての波数成分に対する散逸率を用いることはできない。
【００４８】
そこで、ＬＥＳにおいて、上記コルモゴロフの速度スケールを用いた上記ダンピング関数を用いる場合に、乱流エネルギの全ての波数成分に対する散逸率に代えて、格子フィルタサイズ以下の波数成分に対するSGS乱流エネルギの散逸率を用いることが提案されている。
【００４９】
しかしながら、格子フィルタサイズ以下の波数成分に対するSGS乱流エネルギの散逸率は、乱流エネルギの全ての波数成分に対する散逸率より小さくかつ格子サイズに依存する。従って、上記格子フィルタサイズ以下の波数成分に対するSGS乱流エネルギの散逸率を用いて上記ダンピング関数を規定すると、適切なダンピング効果が得られないとともに計算に用いる格子サイズに大きく依存したダンピング関数となってしまい、計算精度の低下を招く。
【００５０】
そこで、本実施の形態では、コルモゴロフの速度スケールを算出する際に、乱流エネルギの全ての波数成分に対する散逸率の代わりに格子フィルタサイズ以下の波数成分に対するSGS乱流エネルギの散逸率εSGSをそのまま用いるのではなく、乱流をモデル化する際の所定波数（格子フィルタサイズ）及び壁面からの距離により、格子フィルタサイズ以下の波数成分に対するSGS乱流エネルギの散逸率εSGSを乱流エネルギの全ての波数成分に対する散逸率に変換している。より詳細には、長さ比(y/Δ)2を用いたスケール変換を行った上で利用している。なお、yは壁面近傍における乱流（全波数成分を含む）の長さスケールを代表し、Δは乱流のSGS成分の長さスケールを代表する。上記（３）式から理解されるように、この変換は、格子フィルタサイズ以下の波数成分に対するSGS乱流エネルギの散逸率εSGSから、乱流エネルギの全ての波数成分に対する散逸率を予測していることを表している。これにより、計算に用いたダンピング関数は格子サイズに依存せず、適切なダンピング効果を得ることが可能となる。
【００５１】
次に、本実施の形態にかかる流体解析装置１０を説明する。
【００５２】
図１に示すように、本実施の形態にかかる流体解析装置１０は、物理量計算プログラム、ダンピング関数計算プログラム、乱流粘性計算プログラム、及び流体解析プログラムに従っており、パーソナルコンピュータ（パソコン）１２及びプリンタ１４を備えている。
【００５３】
パソコン１４は、ＣＰＵ、ＲＯＭ、ＲＡＭ（メモリ）、及びハードデイスク装置等を備えて構成されている。物理量計算プログラム、ダンピング関数計算プログラム、乱流粘性計算プログラム、及び流体解析プログラムは、ハードデイスク装置内のハードデイスクに記憶されている。
【００５４】
次に、本実施の形態にかかる物理量計算プログラム、ダンピング関数計算プログラム、乱流粘性計算プログラム、及び流体解析プログラムを説明する。
【００５５】
図２には、流体解析プログラムを示すフローチャートが示されている。なお、本実施の形態では、非圧縮性流れ解析の計算手続きを用いたものである。流体解析プログラムにおける流体解析は、以下のステップＳ４以外は周知事項であるので、簡単に説明する。
【００５６】
ステップＳ１では、計算格子や境界条件、計算条件等を、キーボード等を介して入力し、ステップＳ２において、基本物理変数の初期値を設定する。なお、予め用意した各計算格子点上の基本物理変数の値が記述されたデータファイルからデータを読み込むことによって初期値としてもよい。
【００５７】
ステップＳ３において、解析時間tを０として、計算をスタートする。
【００５８】
ステップＳ４において、ＳＧＳ乱流粘性νｔを計算し、メモリ上に記憶しておく。なお、ＳＧＳ乱流粘性νｔの計算内容については後述する。
【００５９】
ステップＳ５では、運動量方程式の対流項、粘性項、ＳＧＳ応力項、圧力項、時間項、外力項などを計算し、メモリ上に記憶する。なお、ＳＧＳ応力項の計算の際には、ステップＳ４の計算により得られたＳＧＳ乱流粘性νｔを用いる。
【００６０】
ステップＳ６では、ステップＳ５において計算された各項の値に基づいて、流体の速度の時間変化の予測値を算出し、速度場の更新を行う。
【００６１】
ステップＳ７では、ステップＳ６において計算された速度の値に基づき、速度場が連続の式を満足するための圧力方程式の各項を算出し、メモリ上に記憶する。
【００６２】
ステップＳ８では、上記ステップＳ７で構成した圧力方程式は連立方程式となるため、ＳＯＲ法や勾配法等の解法により当該連立方程式を解き、圧力場を更新すると共に、更新された圧力の値に基づいて速度場を補正する。
【００６３】
ステップＳ９では、ＳＧＳ応力項のモデル化において、kSGSの輸送方程式といった付加的な方程式の解を用いるか否かを判断する。
【００６４】
付加的な方程式の解を用いる場合には、ステップＳ１０で、kSGSの輸送方程式といった付加的な方程式の各項を計算し、メモリ上に記憶し、ステップＳ１１で、ステップＳ１０で算出した各項の値に基づき、kSGS等の付加的な物理量の時間変化を算出し、算出した時間変化に基づいて、付加的な物理量を更新する。
【００６５】
ステップＳ１２では、必要に応じて途中結果の出力（ハードディスクの記憶、プリントアウトやデイスプレイでの表示）や乱流統計量の算出などを行う。
【００６６】
ステップＳ１３では、解析時間が解析終了時刻ｔmaxを超えたかどうかを判断する。解析終了時刻を越えていない場合は、ステップＳ１４において解析時間を更新し、ステップＳ４に戻って、ステップＳ４からステップＳ１３の処理を繰り返す。解析終了時刻を越えた場合には、ステップＳ１５において最終計算結果を出力（ハードディスクへの記憶、プリントアウトやデイスプレイでの表示）し、終了する。
【００６７】
上記ステップＳ５、６では、以下の空間フィルタ（格子フィルタ）を施した連続の式およびNavier-Stokes方程式（非圧縮性流体解析の場合）を用いる。
【００６８】
【数６】
2005166015.tif 000008
【００６９】
【数７】
2005166015.tif 000009
【００７０】
そして、ステップＳ７、８では、数６で示した２つの式を連立して得られる圧力方程式に対して計算する。
【００７１】
なお、直接計算できないSGS応力τijに対しては、次の渦粘性近似に基づいてモデル化する。
【００７２】
【数８】
2005166015.tif 000010
【００７３】
この場合、乱流粘性νtを与えること(ステップＳ４)によってSGS応力τijがモデル化される。
【００７４】
次に、ステップＳ４における物理量計算プログラム及び乱流粘性計算プログラムを説明する。図３には、物理量計算プログラムを含む乱流粘性計算プログラムを示すフローチャートが示されている。
【００７５】
ステップＡ１では、用いるＳＧＳモデルがＳＧＳ乱流エネルギkSGSの輸送方程式を用いるものであるかどうかを判断する。
【００７６】
ＳＧＳ乱流エネルギkSGSの輸送方程式を用いる場合には、ステップＡ２において、解かれた方程式の解として、ＳＧＳ乱流エネルギkSGSが与えられる。ＳＧＳ乱流エネルギkSGSの輸送方程式を用いない場合には、ステップＡ３において、計算結果として得られる速度のＧＳ成分（格子フィルタサイズより大きなスケールの乱流成分）にフィルタリング操作を施すことにより、ＧＳ成分と比べてより大きな空間フィルタ（テストフィルタ）がかかった速度場
【００７７】
【数９】
2005166015.tif 000011
【００７８】
を算出する。
【００７９】
ステップＡ４では、ステップＡ３で求めたデストフィルタがかかった速度成分とＧＳ成分の差で与えられる、格子フィルタサイズとテストフィルタサイズの間のスケールの乱流成分を用い、その各空間方向成分の２乗平均値からＳＧＳ乱流エネルギkSGSを以下のように見積もる。
【００８０】
【数１０】
2005166015.tif 000012
【００８１】
ステップＡ５では、ステップＡ２又はステップＡ４で求められたＳＧＳ乱流エネルギkSGSの値を用いて、ＳＧＳ乱流エネルギの散逸率εSGSを以下のように見積もる。
【００８２】
【数１１】
2005166015.tif 000013
【００８３】
なお、ＳＧＳ乱流エネルギの散逸率εSGSの輸送方程式も付加方程式として解くタイプのＳＧＳモデルを用いる場合には、その解として得られるＳＧＳ乱流エネルギの散逸率εSGSを用いる場合も考えられる。
【００８４】
ステップＡ６では、上記（３）式に基づくコスモゴロフ速度スケールを計算し、ステップＡ７では、ステップＡ６で計算されたコスモゴロフ速度スケールを用いて、上記（１）式に基づくダンピング関数を計算する。
【００８５】
以上のようにダンピング関数が計算されると、ＳＧＳ乱流粘性νtを算出するための残りの手続きは、予め選択されたＳＧＳモデルによって異なる。
【００８６】
ksgs−方程式型モデルが選択された場合には、ステップＡ２又はステップＡ４で既に求められたＳＧＳ乱流エネルギkSGSの値を用いてＳＧＳ乱流粘性νtを次のように算出する。
【００８７】
【数１２】
2005166015.tif 000014
【００８８】
SmagorinskyモデルやＭＴＳモデルが選択された場合には、ステップＡ９で、ひずみ速度テンソルのＧＳ成分を以下のように求め、
【００８９】
【数１３】
2005166015.tif 000015
【００９０】
ステップＡ１０においてひずみ速度テンソルのＧＳ成分の大きさを次のように求める。
【００９１】
【数１４】
2005166015.tif 000016
【００９２】
Smagorinskyモデルが選択された場合には、ステップＡ１２では、ステップＡ１０において求められたひずみ速度テンソルのＧＳ成分の大きさを用いて、次のようにＳＧＳ乱流粘性νtを計算する。
【００９３】
【数１５】
2005166015.tif 000017
【００９４】
また、ＭＴＳモデルが選択された場合には、ステップＡ１１において計算された時間スケールＴSと、ステップＡ１０において求められたＧＳ成分のひずみ速度テンソルの大きさを用いて、次のようにＳＧＳ乱流粘性νtを計算する。
【００９５】
【数１６】
2005166015.tif 000018
【００９６】
次に、具体例として、上記渦粘性SGSモデル（Ksgs-方程式型モデル）を用いた場合の計算結果を示す。
【００９７】
(1)式におけるモデル係数Ｂ、Ｃは、壁面漸近挙動（渦粘性近似の定義から、乱流粘性（νt）は壁面上でｙの3乗の関数として表現される）を考慮すると、Ksgs-方程式型モデルの場合Ｂ×Ｃ＝４／３を満たす必要がある。
【００９８】
ここでは、Ｂ＝４／３、Ｃ＝１とする。他の各モデル係数についても、相互関係があり任意性を伴うが、一例として本実施の形態で採用したモデル係数を以下に記す。
【００９９】
(5)式を解くことによってksgsを与える場合
【０１００】
【数１７】
2005166015.tif 000019
【０１０１】
以下、これをMethod Ａとする。
【０１０２】
(6)式によってksgsを与える場合
【０１０３】
【数１８】
2005166015.tif 000020
【０１０４】
以下、これをMethod Ｂとする。
【０１０５】
また、比較として、従来技術における、ＬＥＳにおいて、格子フィルタサイズ以下の波数成分に対するSGS乱流エネルギの散逸率をそのまま用いた場合の実験結果を示す。なお、
【０１０６】
【数１９】
2005166015.tif 000021
【０１０７】
をAbe１モデルとし、
【０１０８】
【数２０】
2005166015.tif 000022
【０１０９】
をAbe２モデルとする。
【０１１０】
計算結果は、完全発達したチャネル乱流（図４参照）において、図５および以下の表１に記すように格子解像度（表１中Δｘ+Δy+Δz+が乱流解析における実質的な格子解像度であり、これらの値が小さいほど格子解像度が高い）をCase1からCase5まで５段階変化させた場合のダンピング関数の変化を図６に示す。
【０１１１】
【表１】
2005166015.tif 000023
【０１１２】
図６から分かるように、図６（Ａ）に示すように、Abe1モデルでは格子解像度に対する依存性が非常に高い。図中に実線で記すダンピング関数はチャネル乱流における最適なダンピング関数の一つであるが、この分布とは大きな相違が見られる。特に最も格子解像度が高いCase1では、全領域でほぼダンピング関数f=0となりダンピング関数の影響が壁面から離れたところまで極端に広がりすぎている。図６（Ｂ）に示すように、Abe2モデルでは、ダンピング関数の影響が壁面から遠くまで及ばないようにモデル係数を変更しているが、改善は見られるものの本質的に格子解像度に依存する性質は不変である。
これに対して、本実施の形態におけるMethod Aでは、図６（Ｃ）に示すように、いずれの格子解像度の場合でも最適なダンピング関数の分布にほぼ一致する結果となり、(3)式でのコスモゴロフの速度スケールの算出において上記（ｙ／Δ）1/2を乗じることにより、格子依存性が解決され適切なダンピング効果が得られることが分かる。また、Method Bでも図６（Ｄ）に示すように、ほぼ同様の結果が得られており、本ダンピング関数は乱流エネルギkSGSの与え方にも特に依存せず、良好な結果を与えることがわかる。
以上説明した実施の形態では、ダンピング関数を（１）式により定義しているが、本発明はこれに限定されるものではなく、次のように定義するようにしてもよい。
【０１１３】
まず、RANSで提案されているダンピング関数は他にもある。例えば、W。 Rodiら（W。 Rodi、 N。 N。 Mansour and V。 Michelassi、 One-Equation Near-Wall Turbulence Modeling with the Aid of Direct Simulation Data、 Journal of Fluids Engineering、 vol。 115、 pp。 196-205、 1993）は次のダンピング関数を提案している。
【０１１４】
【数２１】
2005166015.tif 000024
【０１１５】
この場合は、乱流エネルギkの平方根がコルモゴロフの速度スケールの代わりに用いられている。Dejoanら(A。 Dejoan and R。 Schiestel、LES of unsteady turbulence via a one-equation subgrid-scale transport model、 International Journal of Heat and Fluid Flow、 23、 pp。398-412、 2002 )はLES計算においてこれをそのままの形式で用いた。つまり全波数領域に対する乱流エネルギkの代わりにＳＧＳの乱流エネルギkSGSをそのまま用いた。しかし、これは上記従来技術のダンピング関数の場合と同様に、計算に用いる格子サイズに依存した不適切なダンピング効果を与えることになる。この場合は、次式のように変更すればLESで用いることができると考えられる。
【０１１６】
【数２２】
2005166015.tif 000025
【０１１７】
この場合も、ＳＧＳの乱流エネルギkSGSを長さ比（ｙ/Δ）2を用いてスケール変換することによって全波数領域における乱流エネルギkを予測している。数２２の式が適切なダンピング関数となる理由は、以下のようである。
【０１１８】
式(1)〜(3)で記述される本実施の形態において、乱流エネルギの散逸率εSGSを式(4)を用いて乱流エネルギkSGSから与える場合を考える。式(4)の第２項は壁面のごく近傍（壁座標y+≦５）を除くと第１項に対して無視できるほど小さいことは良く知られており、結局y*を次式で見積もっていることになる。
【０１１９】
【数２３】
2005166015.tif 000026
【０１２０】
数２２の式と数２３の式の違いは係数とべき乗の違いだけであり、数２１の式におけるパラメタＡ、Ｂの与え方によって式(1)〜(3)に一致させることができる。例えば、本実施の形態の具体例で用いているＡ＝１０、Ｂ＝４／３に合わせる場合、数２１の式において
【０１２１】
【数２４】
2005166015.tif 000027
【０１２２】
と与えればよい。
【０１２３】
壁面漸近挙動とは、乱流粘性νｔが壁面上でｙの３乗で立ち上がることを意味している。したがって、ダンピング関数ｆの壁面漸近挙動を考えると、用いるＳＧＳモデルに依存する。ｋSGS−方程式型モデルの場合はｆはｙの２乗、Ｓｍａｇｏｒｉｎｓｋｙモデルの場合は１．５乗、ＭＴＳモデルの場合は１乗となる。Ｓｍａｇｏｒｉｎｓｋｙモデルの場合は例えばＢ＝４／３、Ｃ＝３／４あるいはＢ＝１、Ｃ＝１の組み合わせが考えられる。ダンピング関数を構成するパラメータはこのように用いるＳＧＳモデルに対する依存性を持っている。
【０１２４】
以上説明した例では、渦粘性モデルのうち線形モデルに類するＬＥＳモデルを例にとり説明したが、本発明はこれに限定されず、例えば、非線形モデルにも同様に適用することができる。
【実施例】
【０１２５】
以上説明した種々の例は、流体を解析するものであり、流体として、気体や液体を用いてもよい。この場合、液体として例えばガソリンやオイル等を用いてもよい。例えば、管内をガソリンやオイル等の液体が流れる様子を解析する場合に適用できる。管の内面は壁面となる。また、気体としては空気を用いてもよい。例えば、車を配置し、当該車に送風機により空気を送る場合の当該空気の様子を解析する場合にも適用できる。車の外形面が壁面となる。なお、以上は車を例に挙げて説明したものに過ぎず、本発明は、車関係に限定されずに、流体を解析する他の種々の例に適用可能である。
【図面の簡単な説明】
【０１２６】
【図１】流体解析装置の概観図である。
【図２】流体解析流体解析プログラムを示すフローチャートである。
【図３】流体解析流体解析プログラムにおけるステップＳ４の乱流粘性計算プログラムを示すフローチャートである。
【図４】完全発達したチャネル乱流を示す図である。
【図５】格子解像度を示す概念図である。
【図６】ダンピング関数の実験結果であり、（Ａ）及び（Ｂ）は、従来技術におけるダンピング関数の実験結果であり、（Ｃ）及び（Ｄ）は、本実施の形態におけるダンピング関数の実験結果である。
【図７】乱流を規定する波数と乱流変動のエネルギースペクトルとの関係を示す図であり、（Ａ）は、K- εモデルの場合を、（Ｂ）は、ＬＥＳの場合を、それぞれ示している。
【符号の説明】
【０１２７】
１０    流体解析装置
１２    パソコン
１４    プリンタ
【図１】
2005166015.tif 000028
【図２】
2005166015.tif 000029
【図３】
2005166015.tif 000030
【図４】
2005166015.tif 000031
【図５】
2005166015.tif 000032
【図６】
2005166015.tif 000033
【図７】
2005166015.tif 000034