(57)【要約】　　　（修正有）【課題】物体形状を表面データの入力のみで演算を行うと共に、高速解法および高精度解法を利用し、効率よく非定常な流体の挙動および物体に作用する力を見積もること。【解決手段】解析空間の境界を設定する座標位置データと流体特性を設定する特性データを入力データとし、境界条件の下に各境界における流体の吹出し及び循環に対する影響係数マトリクスを算出する機能と、影響係数マトリクスを含む連立方程式から前記吹き出し量及び循環量を算出する機能と、吹き出し量及び循環量から求まる速度場から、各境界近傍に渦層を近似して渦度を算出して離散渦要素を求める機能と、離散渦要素の移動速度を算出する機能と、離散渦要素を代表する渦度の時間変化を渦度輸送方程式に基づいて算出する機能を備え、境界から発生した各離散渦要素の移動及び代表する渦度の時間変化をラグランジュ的に追従して流体現象をシミュレートする。【選択図】図７
2004126925.tif 000002
【特許請求の範囲】
【請求項１】
解析空間の境界を設定する座標位置データと流体特性を設定する特性データを入力データとし、
物体表面の相対速度成分に基づく境界条件の下にビオ・サバールの式から、各境界における流体の吹き出し量及び循環量に対する影響係数マトリクスを算出する機能と、
前記影響係数マトリクスを含む連立方程式から前記吹き出し量及び循環量を算出する機能と、
前記吹き出し量及び循環量から求まる速度場から、各境界近傍に渦層を近似して渦度を算出して離散渦要素を求める機能と、
前記離散渦要素の移動速度をビオ・サバールの式から算出する機能と、
前記離散渦要素の代表する渦度の時間変化を渦度輸送方程式に基づいて算出する機能を備え、
境界から発生した各離散渦要素の移動及び代表する渦度の時間変化をラグランジュ的に追従することにより流体現象をシミュレートすることを特徴とする、流体シミュレーションプログラム。
【請求項２】
前記速度場と圧力ポアソン方程式から任意の位置の圧力を算出する機能を備えることを特徴とする、請求項１に記載の流体シミュレーションプログラム。
【請求項３】
物体に作用する力から座標位置データを変更する機能を備えることを特徴とする、請求項２に記載の流体シミュレーションプログラム。
【請求項４】
座標位置データは時系列データであることを特徴とする、請求項１乃至３のいずれかに記載の流体シミュレーションプログラム。
【請求項５】
前記離散渦要素の移動速度の算出において、各離散渦要素間の相互作用が対称であることを利用し、任意の二つの離散渦要素において、互いの離散渦要素に誘起しあう速度を同時に求めることを特徴とする、請求項１に記載の流体シミュレーションプログラム。
【請求項６】
前記離散渦要素の移動及び代表する渦度の変化を算出するにあたり、その演算頻度を自動的に最適化することを特徴とする、請求項１に記載の流体シミュレーションプログラム。
【請求項７】
前記渦度の生成モデルにおいて、渦層の近似高さ、パネル分割長さ、および時間刻みを自動的に最適化することを特徴とする、請求項１に記載の流体シミュレーションプログラム。
【請求項８】
前記渦層の近似高さとパネル分割長さは、摩擦速度に基づいて自動的に最適化することを特徴とする、請求項７に記載の流体シミュレーションプログラム。
【請求項９】
前記時間刻みは、離散渦要素の移動速度と計算対象流れ場中の局所スケールを表す長さに基づいて自動的に最適化することを特徴とする、請求項７又は８に記載の流体シミュレーションプログラム。
【請求項１０】
請求項１乃至９に記載の流体シミュレーションプログラムを記録した記録媒体。
【発明の詳細な説明】
【０００１】
【発明の属する技術分野】
本発明は、数値流体力学　（　Ｃｏｍｐｕｔａｔｉｏｎａｌ　Ｆｌｕｉｄ　Ｄｙｎａｍｉｃｓ：　ＣＦＤ）　の手法を使用して、物体周りの流れを予測推定するための解析手段に関するものである。
【０００２】
【従来の技術】
一般に、コンピュータ技術の発達に伴って、あらゆる流体の挙動をその支配方程式から、計算によって予測しようとする試みが盛んに行われており、それらのアプローチを総称して数値流体力学（Ｃｏｍｐｕｔａｔｉｏｎａｌ　Ｆｌｕｉｄ　Ｄｙｎａｍｉｃｓ：ＣＦＤ）と呼んでいる。ＣＦＤの手法は、機械、航空、船舶ほか幅広い分野で起こる流体現象の予測推定に役立つものと期待されている。従来からＣＦＤ手法には、現象を表す偏微分方程式を数値的に計算するために、計算対象空間を微小領域に分割する数理的モデル化によって計算を行う、差分法や有限要素法が一般に用いられている。このように偏微分方程式を支配方程式とし、対象領域に計算格子を設けて解析する手法を総称して領域型解法という。この領域型解法は、解析する領域を有限個の離散化格子によって分割し、各格子点上で定義される物理量（速度、圧力、温度等）を、各格子点において支配方程式に基づいて数値計算で求めるものである。ここで、格子点とは、計算領域内で未知数が定義された点である。
【０００３】
上記従来の解析手法においては、対象空間を分割する計算格子を設ける必要があり、複雑な形状を伴う物体や複数の物体が移動するような対象を解析するためには、格子を形成するために多くの手間と時間を要するという問題がある。
【０００４】
このような計算格子を用いない解析方法として渦法を用いた解析方法が提案されており、物体壁面近傍に生成する渦度の挙動を流れの物理に基づいて直接的に取り扱うことにより、高精度かつ安定な流れ解析を実現する離散渦要素導入法も合わせて提案されている（例えば、非特許文献１、２参照。）。
【０００５】
【非特許文献１】
亀本喬司，　祝宝山，　小島成，　最新渦法の注目すべき特性とラグランジュＬＥＳ法としての課題，　第１４回数値流体力学シンポジウム，　Ｂ０６−４，　２０００．
【０００６】
【非特許文献２】
Ｏｊｉｍａ，　Ａ．　ａｎｄ　Ｋａｍｅｍｏｔｏ，　Ｋ．，　Ｎｕｍｅｒｉｃａｌ　ｓｉｍｕｌａｔｉｏｎ　ｏｆ　ｕｎｓｔｅａｄｙ　ｆｌｏｗｓ　ａｒｏｕｎｄ　ｔｈｒｅｅ　ｄｉｍｅｎｓｉｏｎａｌ　ｂｌｕｆｆ　ｂｏｄｉｅｓ　ｂｙ　ａｎ　ａｄｖａｎｃｅｄ　ｖｏｒｔｅｘ　ｍｅｔｈｏｄ，　ＪＳＭＥ　Ｉｎｔ．　Ｊｏｕｒｎａｌ，　Ｂ　４３−２　（２０００），　ｐｐ．１２７−１３５．
【０００７】
【発明が解決しようとする課題】
しかし、前記渦法を用いた解析方法では、渦層の近似高さと物体を形成するパネルの分割精度および時間精度との関係は明確にされておらず、最適な計算条件を設定することが困難であった。そのため、最適な流体のシミュレーションを得るには、渦層の近似高さ、物体を形成するパネルの分割、演算の時間間隔等のパラメータを変更しながら試行錯誤的に行わざるを得ないという問題がある。
【０００８】
そこで、本発明は、上記従来の問題点を解決し、対象区間を分割する計算格子の設定を不要とし、効率よく非定常な流体の挙動および物体に作用する力を見積もることのできる非定常流体特性のシミュレーションを行うことを目的とし、渦法を用いた流体解析において、渦層の近似高さおよびパネルサイズを自動的に最適化することを目的とする。
【０００９】
また、渦法を用いた流体シミュレーションにおいて、演算に要する計算時間は、殆どが離散渦要素の移動速度の算出時間が占めている。高精度の解析を行うためには多数の渦要素が必要となるが、渦要素の増加は計算時間を長期間化する要因となる。
【００１０】
そこで、本発明は、離散渦要素の移動速度の計算に要する演算回数を減らし、効率よく非定常な流体の挙動および物体に作用する力を見積もることができる非定常流体特性のシミュレーションを行うことを他の目的とする。
【００１１】
【課題を解決するための手段】
本発明は、前記目的を達成するために、空間中に配置される物体形状を境界要素法に基づいて考慮し、連続的な渦度の分布を多数の微小渦要素によって流れ場を離散的に表し、渦度輸送方程式を数値的に計算し、各離散渦要素の移動および代表する渦度の変化をラグランジュ的に追従することにより、流れをタイムスプリット的に非定常解析する演算を実行し、解析結果を表示装置へ表示することを特徴とするものであり、物体形状の表面データのみを入力することにより演算が可能であり、高速計算法、及び高精度解析法を適用することができる。また、解析結果を離散渦要素分布図、ベクトル図、あるいはカラーコンター図として出力することが好適である。ここで、離散渦要素分布図は瞬時の離散渦要素の位置もしくは離散渦要素が描く軌跡を表現したものであり、ベクトル図は空間の任意の位置における速度をベクトルにより表現したものである。また、カラーコンター図は空間の圧力分布を明るい色から暗い色への変化で表現したものである。
【００１２】
本発明の流体シミュレーションプログラムは、解析空間の境界を設定する座標位置データと流体特性を設定する特性データを入力データとし、物体表面の相対速度成分に基づく境界条件の下にビオ・サバールの式から、各境界における吹き出し量及び循環量に対する影響係数マトリクスを算出する機能と、影響係数マトリクスを含む連立方程式から前記吹き出し量及び循環量を算出する機能と、前記吹き出し量及び循環量から求まる速度場から、各境界近傍に渦層を近似して渦度を算出して離散渦要素を求める機能と、離散渦要素の移動速度をビオ・サバールの式から算出する機能と、離散渦要素を代表する渦度の時間変化を渦度輸送方程式に基づいて算出する機能とを備え、境界から発生した各離散渦要素の移動及び代表する渦度の時間変化をラグランジュ的に追従することにより流体現象をシミュレートする。さらに、前記流体シミュレーションから得られた速度場の情報に基づき圧力ポアソン方程式から任意の位置の圧力を算出する機能を備えることもできる。
【００１３】
また、座標位置データを変更する機能を有しており、この変更は物体に作用する力に基づいて行うことができる。座標位置データは時系列データとすることができ、これにより、移動・変形を伴う物体周りの流れをシミュレートすることができる。
【００１４】
本発明の流体シミュレーションプログラムは、離散渦要素の移動速度の算出において各離散渦要素間の相互作用を利用したり、離散渦要素の移動速度と空間スケールとの関係を考慮し演算頻度を自動的に最適化することにより、渦要素の移動速度の計算に要する演算回数を減らして演算速度を高めることができる。
【００１５】
また、本発明の流体シミュレーションプログラムは、渦度の生成モデルにおいて、渦層の近似高さ、パネル分割長さ、および時間刻みを自動的に最適化することができ、渦層の近似高さとパネル分割長さは、摩擦速度に基づいて自動的に最適化し、また、時間刻みは、離散渦要素の移動速度と計算対象流れ場中の局所スケールを表す長さに基づいて自動的に最適化する。
【００１６】
本発明によれば、単に物体形状を表す表面データのみを入力するだけで、物体による流体周りの挙動を解析することができる。解析結果は、離散渦要素分布図、ベクトル図またはカラーコンター図として表示することができ、視覚表示が容易である。
【００１７】
さらに、物体に作用する力の時系列データも出力し表示することができる。これらの簡易な入力、及び表示により、流体力学の専門知識のない技術者であっても流体解析シミュレーション結果を容易に理解することができる。
【００１８】
また、解析結果を検討した上で、記憶装置に保存されている形状データを用い、空間中に配置する物体の位置の修正や、パラメータの変更をデータ入力部より行うことにより効率良く流体計算をシミュレートする作用を営むことができる。
【００１９】
また、解析結果を、記憶装置に保存されている形状データにフィードバックし、空間中に配置する物体位置の修正や、パラメータの変更をデータ入力部より行うことにより効率よく流体現象をシミュレーションすることができる。
【００２０】
【発明の実施の形態】
以下、図を用いて、本発明に係る流体シミュレーション装置について詳細に説明する。
【００２１】
図１は、本発明が適用された流体シミュレーションプログラムの構成を示したブロック図である。図１に示すように、本プログラムは、データ入力部１、データ処理部２、および表示部３を備える。
【００２２】
データ処理部２は、データ入力部１で入力したデータを用いて流体の挙動を求め、表示部３に表示する。表示結果等により得られたデータを検討し、検討結果をデータ入力部１にフィードバックすることにより、最適な流体シミュレーションを行う。
【００２３】
データ入力部１は、流体解析に必要なデータを入力する。入力データは、流体解析を行う外形形状を特定する表面形状データである。表示形状データは、三次元データあるいは二次元データとすることができ、二次元データを入力データとする場合には、平面状の解析データが得られる。データ入力は、例えば、表示部３に表示された入力案内に基づきユーザーがキーボード等から行い、入力したデータをデータ処理部２に与える。
【００２４】
データ入力は、キーボード等の入力手段に限らず、記録媒体中に記録しておいたデータを読み込ませたり、あるいは、他の処理手段や制御手段から入力することもできる。この場合には、直接行った処理手段や制御手段から入力することも、ネットワークを介して行うこともできる。
【００２５】
データ処理部２は、形状設定部２ａ、境界条件設定部２ｂ、演算部２ｃ、シミュレーション部２ｄ等の各種処理部を備え、これら処理部によって流体解析を行う。その解析結果は、表示部３にグラフィック等により表示することができる。
【００２６】
図２は、データ処理部２の構成を説明するためのブロック図である。データ処理部２において、形状設定部２ａは、流体解析において空間に配置される物体の配置をデータ入力部１から与えられたデータに基づいて設定する。例えば、円柱周りの流れを解析するものとすると、その円柱の形状と円柱を配置する空間座標位置を設定する。
【００２７】
境界条件設定部２ｂは、形状設定部２ａによって設定された物体周りの流体の物性等をデータ入力部から与えられたデータに基づいて設定する。例えば、ユーザーは、物体の大きさ、流れの速度および流体の物性値に基づいて任意に設定することができ、例えば、計算対象流れ場の流入速度を設定する。
【００２８】
演算部２ｃは、流れの基礎方程式に基づき計算を行う。
【００２９】
以下、演算部２ｃで行う演算に用いる流れの基礎方程式について説明する。
【００３０】
本発明に係る流体解析方法の原理を説明する。流体解析は、図３に示すように、流れ場１０内に連続的に分布する渦度を微小な離散渦要素として取扱い、それら離散渦要素の挙動を追従することにより流れを表現する。
【００３１】
図３において、ある時刻ｔに物体１１の表面近傍で発生した渦度を離散的に取扱った離散渦要素２０は、時刻（ｔ＋１）において他の位置に移動し、時刻（ｔ＋２）においてさらに別の位置に移動する。図３では、時刻ｔにおいて渦度の生成に伴い流れ場に導入された離散渦要素の挙動を、それぞれ２０ａ〜２０ｅで表している。また、離散渦要素の移動と同時に、離散渦要素を代表する渦度も時間的に変化する。
【００３２】
これら各離散渦要素の位置および代表する渦度が時間的に変化することによって、流れ場の時間的な状態の変化を解析することができる。
【００３３】
ここで、渦度とは、局所的な流体の回転運動の強さを表す物理量であり、渦度をωとしたとき数（１）で表される。
ω＝ｒｏｔ　ｕ　　　　　　　　　　　　　　　　…（１）
なお、ｕは流速ベクトルを表している。
【００３４】
また、離散渦要素２０（２０ａ〜２０ｅ）は、図４に示すように、流れ場１０を分割した要素内の代表点２１およびその点を代表する渦度（図４では２２の矢印で模式的に表す）を表す要素である。この多数の離散渦要素２０の分布により流れ場１０は表される。
【００３５】
この流れ場１０において、渦度の分布は、零でない渦度ωを有する離散渦要素と渦度ωが零の離散渦要素とにより表すことができる。図４中の閉曲線２３は、流れ場１０において渦度ωを有する離散渦要素の範囲を示しており、流れ場の変化は、この離散渦要素の移動および代表する渦度の変化により表される。
【００３６】
以下、本発明に用いられる流体解析方法の基礎式について説明する。計算対象とする流体は、非圧縮粘性流体であり、以下の式（２）、および式（３）によって表される。
∂ｕ／∂ｔ＋（ｕ・ｇｒａｄ）ｕ＝−（ｇｒａｄ　ｐ）／ρ＋ν∇２ｕ　　　　…（２）
ｄｉｖ　ｕ＝０　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　…（３）
【００３７】
上記式において、ｐは圧力、ｔは時間、νは動粘性係数、ρは密度をそれぞれ表している。式（２）はナビエストークス方程式と呼ばれ、式（３）は連続の式と呼ばれる。
【００３８】
ここで用いる流体解析方法は、上述したように、渦度ωの時間変化を追従することにより流れを表現する。渦度ωの時間変化を追従することによる流体解析は、上記式（２）の回転をとり、式（３）を考慮することにより、以下の式（４）に示される渦度に関する偏微分方程式を導くことにより行うことができる。
∂ω／∂ｔ＋（ｕ・ｇｒａｄ）ω＝（ω・ｇｒａｄ）ｕ＋ｖ∇２ω　　　…（４）
【００３９】
式（４）は渦度輸送方程式と呼ばれ、渦度の時間的な変化を表している。式（１）および式（４）に注目すると、それぞれの式に圧力項は含まれておらず、渦度輸送方程式の未知数は、渦度ωを媒介変数とすれば流速ｕであり、流速ｕを媒介変数とすれば渦度ωとなる。つまり、式（１）および　　式（４）を与えられた境界条件のもとに解くことができれば、流速ｕあるいは渦度ωによって流れ場を規定することができる。
【００４０】
圧力に関しては、式（２）の発散をとることにより、式（５）に示す圧力に関する式が導かれる。
∇２ｐ＝−ρ　ｄｉｖ（ｕ・ｇｒａｄ　ｕ）　　　　　　　　　　　　　　　　　　　　　　　…（５）
【００４１】
式（５）は圧力ポアッソン方程式と呼ばれ、圧力は速度場より従属的に求めることが可能であることを意味する。
【００４２】
渦法の基礎式である式（１）は、流速ｕによって渦度ωを定義するものであるが、式（１）をベクトル的に積分することによりこの関係が逆転し、渦度ωによって流速ｕを定義する式（６）に示す積分方程式となる。
ｕ＝∫ｖ　ω×∇Ｇ　ｄｖ＋∫ｓ（（ｕ・ｎ）∇Ｇ−（ｕ×ｎ）×∇Ｇ）ｄｓ…（６）
【００４３】
なお、Ｇはラプラス方程式の基本解を表しており、次式により与えられる。
Ｇ＝−（１／２π）ｌｏｇ　Ｒ　　　　　　　　（２次元）　　　　　　　　　　　　…（７）
Ｇ＝１／（４πＲ）　　　　　　　　　　　　　　　　　（３次元）　　　　　　　　　　　　…（８）
【００４４】
式（６）は、一般化されたビオ・サバールの式と呼ばれる。式（６）の右辺第１項の空間積分は、流れ場に分布する渦度ωから誘起される流速を表しており、右辺第２項の境界積分は、境界上の流速条件によって決定される特異性をもつ関数から誘起される流速を表している。ここで、右辺第２項に示す分布関数は、境界上の法泉方向速度（ｕ・ｎ）および接線方向速度（ｕ×ｎ）にそれぞれに対応する吹出しσおよび循環γと等価である。したがって、流れ場の流速ｕは、図５に示すように空間に分布する渦度ωから誘起される流速ｕａと、境界に分布する吹出しσおよび循環γから誘起される流速ｕｂ，ｕｃの線形和からなっていることを意味する。
【００４５】
渦法における流れのモデル化は、式（６）の右辺第１項の渦度ωに関する空間積分を離散化して図４に示す多数の離散渦要素からの和として表し、第２項の境界積分を離散化して図６に示す吹出しσおよび循環γが分布する多数のパネル要素からの誘起速度として表現する。各パネル要素において、それぞれ吹出しσおよび循環γが設定される。
【００４６】
ここで、境界に分布する吹出しσおよび循環γの強さは、境界条件に基づいて決定される。このように境界に特異性をもつ関数を分布させ、その強さを境界条件によって決定する手法は、境界要素法の一種でパネル法と呼ばれる。
【００４７】
流れ場の時間変化は、図７に示すフローチャートの計算ステップを順次実行することにより表される。
【００４８】
ステップ１で、パネル法を用いて境界に分布する吹出しσあるいは循環γを算出する。ここで、図８に示すように、境界条件として境界上で流れの法線方向速度成分ｕ０・ｎを式（９）により規定し、これを境界条件とすることにより連立１次方程式を解き未知数σあるいはγを算出する。
ｕ０・ｎ＝ｖＢ・ｎ　　　　　　　　　　　　　　　　…（９）
【００４９】
ここで、ｕ０は境界上の速度、ｖＢは境界が移動する速度をそれぞれ表す。例えば、物体が静止している場合、境界は移動しないので右辺は零となる。境界上の点に誘起される速度は、式（６）より計算することができる。式（６）において、右辺第１項は渦度場からの影響であり既知である。
【００５０】
また、循環γについては、ケルビンの循環定理に従い、物体周りに分布する循環量の総和と流れ場に分布する総循環量が等しくなる条件を付加する必要がある。
【００５１】
物体周りに分布する吹出しσおよび循環γが規定されることにより、式（６）から速度場が既知となる。
【００５２】
続いて、ステップ２では、式（６）の右辺第１項に示す流れ場に存在する渦度ωの分布を表現するために、流れ場に離散渦要素を導入する。物体周りの流れは、物体から遠く離れた位置では主流速度で流れているのに対し、物体壁面上では速度が零となる。その間には、必ず速度の変化が存在する。言い換えれば、そこには渦度ωが存在するといえる。
【００５３】
図９は、二次元境界層の速度変化及び式（１）に示す渦度変化を説明するための図である。図９に示すように、物体壁面近傍に存在する流体粘性効果の大きな領域は境界層と呼ばれ、境界層内では急激な速度変化が生じているとともに（図９（ａ））、強い渦度が存在する（図９（ｂ））。
【００５４】
図１０は、離散渦要素の流れ場への導入を説明するための図である。本発明では、図１０に示すように、物体のごく近傍の速度分布を近似することにより、渦度の層（近似渦層）を仮定し、その渦層内の渦度の移流および拡散といった挙動を直接的に取扱い、離散渦要素を流れ場へ導入する手法を用いる。
【００５５】
図１０は、簡単のため最も基本的な渦層内を直線速度分布で近似したモデルを示している。図１０において、離散渦要素の流れ場への導入は、物体の表面から高さｈ（ｓ）における速度分布がｕ（ｓ，　ｎ）で近似される近似渦層をモデルとして想定する。近似渦層から流れ場への渦度の流出は、近似渦層内においてＡおよびＢに挟まれる微小空間に分布する流体塊の移動および粘性拡散に伴う拡がりにより表すことができる。
【００５６】
この手法を用いることにより、（１）壁面近傍の渦層の厚さｈ（ｓ）を調節することが可能となり、境界層の近似精度を自由に選択することができる。また、（２）壁面近傍の速度分布ｕ（ｓ，　ｎ）を近似することにより、壁面せん断応力を容易に見積もることができる。ここで用いる渦要素導入法は、上記の利点を備えていることから、計算機の容量、計算時間および計算精度に対する要求に柔軟に対応することが可能となる。
【００５７】
この離散渦要素の流れ場への導入において、高い計算精度、時間精度が得られる最適な計算条件とするために、渦層の近似高さおよびパネルサイズを最適化する。
【００５８】
はじめに、渦層の近似高さの最適化について説明する。
【００５９】
物体壁面のごく近傍では、層流・乱流に関わらず分子粘度が支配的な薄い粘性層が形成される。この薄い粘性層内では、式（１０）を満足するとみなすことができる。
τ＝μ・∂ｕ／∂ｎ＝ｃｏｎｓｔ．　　　　　　　　　　　　　　　　　　　　　　　　　…（１０）
【００６０】
ここで、τはせん断応力、μは粘性係数をそれぞれ表す。つまり、法線方向（ｎの方向）への速度勾配は一定であり、直線速度分布とみなすことができる。図１０中では、近似渦層において速度勾配を破線の矢印で示している。
【００６１】
ここで、粘性層の領域は、式（１１）に示すｙ＋により表すことができる。
ｙ＋＝ｖ＊ｙ／ν　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　…（１１）
【００６２】
ここで、ｖ＊は摩擦速度と呼ばれ、式（１２）によって定義される。
ｖ＊＝（τ０／ρ）１／２　　　　　　　　　　　　　　　　　　　　　　　　　　　　　…（１２）
ここで、τ０は、壁面上のせん断応力である。
【００６３】
図１０において、渦層内の速度分布を近似した際に算出された速度勾配を用いることにより、式（１３）に基づいて渦層高さｈを決定する。
ｈ＝ｙ＋ν／ｖ＊　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　…（１３）
【００６４】
特に流れが乱流のとき、前記壁面近傍の薄い粘性層は粘性底層と呼ばれ、その領域はｙ＋＜５とされる。高精度の解を求める場合、近似する渦層高さは粘性底層ほどの高さに設定することが望ましいが、求める解析精度に応じて任意に設定することができる。例えば、高い解析精度を求める場合には小さなｙ＋を設定し、高い解析精度が不要である場合には大きなｙ＋を設定する。
【００６５】
したがって、ｙ＋を解析精度の指標として設定することにより、式（１３）から渦層高さｈを自動的に決定することができる。
【００６６】
次に、パネルサイズの最適化について説明する。パネルサイズは、図１０において、局所的なパネル長さΔｓの長短により定めることができ、パネル長さΔｓを短く設定することにより解析精度は高めることができる。
【００６７】
渦層高さｈが設定された状態において、解析精度を保つためには、近似した渦層高さｈと同じ程度のスケールでパネル分割を行う必要があり、局所的なパネル長さΔｓを、Δｓ＝ｈとなるように設定する。
【００６８】
したがって、ｙ＋を解析精度の指標として設定することにより、Δｓ＝ｈとしてパネル長さΔｓを自動的に決定することができる。ここでは、Δｓ＝ｈに設定しているが、Δｓ＝ｋｈとして任意の係数ｋにより調整するようにしてもよい。
【００６９】
なお、渦層高さｈ及びパネル長さΔｓを定める物体壁面のごく近傍の位置は、入力した物体の形状データに基づいて定めることができる。また、物体壁面のごく近傍の各位置に応じてｙ＋を設定し、渦層高さｈ及びパネル長さΔｓを定めるようにしてもよい。
【００７０】
このように、本発明では、ｙ＋を解析精度の指標として設定することにより、渦層高さｈ、及びパネル長さΔｓを自動的に決定することができ、流れの状況にあった渦層近似を行うと同時に、自動的に最適なパネル分割が可能となる。
【００７１】
次に、時間刻みの設定について説明する。
離散渦要素による流れ解析において、流れ場の時間変化を求めるには最適な時間精度で離散渦要素の移動を計算する必要がある。
【００７２】
従来、各計算ステップにおいて全離散渦要素の移動速度を計算し、同一の時間刻みで移動させていた。流れ場は、局所的に様々な速度分布となるため、同一の時間刻みで離散渦要素を移動すると、局所的な流れの状況により解析精度へ影響を与える可能性が考えられる。つまり、局所的に速度変化の激しい領域では、時間刻みが大き過ぎるために離散渦要素の正確な挙動が表現できず、解析精度を著しく損なう恐れがある。また、ある離散渦要素が、毎ステップ同じ速度で移動していたとしても、離散渦要素の移動速度を算出する演算は、毎ステップで行われることになり、無駄に演算を重ねることにもなる。
【００７３】
そこで、局所的な流れの状況に応じて時間刻みを調節し、最適な時間刻みで離散渦要素の挙動を追跡することより、局所的に流れの変化の激しい領域での解析精度を保持するとともに、離散渦要素の移動速度を計算する演算回数を調節することを目的とする。
【００７４】
この局所的な流れの様相に合わせて時間刻みを決定するために、流れ場の任意の位置に置かれた離散渦要素の移動速度ｕが、式（１４）の関係を満足するように時間刻みΔｔを決定する。
ｕΔｔ≦Ｌ　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　…（１４）
【００７５】
なお、Ｌは計算対象とする流れ場中の局所スケールを表す長さとするが、求める解析精度に応じて任意に設定することができる。例えば、Ｌはパネル長さと設定することができる。
【００７６】
次に、ステップ３において、離散渦要素の位置および代表する渦度を時間変化させる。
【００７７】
離散渦要素の移動する速度は、式（６）に示すビオ・サバールの式より求めることができる。また、離散渦要素の代表する渦度の変化は、式（４）に示す渦度輸送方程式に基づいた数理的モデルを用いて計算する。例えば、前記数理モデルのうち粘性拡散に対してはコアスプレディング法と呼ばれる手法があり、これは時間の経過とともに離散渦要素の粘性核半径を拡大することによって離散渦要素が代表する渦度の減衰を表現するものである。なお、コアスプレディング法は、Ｌｅｏｎａｒｄ，　Ａ．，　Ｖｏｒｔｅｘ　ｍｅｔｈｏｄｓ　ｆｏｒ　ｆｌｏｗ　ｓｉｍｕｌａｔｉｏｎｓ，　Ｊ．　Ｃｏｍｐ．　Ｐｈｙｓ．　３７，　（１９８０），　ｐｐ．２８９−３３５．に知られている。
【００７８】
この離散渦要素の移動速度を計算するための演算回数が渦法の計算において多くの時間を要することとなる。
【００７９】
そこで、本発明は離散渦要素の相互作用を利用することにより、離散渦要素の移動速度の算出に関する演算回数を削減する。以下、この演算回数を削減するための方法を、図１１の離散渦要素に誘起される速度を示した図を用いて説明する。
【００８０】
図１１（ａ）に示すように、ある離散渦要素に誘起される速度は、その他の離散渦要素から誘起される速度の積分値として算出される。この関係は以下の式（１５）で表される。
ｕ＝Σ∫ｖω×∇Ｇｄｖ　　　　　　　　　　　　　　　…（１５）
【００８１】
したがって、ある１つの離散渦要素に誘起される速度を求めるためには、離散渦要素がＮ個存在するならばＮ回の演算を要することになる。Ｎ個の全離散渦要素の速度を求めるためにはＮ２回の演算を要することになる。
【００８２】
しかし、式（１５）の∇Ｇは、式（７）および式（８）に示すように位置の関数であることより、各離散渦要素の相互作用は対称となる（図１１（ｂ））。したがって、ある離散渦要素から誘起される速度を計算するのと同時に相手の離散渦要素へ誘起する速度も計算することが可能である。例えば、２つの離散渦要素として離散渦要素Ｐ，Ｑがある場合、離散渦要素Ｑから離散渦要素Ｐに誘起される速度を求める際に算出した∇Ｇと、離散渦要素Ｐから離散渦要素Ｑに誘起される速度を求める際に算出される∇Ｇは、方向が逆で同じ大きさとなる。そこで、例えば、離散渦要素Ｑから離散渦要素Ｐへ向かうベクトルを規定し∇Ｇを算出すれば、この∇Ｇの値は、離散渦要素Ｑから離散渦要素Ｐに誘起される速度を算出するために利用するだけではなく、離散渦要素Ｐから離散渦要素Ｑに誘起される速度の計算に利用することができる。
【００８３】
本発明は、この計算方法を用いることにより、離散渦要素数Ｎ個について要する演算回数はＮ２／２回に低減することができる。
【００８４】
次に、ステップ４において、必要に応じて式（４）から任意の位置における圧力を計算する。一般的には、式（５）を差分的に解く方法が用いられているが、その差分法による場合には、流れ場に微係数を定義するための計算格子を設ける必要があり、本発明の計算格子が不要であるという利点が損なわれてしまう。
【００８５】
そこで、本発明は、式（５）を積分方程式に変換することにより、流れ場に計算格子を設けることなく圧力を計算する手法を用いる。ベルヌーイ関数を式（１６）のように定義し、式（５）を流れ場で積分すると式（１７）が導かれる。
Ｈ＝ｐ／ρ＋ｕ２／２　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　…（１６）
βＨ＋∫ｓＨ（∂Ｇ／∂ｎ）
＝∫ｖ　∇Ｇ（ｕ×ω）ｄｖ
−∫ｓ（Ｇｎ・∂ｕ／∂ｎ＋νｎ・（　∇Ｇ×ω））ｄｓ　　…（１７）
【００８６】
なお、βは、領域内で１であり、境界面上で１／２である。また、Ｇ　は式（７）および式（８）に示すラプラス方程式の基本解である。
【００８７】
式（１７）の右辺の各項は、流れ場に分布する渦度および速度が既知であれば計算することが可能である。したがって式（１７）について、ベルヌーイ関数Ｈを未知数とした連立１次方程式として解くことが可能である。ベルヌーイ関数Ｈが求められたならば、式（１６）よりただちに圧力を算出することができる。
【００８８】
なお、上記式の導出に関しては、　Ｕｈｌｍａｎ，　Ｊ．Ｓ．，　Ａｎ　ｉｎｔｅｇｒａｌ　ｅｑｕａｔｉｏｎ　ｆｏｒｍｕｌａｔｉｏｎ　ｏｆ　ｔｈｅ　ｅｑｕａｔｉｏｎ　ｏｆ　ｍｏｔｉｏｎ　ｏｆ　ａｎ　ｉｎｃｏｍｐｒｅｓｓｉｂｌｅ　ｆｌｕｉｄ，　Ｎａｖａｌ　Ｕｎｄｅｒｓｅａ　Ｗａｒｆａｒｅ　Ｃｅｎｔｅｒ　Ｔ．Ｒ．，　１０−０８６，　（１９９２）．に示されている。
【００８９】
以上の第１から第４の各操作ステップを繰り返し実行することにより、計算格子を生成することなく流体の挙動が計算できる。
【００９０】
次に、移動・変形を伴う物体周りの流れについて説明する。
空気などの流体が流れている中にある物体は、流体から力を受ける。このとき、物体の構造特性と流れ場の状態の相互作用によって振動等が発生することがある。機械を設計する上で振動の発生を事前に予測することは非常に重要である。本発明は、流体解析によって物体に作用する力を算出し、それを用いて構造解析を行うことにより流体現象に起因する振動現象の予測を行うことができる。
【００９１】
前記した方法により物体表面に作用する圧力ｐおよび式（１０）に基づきせん断力τを求め、物体表面に作用する圧力ｐおよびせん断力τを表面積分することにより、各方向に作用する力を算出することができる。
【００９２】
物体の移動および変形は、前記計算により求められた物体に作用する力を用い、運動方程式を解くことにより計算することができる。本発明では、運動方程式より得られた解に基づいて物体の位置もしくは形状を変化させ、その時間変化率を境界条件として流れを解析する。
【００９３】
具体的には、物体周りに選択された各選点の移動速度ｖＢの法線方向成分が、式（９）の右辺に示す項となる。
【００９４】
物体の移動・変形は、図１２のフローチャートに示す計算ステップを順次実行することにより行われる。
【００９５】
ステップ１１で、前記の流れ解析を実施して得られた壁面での圧力及びせん断力を物体周りに積分することにより、物体に作用する力を算出する。図１３は、物体に作用する力を説明するための図である。図１３に示すように、圧力ｐは壁面に垂直に作用し、せん断力τは壁面に平行に作用することより、物体周りで選んだ各選点に作用する単位面積あたりの力ｆは、式（１８）のようになる。
ｆ＝ｐｎ　＋τｓ　　　　　　　　　　　　　　　　　　　…（１８）
【００９６】
したがって、物体に作用する力Ｆは、ｆを物体周りに積分することにより算出することができる。
【００９７】
ステップ１２では、ステップ１１で求められた力が物体に作用することにより、物体が移動もしくは変形する変化量を式（１９）で算出する。
Ｍ（ｄ２Ｘ／ｄｔ２）　＋　Ｃ（ｄＸ／ｄｔ）＋　ＫＸ　＝Ｆ　…（１９）
【００９８】
式（１９）は、運動方程式と呼ばれ、図１４に示すばねとダンパーを含む基本的な運動モデルで表される。ここで、Ｍ，Ｃ，およびＫは、質量、減衰係数およびばね係数をそれぞれ表し、Ｆは物体に作用する外力を表す。Ｍ、Ｃ、およびＫは物体固有の量であるので、外力Ｆが流体解析より求められた場合には、式（１９）は数値的に解くことができ、変位Ｘを算出することができる。
【００９９】
ステップ１３では、ステップ１２で求められた変位に従い、図１５に示すように物体を表現する点列データの各位置座標を移動する。点列データの各位置座標の移動速度が、次ステップで流れ解析を行う際の境界条件となる。
【０１００】
以上の各ステップを繰り返し実行することにより、移動・変形を伴う物体周りの流れをシミュレートすることができる。
【０１０１】
図２において、以上の演算処理が終了すると、シミュレーション処理部２ｄは、記憶装置の出力データ領域のデータを参照しながら表示部３に表示する画像データを作成する。
【０１０２】
図１６はシミュレーション結果の一例であり、表示部３にシミュレーション結果を表示することにより確認することができる。これらのシミュレーション結果を時間間隔毎に、カラープリンターなどに印字して設計資料として保存することもできる。
【０１０３】
以上のように、本実施例の流体特性シミュレーションによれば、被計算対象物の干渉チェックが可能なため、正確な流体計算が行えると共に、誤って形状作成を行うことを抑制することができ、オペレータは形状チェックなどを行う必要がなく負担軽減が図れる。なお、本実施例においては渦要素分布図で演算結果を画面表示したがベクトル図やカラーコンター図あるいは他の表示であってもかまわない。また、形状データの入力は点列データ入力方式　としたがＣＡＤなどの作図データを直接入力するものであってもよい。すなわち本発明は、上記実施例に限定されるものではなく、本発明の趣旨に基づき種々の変形が可能であり、これらを本発明の範囲から排除するものではない。
【０１０４】
【発明の効果】
本発明によれば、物体形状を表面データの入力のみで演算を行うと共に、高速解法および高精度解法を利用し、効率よく非定常な流体の挙動および物体に作用する力を見積もることのできる非定常流体特性のシミュレーション装置を提供できる。
【図面の簡単な説明】
【図１】本発明が適用された流体シミュレーションプログラムの構成を示したブロック図である。
【図２】本発明のデータ処理部の構成を説明するためのブロック図である。
【図３】本発明に係る流体解析方法の原理を説明するための図である。
【図４】本発明の離散渦要素を説明するための図である。
【図５】流れ場の流速を説明するための図である。
【図６】パネル要素における吹出しσおよび循環γの分布を説明するための図である。
【図７】本発明の流れ場の時間変化の計算手順を説明するためのフローチャートである。
【図８】本発明の流れ場の時間変化の計算に用いる境界条件を説明するための図である。
【図９】本発明の境界層の速度変化及び渦度変化を説明するための図である。
【図１０】本発明の離散渦要素の流れ場への導入を説明するための図である。
【図１１】本発明の離散渦要素に誘起される速度を説明ための図である。
【図１２】本発明の物体の移動・変形の計算手順を説明するためのフローチャートである。
【図１３】物体に作用する力を説明するための図である。
【図１４】ばねとダンパーを含む基本的な運動モデルを示す図である。
【図１５】本発明の物体を表現する点列データの各位置座標の移動を示す図である。
【図１６】本発明のシミュレーション結果の一例を示す図である。
【符号の説明】
１　　データ入力部
２　　データ処理部
３　　表示部
２ａ　　形状設定部
２ｂ　　境界条件設定部
２ｃ　　演算部
２ｄ　　シミュレーション部
ｕ　　速度
γ　　循環
σ　　吹出し
ω　　渦度
ν　　動粘性係数
μ　　粘性係数
ρ　　密度
ｐ　　圧力
【図１】
2004126925.tif 000003
【図２】
2004126925.tif 000004
【図３】
2004126925.tif 000005
【図４】
2004126925.tif 000006
【図５】
2004126925.tif 000007
【図６】
2004126925.tif 000008
【図７】
2004126925.tif 000009
【図８】
2004126925.tif 000010
【図９】
2004126925.tif 000011
【図１０】
2004126925.tif 000012
【図１１】
2004126925.tif 000013
【図１２】
2004126925.tif 000014
【図１３】
2004126925.tif 000015
【図１４】
2004126925.tif 000016
【図１５】
2004126925.tif 000017
【図１６】
2004126925.tif 000018