JPWO2004061723A1 - Ｖ−ｃａｄデータを直接用いた非圧縮性粘性流体の流れ場の数値解析方法と装置 - Google Patents
Ｖ−ｃａｄデータを直接用いた非圧縮性粘性流体の流れ場の数値解析方法と装置 Download PDFInfo
- Publication number
- JPWO2004061723A1 JPWO2004061723A1 JP2004564488A JP2004564488A JPWO2004061723A1 JP WO2004061723 A1 JPWO2004061723 A1 JP WO2004061723A1 JP 2004564488 A JP2004564488 A JP 2004564488A JP 2004564488 A JP2004564488 A JP 2004564488A JP WO2004061723 A1 JPWO2004061723 A1 JP WO2004061723A1
- Authority
- JP
- Japan
- Prior art keywords
- boundary
- cell
- data
- flow field
- equation
- Prior art date
- Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
- Granted
Links
Images
Classifications
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06F—ELECTRIC DIGITAL DATA PROCESSING
- G06F30/00—Computer-aided design [CAD]
- G06F30/20—Design optimisation, verification or simulation
- G06F30/23—Design optimisation, verification or simulation using finite element methods [FEM] or finite difference methods [FDM]
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06F—ELECTRIC DIGITAL DATA PROCESSING
- G06F2111/00—Details relating to CAD techniques
- G06F2111/10—Numerical modelling
Landscapes
- Engineering & Computer Science (AREA)
- Physics & Mathematics (AREA)
- Theoretical Computer Science (AREA)
- Computer Hardware Design (AREA)
- Evolutionary Computation (AREA)
- Geometry (AREA)
- General Engineering & Computer Science (AREA)
- General Physics & Mathematics (AREA)
- Management, Administration, Business Operations System, And Electronic Commerce (AREA)
- Complex Calculations (AREA)
Abstract
非圧縮性粘性流体と接する対象物の境界データからなる外部データ１２を境界が直交する複数のセル１３に分割する分割ステップ（Ａ）と、分割された各セルを対象物の内側に位置する内部セル１３ａと境界データを含む境界セル１３ｂとに区分するセル区分ステップ（Ｂ）と、境界データによる境界セル１３ｂの稜線の切断点を求める切断点決定ステップ（Ｃ）と、求めた切断点を結ぶ多角形を境界面のセル内部データとする境界面決定ステップＤと、流れ場の境界に、ＶＯＦ法を併用したカットセル有限体積法を適用して解析する解析ステップ（Ｅ）とを備える。
Description
発明の背景
発明の技術分野
本発明は、形状と物理量を統合した実体データを記憶するＶ−ＣＡＤデータを直接用いた非圧縮性粘性流体の流れ場の数値解析方法と装置に関する。
関連技術の説明
形状と物性を統合した実体データを小さい記憶容量で記憶することができ、これにより、物体の形状・構造・物性情報・履歴を一元的に管理し、設計から加工、組立、試験、評価など一連の工程に関わるデータを同じデータで管理することができ、ＣＡＤとシミュレーションを一元化することできる実体データの記憶方法として、［特許文献１］が開示されている。
特開２００２−２３００５４号公報
発明の技術分野
本発明は、形状と物理量を統合した実体データを記憶するＶ−ＣＡＤデータを直接用いた非圧縮性粘性流体の流れ場の数値解析方法と装置に関する。
関連技術の説明
形状と物性を統合した実体データを小さい記憶容量で記憶することができ、これにより、物体の形状・構造・物性情報・履歴を一元的に管理し、設計から加工、組立、試験、評価など一連の工程に関わるデータを同じデータで管理することができ、ＣＡＤとシミュレーションを一元化することできる実体データの記憶方法として、［特許文献１］が開示されている。
［特許文献１］の「形状と物性を統合した実体データの記憶方法」は、図１に示すように、外部データ入力ステップ（Ａ）、八分木分割ステップ（Ｂ）、及びセルデータ記憶ステップ（Ｃ）からなり、外部データ入力ステップ（Ａ）では、外部データ取得ステップＳ１で取得した対象物の境界データからなる外部データ１２を本発明の方法を記憶したコンピュータ等に入力し、八分木分割ステップ（Ｂ）では、外部データ１２を八分木分割により境界平面が直交する直方体のセル１３に分割し、セルデータ記憶ステップ（Ｃ）では、各セル毎に種々の物性値を記憶するものである。
上述した［特許文献１］の発明は、対象物の形状データからなる外部データを、八分木分割により境界平面が直交する直方体のセルに分割し、各セル毎に種々の物性値を記憶するものである。分割された各セルは対象物の内側に位置する内部セルと、境界面を含む境界セルとからなる。また、内部セルは、属性として１種の物性値を持ち、境界セルは、対象物の内側と外側の２種の物性値をもつものである。
以下、この方法によるデータを「Ｖ−ＣＡＤデータ」と呼び、このデータを用いた設計やシミュレーションを「ボリュームＣＡＤ」又は「Ｖ−ＣＡＤ」と呼ぶ。図１において１４がＶ−ＣＡＤデータである。
ＣＦＤ（Ｃｏｍｐｕｔａｔｉｏｎａｌ Ｆｌｕｉｄ Ｄｙｎａｍｉｃｓ）が実用化するに従って格子生成に手間や時間がかかり、複雑形状では計算時間よりも格子生成の時間の方が長くなってきている。このため、近年、直交格子による流体解析が再び話題となっている。直交格子による流体解析に関しては、［非特許文献１］〜［非特許文献１７］が知られている。
また、「強制振動円柱まわりの流れ」に関しては、［非特許文献１８］に実験結果が開示され、「円柱の渦放出による自励振動」に関しては、［非特許文献１９］にＡＬＥ有限要素法による計算結果が開示されている。
Ｓａｉｋｉ，Ｅ．Ｍ．，Ｂｉｒｉｎｇｅｎ，Ｓ．，１９９６，Ｎｕｍｅｒｉｃａｌ Ｓｉｍｕｌａｔｉｏｎ ｏｆ ａ Ｃｙｌｉｎｄｅｒ ｉｎ Ｕｎｉｆｏｒｍ ｆｌｏｗ：Ａｐｐｌｉｃａｔｉｏｎ ｏｆ ａ Ｖｉｒｔｕａｌ Ｂｏｕｎｄａｒｙ Ｍｅｔｈｏｄ，Ｊ．Ｃｏｍｐｕｔ．Ｐｈｙｓ．１２３，４５０−４６５． 矢部孝・肖鋒，１９９７，固体・液体・気体の統一解法とＣＩＰ法（２），数値流体力学，７，１０３−１１４． Ｙｅ，Ｔ．，Ｍｉｔｔａｌ，Ｒ．，Ｕｄａｙｋｕｍａｒ，Ｈ．Ｓ．，＆ Ｓｈｖｙ，Ｗ．，１９９９，Ａ Ｃａｒｔｅｓｉａｎ Ｇｒｉｄ Ｍｅｔｈｏｄ ｆｏｒ Ｖｉｓｃｏｕｓ Ｉｎｃｏｍｐｒｅｓｓｉｂｌｅ Ｆｌｏｗｓ ｗｉｔｈ Ｃｏｍｐｌｅｘ Ｉｍｍｅｒｓｅｄ Ｂｏｕｎｄａｒｉｅｓ，ＡＩＡＡ−９９−３３１２，５４５−５５７． 中村明・下村信雄・里深信行，１９９５，デカルト格子系による任意形状物体周りの圧縮性粘性流計算，日本機械学会論文集，６１Ｂ−５９２，４３１９−４３２６． 市川治・藤井孝蔵，２００２，直交格子を使用した三次元の任意形状物体まわりの流体シミュレーション，日本機械学会論文集，６８Ｂ−６６９，１３２９−１３３６． 朴炳湖・黒田成昭，２０００，非圧縮性粘性流れの直交格子解法，ながれ，１９，３７−４６． Ｏｎｏ，Ｋ．，Ｔｏｍｉｔａ，Ｎ．，Ｆｕｊｉｔａｎｉ，Ｋ．，＆ Ｈｉｍｅｎｏ，Ｒ．，１９９８，Ａｎ Ａｐｐｌｉｃａｔｉｏｎ ｏｆ Ｖｏｘｅｌ Ｍｏｄｅｌｉｎｇ Ａｐｐｒｏａｃｈ ｔｏ Ｐｒｅｄｉｃｔｉｏｎ ｏｆ Ｅｎｇｉｎｅ Ｃｏｏｌｉｎｇ Ｆｌｏｗ，Ｓｏｃｉｅｔｙ ｏｆ Ａｕｔｏｍｏｔｉｖｅ Ｅｎｇｉｎｅｅｒｓ ｏｆ Ｊａｐａｎ，Ｓｐｒｉｎｇ Ｃｏｎｖｅｎｔｉｏｎ，Ｎｏ．９８４，１６５−１６８． ｈｔｔｐ：／／ｋｕｗａｈａｒａ．ｉｓａｓ．ａｃ．ｊｐ／ｉｎｄｅｘ．ｈｔｍｌ 寺本進・藤井孝蔵，１９９８，直交格子法による三次元物体周りの流れ解析，第１２回数値流体力学シンポジウム講演論文集，２９９−３００． Ｑｕｉｒｋ，Ｊ．Ｊ．，１９９４，Ａｎ Ａｌｔｅｒｎａｔｉｖｅ ｔｏ Ｕｎｓｔｒｕｃｔｕｒｅｄ Ｇｒｉｄｓ ｆｏｒ Ｃｏｍｐｕｔｉｎｇ Ｇａｓ Ｄｙｎａｍｉｃ Ｆｌｏｗｓ Ａｒｏｕｎｄ Ａｒｂｉｔｒａｒｉｌｙ Ｃｏｍｐｌｅｘ Ｔｗｏ−Ｄｉｍｅｎｓｉｏｎａｌ Ｂｏｄｉｅｓ，Ｃｏｍｐｕｔｅｒｓ Ｆｌｕｉｄｓ，２３，１２５−１４２． Ｋａｒｍａｎ，Ｓ．Ｌ．Ｊｒ．，１９９５，ＳＰＬＩＴＦＬＯＷ：Ａ ３Ｄ Ｕｎｓｔｒｕｃｔｕｒｅｄ Ｃａｒｔｅｓｉａｎ／Ｐｒｉｓｍａｔｉｃ Ｇｒｉｄ （１２） ｙｎａｍｉｃｓ ｏｆ ＣＦＤ Ｃｏｄｅ ｆｏｒ Ｃｏｍｐｌｅｘ Ｇｅｏｍｅ−ｔｒｉｅｓ，ＡＩＡＡ ９５−０３４３． Ｈｉｒｔ，Ｃ．Ｗ．，＆ Ｎｉｃｈｏｌｓ，Ｂ．Ｄ．，１９８１，Ｖｏｌｕｍｅ ｏｆ Ｆｌｕｉｄ（ＶＯＦ） Ｍｅｔｈｏｄ ｆｏｒ ｔｈｅ Ｄ Ｆｒｅｅ Ｂｏｕｎｄａｒｉｅｓ，Ｊ．Ｃｏｍｐｕｔ．Ｐｈｙｓ．３９，２０１−２２５． Ｈｉｒｔ，Ｃ．Ｗ．，＆ Ｃｏｏｋ，Ｊ．Ｌ．，１９７２，Ｃａｌｃｕｌａｔｉｎｇ Ｔｈｒｅｅ−ｄｉｍｅｎｓｉｏｎａｌ Ｆｌｏｗｓ Ａｒｏｕｎｄ Ｓｔｒｕｃｔｕｒｅｓ ａｎｄ Ｏｖｅｒ Ｒｏｕｇｈ Ｔｅｒｒａｉｎ，Ｊ．Ｃｏｍｐｕｔ．ｐｈｙｓ．１０，３２４−３４０． 加瀬究・手嶋吉法，２００１，ボリュームＣＡＤの開発，理研シンポジウム・ものつくり情報技術統合化研究，第１回，６−１１． 豊田郁夫・荒川忠一，１９９９，直交適合格子法による円柱周り流れの解析，第１３回数値流体力学シンポジウム，Ｆ０３−１，ＣＤ−ＲＯＭ． 松宮輝・木枝香織・谷口伸行・小林敏雄，１９９３，三次精度風上差分法による二次元円柱周り流れの数値シミュレーション，日本機械学会論文集，５９Ｂ−５６６，２９３７−２９４３． Ｂｏｕａｒｄ，Ｒ．，＆ Ｃｏｕｔａｎｃｅａｕ，Ｍ．，１９８０，Ｔｈｅ ｅａｒｌｙ ｓｔａｇｅ ｏｆ ｄｅｖｅｌｏｐｍｅｎｔ ｏｆ ｔｈｅ ｗａｋｅ ｂｅｈｉｎｄ ａｎ ｉｍｐｕｌｓｉｖｅｌｙ ｓｔａｒｔｅｄ ｃｙｌｉｎｄｅｒ ｆｏｒ ４０＜Ｒｅ＜１０＾４，Ｊ．Ｆｌｕｉｄ Ｍｅｃｈ．，１０１−３，５８３−６０７． Ｏｋａｍｏｔｏ，Ｓ．，Ｕｅｍａｔｓｕ，Ｒ．，ａｎｄ Ｔａｇｕｗａ，Ｙ．，Ｆｌｕｉｄ ｆｏｒｃｅ ａｃｔｉｎｇ ｏｎ ｔｗｏ−ｄｉｍｅｎｓｉｏｎａｌ ｃｉｒｃｕｌａｒ ｃｙｌｉｎｄｅｒ ｉｎ Ｌｏｋ−ｉｎ ｐｈｅｎｏｍｅｎｏｎ，ＪＳＭＥ Ｉｎｔｅｒｎａｔｉｏｎａｌ Ｊｏｕｒｎａｌ，Ｂ４５，Ｎｏ．４，（２００２），８５０−８５６． 近藤典夫，円柱の空力挙動に関する数値シミュレーション，第１５回数値流体力学シンポジウム，Ｅ０９−２，（２００１），ＣＤ−ＲＯＭ．
上述した［特許文献１］の発明は、対象物の形状データからなる外部データを、八分木分割により境界平面が直交する直方体のセルに分割し、各セル毎に種々の物性値を記憶するものである。分割された各セルは対象物の内側に位置する内部セルと、境界面を含む境界セルとからなる。また、内部セルは、属性として１種の物性値を持ち、境界セルは、対象物の内側と外側の２種の物性値をもつものである。
以下、この方法によるデータを「Ｖ−ＣＡＤデータ」と呼び、このデータを用いた設計やシミュレーションを「ボリュームＣＡＤ」又は「Ｖ−ＣＡＤ」と呼ぶ。図１において１４がＶ−ＣＡＤデータである。
ＣＦＤ（Ｃｏｍｐｕｔａｔｉｏｎａｌ Ｆｌｕｉｄ Ｄｙｎａｍｉｃｓ）が実用化するに従って格子生成に手間や時間がかかり、複雑形状では計算時間よりも格子生成の時間の方が長くなってきている。このため、近年、直交格子による流体解析が再び話題となっている。直交格子による流体解析に関しては、［非特許文献１］〜［非特許文献１７］が知られている。
また、「強制振動円柱まわりの流れ」に関しては、［非特許文献１８］に実験結果が開示され、「円柱の渦放出による自励振動」に関しては、［非特許文献１９］にＡＬＥ有限要素法による計算結果が開示されている。
Ｓａｉｋｉ，Ｅ．Ｍ．，Ｂｉｒｉｎｇｅｎ，Ｓ．，１９９６，Ｎｕｍｅｒｉｃａｌ Ｓｉｍｕｌａｔｉｏｎ ｏｆ ａ Ｃｙｌｉｎｄｅｒ ｉｎ Ｕｎｉｆｏｒｍ ｆｌｏｗ：Ａｐｐｌｉｃａｔｉｏｎ ｏｆ ａ Ｖｉｒｔｕａｌ Ｂｏｕｎｄａｒｙ Ｍｅｔｈｏｄ，Ｊ．Ｃｏｍｐｕｔ．Ｐｈｙｓ．１２３，４５０−４６５． 矢部孝・肖鋒，１９９７，固体・液体・気体の統一解法とＣＩＰ法（２），数値流体力学，７，１０３−１１４． Ｙｅ，Ｔ．，Ｍｉｔｔａｌ，Ｒ．，Ｕｄａｙｋｕｍａｒ，Ｈ．Ｓ．，＆ Ｓｈｖｙ，Ｗ．，１９９９，Ａ Ｃａｒｔｅｓｉａｎ Ｇｒｉｄ Ｍｅｔｈｏｄ ｆｏｒ Ｖｉｓｃｏｕｓ Ｉｎｃｏｍｐｒｅｓｓｉｂｌｅ Ｆｌｏｗｓ ｗｉｔｈ Ｃｏｍｐｌｅｘ Ｉｍｍｅｒｓｅｄ Ｂｏｕｎｄａｒｉｅｓ，ＡＩＡＡ−９９−３３１２，５４５−５５７． 中村明・下村信雄・里深信行，１９９５，デカルト格子系による任意形状物体周りの圧縮性粘性流計算，日本機械学会論文集，６１Ｂ−５９２，４３１９−４３２６． 市川治・藤井孝蔵，２００２，直交格子を使用した三次元の任意形状物体まわりの流体シミュレーション，日本機械学会論文集，６８Ｂ−６６９，１３２９−１３３６． 朴炳湖・黒田成昭，２０００，非圧縮性粘性流れの直交格子解法，ながれ，１９，３７−４６． Ｏｎｏ，Ｋ．，Ｔｏｍｉｔａ，Ｎ．，Ｆｕｊｉｔａｎｉ，Ｋ．，＆ Ｈｉｍｅｎｏ，Ｒ．，１９９８，Ａｎ Ａｐｐｌｉｃａｔｉｏｎ ｏｆ Ｖｏｘｅｌ Ｍｏｄｅｌｉｎｇ Ａｐｐｒｏａｃｈ ｔｏ Ｐｒｅｄｉｃｔｉｏｎ ｏｆ Ｅｎｇｉｎｅ Ｃｏｏｌｉｎｇ Ｆｌｏｗ，Ｓｏｃｉｅｔｙ ｏｆ Ａｕｔｏｍｏｔｉｖｅ Ｅｎｇｉｎｅｅｒｓ ｏｆ Ｊａｐａｎ，Ｓｐｒｉｎｇ Ｃｏｎｖｅｎｔｉｏｎ，Ｎｏ．９８４，１６５−１６８． ｈｔｔｐ：／／ｋｕｗａｈａｒａ．ｉｓａｓ．ａｃ．ｊｐ／ｉｎｄｅｘ．ｈｔｍｌ 寺本進・藤井孝蔵，１９９８，直交格子法による三次元物体周りの流れ解析，第１２回数値流体力学シンポジウム講演論文集，２９９−３００． Ｑｕｉｒｋ，Ｊ．Ｊ．，１９９４，Ａｎ Ａｌｔｅｒｎａｔｉｖｅ ｔｏ Ｕｎｓｔｒｕｃｔｕｒｅｄ Ｇｒｉｄｓ ｆｏｒ Ｃｏｍｐｕｔｉｎｇ Ｇａｓ Ｄｙｎａｍｉｃ Ｆｌｏｗｓ Ａｒｏｕｎｄ Ａｒｂｉｔｒａｒｉｌｙ Ｃｏｍｐｌｅｘ Ｔｗｏ−Ｄｉｍｅｎｓｉｏｎａｌ Ｂｏｄｉｅｓ，Ｃｏｍｐｕｔｅｒｓ Ｆｌｕｉｄｓ，２３，１２５−１４２． Ｋａｒｍａｎ，Ｓ．Ｌ．Ｊｒ．，１９９５，ＳＰＬＩＴＦＬＯＷ：Ａ ３Ｄ Ｕｎｓｔｒｕｃｔｕｒｅｄ Ｃａｒｔｅｓｉａｎ／Ｐｒｉｓｍａｔｉｃ Ｇｒｉｄ （１２） ｙｎａｍｉｃｓ ｏｆ ＣＦＤ Ｃｏｄｅ ｆｏｒ Ｃｏｍｐｌｅｘ Ｇｅｏｍｅ−ｔｒｉｅｓ，ＡＩＡＡ ９５−０３４３． Ｈｉｒｔ，Ｃ．Ｗ．，＆ Ｎｉｃｈｏｌｓ，Ｂ．Ｄ．，１９８１，Ｖｏｌｕｍｅ ｏｆ Ｆｌｕｉｄ（ＶＯＦ） Ｍｅｔｈｏｄ ｆｏｒ ｔｈｅ Ｄ Ｆｒｅｅ Ｂｏｕｎｄａｒｉｅｓ，Ｊ．Ｃｏｍｐｕｔ．Ｐｈｙｓ．３９，２０１−２２５． Ｈｉｒｔ，Ｃ．Ｗ．，＆ Ｃｏｏｋ，Ｊ．Ｌ．，１９７２，Ｃａｌｃｕｌａｔｉｎｇ Ｔｈｒｅｅ−ｄｉｍｅｎｓｉｏｎａｌ Ｆｌｏｗｓ Ａｒｏｕｎｄ Ｓｔｒｕｃｔｕｒｅｓ ａｎｄ Ｏｖｅｒ Ｒｏｕｇｈ Ｔｅｒｒａｉｎ，Ｊ．Ｃｏｍｐｕｔ．ｐｈｙｓ．１０，３２４−３４０． 加瀬究・手嶋吉法，２００１，ボリュームＣＡＤの開発，理研シンポジウム・ものつくり情報技術統合化研究，第１回，６−１１． 豊田郁夫・荒川忠一，１９９９，直交適合格子法による円柱周り流れの解析，第１３回数値流体力学シンポジウム，Ｆ０３−１，ＣＤ−ＲＯＭ． 松宮輝・木枝香織・谷口伸行・小林敏雄，１９９３，三次精度風上差分法による二次元円柱周り流れの数値シミュレーション，日本機械学会論文集，５９Ｂ−５６６，２９３７−２９４３． Ｂｏｕａｒｄ，Ｒ．，＆ Ｃｏｕｔａｎｃｅａｕ，Ｍ．，１９８０，Ｔｈｅ ｅａｒｌｙ ｓｔａｇｅ ｏｆ ｄｅｖｅｌｏｐｍｅｎｔ ｏｆ ｔｈｅ ｗａｋｅ ｂｅｈｉｎｄ ａｎ ｉｍｐｕｌｓｉｖｅｌｙ ｓｔａｒｔｅｄ ｃｙｌｉｎｄｅｒ ｆｏｒ ４０＜Ｒｅ＜１０＾４，Ｊ．Ｆｌｕｉｄ Ｍｅｃｈ．，１０１−３，５８３−６０７． Ｏｋａｍｏｔｏ，Ｓ．，Ｕｅｍａｔｓｕ，Ｒ．，ａｎｄ Ｔａｇｕｗａ，Ｙ．，Ｆｌｕｉｄ ｆｏｒｃｅ ａｃｔｉｎｇ ｏｎ ｔｗｏ−ｄｉｍｅｎｓｉｏｎａｌ ｃｉｒｃｕｌａｒ ｃｙｌｉｎｄｅｒ ｉｎ Ｌｏｋ−ｉｎ ｐｈｅｎｏｍｅｎｏｎ，ＪＳＭＥ Ｉｎｔｅｒｎａｔｉｏｎａｌ Ｊｏｕｒｎａｌ，Ｂ４５，Ｎｏ．４，（２００２），８５０−８５６． 近藤典夫，円柱の空力挙動に関する数値シミュレーション，第１５回数値流体力学シンポジウム，Ｅ０９−２，（２００１），ＣＤ−ＲＯＭ．
流体解析では現在、複雑な三次元形状の流れ場も重合格子や非構造格子技法を使って計算できるようになったが、メッシュ生成がシミュレーション全体の大きな部分を占めるようになった。このため、完全自動化できるメッシュ生成法として、直交格子を用いることが有望視されている。
直交格子系による任意形状の数値解析では、物体境界の取り扱いの難しさがよく知られている。近年、流れ場の境界近くの離散化や境界条件の扱う方法によって、いくつかの直交格子法が提案されている。
例えば、仮想境界法［非特許文献１］（Ｖｉｒｔｕａｌ Ｂｏｕｎｄａｒｙ）、ＣＩＰ［非特許文献２］（Ｃｕｂｉｃ−Ｉｎｔｅｒｐｏｌａｔｅｄ Ｐｒｏｐａｇａｔｉｏｎ）密度関数法、Ｉｍｍｅｒｓｅｄ Ｂｏｕｎｄａｒｙ ｍｅｔｈｏｄ法［非特許文献３］，ＮＰＬＣ（Ｎｅｉｇｈｂｏｒｉｎｇ Ｐｏｉｎｔ Ｌｏｃａｌ Ｃｏｌｌｏｃａｔｉｏｎ）法［非特許文献４］，格子点間に位置する境界までの距離を差分スキームに取り込む法［非特許文献５］、部分境界適合直交格子法［非特許文献６］などがある。
これらの方法は物体境界が厳密に扱われるが、その分計算処理が複雑であり、必ずしも任意形状の三次元問題に向いているとはいえない。
一方、実用化の観点から、直交格子法は基本的に直交格子だけで階段状の境界を生成し、物体形状を近似する方法（例えば、小野［非特許文献７］・日産自動車、桑原［非特許文献８］・計算流体研）と、カットセルを導入して境界形状を取り扱う近似度を高める方法（例えば、藤井［非特許文献９］・宇宙研、Ｑｕｉｒｋ［非特許文献１０］．Ｊ．Ｊ．，ＮＡＳＡ）の二種類が有望である。
しかし、カットセルによる方法では、直交格子内を境界が任意の場所を通るため、境界上で隣り合ったセルの大きさに大きな差が生じることがあり、カットセル直交格子では粘性流れの解析が難しいという報告［非特許文献１１］もあった。
上述したように、従来の重合格子や非構造格子を用いて、非圧縮粘性流体の流れ場の数値解析を行う場合には、格子生成を完全には自動化できず、そのため格子生成がシミュレーション時間全体に占める割合が高く、シミュレーション時間の短縮が困難な問題点があった。
一方、直交格子を用いた流れ場の数値解析は、格子生成を自動化できるものの、直交格子により物体境界を表現することが難しく、結果としてシミュレーション精度がわるい問題点があった。特に、移動境界に伴う流れの数値計算では、境界の移動距離が一定の大きさのメッシュの倍数に制約されるため、計算不安定に落ちる場合も多かった。
また特に、カットセルによる方法では、直交格子内を境界が任意の場所を通るため、境界上で隣り合ったセルの大きさに大きな差が生じることがあり、カットセル直交格子では粘性流れの解析が難しかった。
発明の要約
本発明は、かかる問題点を解決するために創案されたものである。すなわち、本発明の目的は、格子生成を完全に自動化することができ、物体境界の表現が容易であり、比較的簡単な計算処理により、計算不安定に落ちいることなく、精度の高いシミュレーションを短時間で行うことができるＶ−ＣＡＤデータを直接用いた非圧縮性粘性流体の流れ場の数値解析方法と装置を提供することにある。
本発明によれば、非圧縮性粘性流体と接する対象物の境界データからなる外部データ（１２）を境界が直交する複数のセル（１３）に分割する分割ステップ（Ａ）と、分割された各セルを対象物の内側に位置する内部セル（１３ａ）と境界データを含む境界セル（１３ｂ）とに区分するセル区分ステップ（Ｂ）と、前記境界データによる境界セル（１３ｂ）の稜線の切断点を求める切断点決定ステップ（Ｃ）と、求めた切断点を結ぶ多角形を境界面のセル内部データとする境界面決定ステップ（Ｄ）と、流れ場の境界に、ＶＯＦ法を併用したカットセル有限体積法を適用して解析する解析ステップ（Ｅ）と、を備えたことを特徴とするＶ−ＣＡＤデータを直接用いた非圧縮性粘性流体の流れ場の数値解析方法が提供される。
また、本発明によれば、非圧縮性粘性流体と接する対象物（１）の境界データからなる外部データ（１２）を入力する入力装置（２）と、形状と物理量を統合した実体データとその記憶演算プログラムを記憶する外部記憶装置（３）と、前記記憶プログラムを実行するための内部記憶装置（４）及び中央処理装置（５）と、実行結果を出力する出力装置（６）とを備え、
前記外部データを境界が直交する複数のセル（１３）に分割し、分割された各セルを対象物の内側に位置する内部セル（１３ａ）と境界データを含む境界セル（１３ｂ）とに区分し、前記境界データによる境界セル（１３ｂ）の稜線の切断点を求め、求めた切断点を結ぶ多角形を境界面のセル内部データとし、流れ場の境界に、ＶＯＦ法を併用したカットセル有限体積法を適用して解析する、ことを特徴とするＶ−ＣＡＤデータを直接用いた非圧縮性粘性流体の流れ場の数値解析装置が提供される。
さらに、本発明によれば、コンピュータに、非圧縮性粘性流体と接する対象物の境界データからなる外部データ（１２）を境界が直交する複数のセル（１３）に分割する分割ステップ（Ａ）と、分割された各セルを対象物の内側に位置する内部セル（１３ａ）と境界データを含む境界セル（１３ｂ）とに区分するセル区分ステップ（Ｂ）と、前記境界データによる境界セル（１３ｂ）の稜線の切断点を求める切断点決定ステップ（Ｃ）と、求めた切断点を結ぶ多角形を境界面のセル内部データとする境界面決定ステップ（Ｄ）と、流れ場の境界に、ＶＯＦ法を併用したカットセル有限体積法を適用して解析する解析ステップ（Ｅ）と、を実行させるための非圧縮性粘性流体の流れ場の数値解析プログラムが提供される。
上記本発明の方法及び装置により、精度と安定性を有し、計算コストがあまりかからず、かつ、格子生成を完全に自動化することができ、物体境界の表現が容易であり、比較的簡単な計算処理により、精度の高いシミュレーションを短時間で行うことができることが分かった。
本発明の好ましい実施形態によれば、前記解析ステップ（Ｅ）において、空間積分について、対流項に二次元のＱＵＩＣＫ補間スキームを使用し、拡散項に二次精度の中心差分を用い、時間進行について、対流項と拡散項を合わせて二次精度のＡｄａｍｓ−Ｂａｓｈｆｏｒｔｈ法を適用し、圧力勾配項に一次精度のＥｕｌｅｒ陰解法を用いる。
この方法により、格子生成を完全に自動化することができ、かつ安定性と離散化精度を確保することができる。
二次元境界セルにおいて、有限体積法における支配方程式を式（７）であらわす。
この式（７）は、非圧縮粘性流体の基礎支配方程式（１）をテンソル形の発散型に書き換えた式（６）を二次元境界セルの流体部分をコントロールボリューム（ＣＶ）Ｖｉ，ｊとして、空間積分したものであり、非圧縮粘性流体の基礎支配方程式（１）を満たすことができる。
有限体積法における支配方程式における対流項、圧力勾配項および拡散項を式（８）、（９）、（１０）であらわす。
式（８）〜（１０）には、稜線の切断点を結ぶ対象物の境界データが含まれるため、境界における流れ場の非圧縮性粘性流体の数値解析ができる。
固体境界の積分において、固体境界でノースリップ境界条件の場合に、対流項は零とし、圧力勾配項と拡散項に対しては、切断線の中間点Ｐの値を平均値として用いて積分し、空間積分に対してはすべての項に開口率を適用する。
この方法により、直交格子の完全自動化が保たれ、かつ流体計算におけるコントロールボリュームで保存則を厳密に満たすことができる。
カットセルにおいて算出される物理変数の定義点を、ＶＯＦ＝０．０１を閾値としてそれより小さい境界セルは完全な固体とみなし、それより大きい境界セルにおいて算出される変数をセルの中心に置く、また、稜線の変数定義点をセル稜線の中心に定義し、更に、線分４の中心点の変数値を線形補間により求める。
この方法によっても、直交格子の完全自動化が保たれ、かつ流体計算におけるコントロールボリュームで保存則を厳密に満たすことができる。
物体に働く抗力（流れ方向の力）と揚力（流れの垂直方向の力）を、式（１２）（１３）であらわす。
この式により、直交格子において、抗力と揚力を容易かつ精度よく求めるもとができる。
移動境界を伴う流体構造連成解析において、所定の時間刻み毎に、流体系と構造系を別々に解析し、それぞれの境界条件を陽的に受け渡す、ことが好ましい。
この解析方法により、流体系と構造系を連立して統一的に扱う強連成に比較して計算処理を簡略化でき、計算時間を短縮化できる。また境界の移動距離が一定の大きさのメッシュの倍数に制約されないので、計算不安定に落ちいることがほとんどない。
円柱の強制振動の流体構造連成解析において、円柱が弾性支持された剛体構造物として、流れ方向に対して直角方向に振動する１質点１自由度系モデルと仮定し、円柱中心のＹ方向の変位を式（１７）で与え、円柱表面のＹ方向の速度境界条件を、式（１８）で与える。
また、式（１８）から求められた振動する円柱の運動速度を、流れ場の速度境界条件として計算時間ステップごとに変化して与える。
この解析方法により、（１）振動する円柱に作用する抗力および揚力が静止円柱のものより大きい、（２）Ｓｔｒｏｕｈａｌ数０．２付近の０．１５＜ｆｃ＜０．２５において、強制振動による渦放出周波数のロックイン現象がよく再現される、（３）ロックイン領域での抗力や揚力が増大する、等、実験データと整合性の高い解析結果が確認された。
また円柱の渦放出による自励振動の流体構造連成解析において、、有次元の振動方程式を１質点１自由度のダンパ・ばねモデルとて、式（１９）又は式（２０）であらわす。
さらに式（２０）から求められた振動する円柱の運動速度を、流れ場の速度境界条件として計算時間ステップごとに変化して与える。
また円柱の初期変位と初期速度を０とし、かつ式（２０）の揚力は現在の値を用いて陽的に与え、式（２０）の振動方程式をニューマークのβ法により積分し、円柱の振動変位と振動速度を求める。
この解析方法により、（１）無次元流速Ｖｒ＝４（円柱の固有振動数ｆｏ＝１／４＝０．２５）の場合に円柱の振幅がもっとも大きくなる、（２）最大変位ｙが、Ｖｒ＝２、３、４ではｙの最大変位は次第に大きくなり、Ｖｒ＝４で最大になる、（３）その後、Ｖｒが大きくなるにつれて、ｙの最大変位は波状に小さくなっていく、等、ＡＬＥ有限要素法による計算結果と定性的に一致していることが確認された。
本発明のその他の目的及び有利な特徴は、添付図面を参照した以下の説明から明らかになろう。
直交格子系による任意形状の数値解析では、物体境界の取り扱いの難しさがよく知られている。近年、流れ場の境界近くの離散化や境界条件の扱う方法によって、いくつかの直交格子法が提案されている。
例えば、仮想境界法［非特許文献１］（Ｖｉｒｔｕａｌ Ｂｏｕｎｄａｒｙ）、ＣＩＰ［非特許文献２］（Ｃｕｂｉｃ−Ｉｎｔｅｒｐｏｌａｔｅｄ Ｐｒｏｐａｇａｔｉｏｎ）密度関数法、Ｉｍｍｅｒｓｅｄ Ｂｏｕｎｄａｒｙ ｍｅｔｈｏｄ法［非特許文献３］，ＮＰＬＣ（Ｎｅｉｇｈｂｏｒｉｎｇ Ｐｏｉｎｔ Ｌｏｃａｌ Ｃｏｌｌｏｃａｔｉｏｎ）法［非特許文献４］，格子点間に位置する境界までの距離を差分スキームに取り込む法［非特許文献５］、部分境界適合直交格子法［非特許文献６］などがある。
これらの方法は物体境界が厳密に扱われるが、その分計算処理が複雑であり、必ずしも任意形状の三次元問題に向いているとはいえない。
一方、実用化の観点から、直交格子法は基本的に直交格子だけで階段状の境界を生成し、物体形状を近似する方法（例えば、小野［非特許文献７］・日産自動車、桑原［非特許文献８］・計算流体研）と、カットセルを導入して境界形状を取り扱う近似度を高める方法（例えば、藤井［非特許文献９］・宇宙研、Ｑｕｉｒｋ［非特許文献１０］．Ｊ．Ｊ．，ＮＡＳＡ）の二種類が有望である。
しかし、カットセルによる方法では、直交格子内を境界が任意の場所を通るため、境界上で隣り合ったセルの大きさに大きな差が生じることがあり、カットセル直交格子では粘性流れの解析が難しいという報告［非特許文献１１］もあった。
上述したように、従来の重合格子や非構造格子を用いて、非圧縮粘性流体の流れ場の数値解析を行う場合には、格子生成を完全には自動化できず、そのため格子生成がシミュレーション時間全体に占める割合が高く、シミュレーション時間の短縮が困難な問題点があった。
一方、直交格子を用いた流れ場の数値解析は、格子生成を自動化できるものの、直交格子により物体境界を表現することが難しく、結果としてシミュレーション精度がわるい問題点があった。特に、移動境界に伴う流れの数値計算では、境界の移動距離が一定の大きさのメッシュの倍数に制約されるため、計算不安定に落ちる場合も多かった。
また特に、カットセルによる方法では、直交格子内を境界が任意の場所を通るため、境界上で隣り合ったセルの大きさに大きな差が生じることがあり、カットセル直交格子では粘性流れの解析が難しかった。
発明の要約
本発明は、かかる問題点を解決するために創案されたものである。すなわち、本発明の目的は、格子生成を完全に自動化することができ、物体境界の表現が容易であり、比較的簡単な計算処理により、計算不安定に落ちいることなく、精度の高いシミュレーションを短時間で行うことができるＶ−ＣＡＤデータを直接用いた非圧縮性粘性流体の流れ場の数値解析方法と装置を提供することにある。
本発明によれば、非圧縮性粘性流体と接する対象物の境界データからなる外部データ（１２）を境界が直交する複数のセル（１３）に分割する分割ステップ（Ａ）と、分割された各セルを対象物の内側に位置する内部セル（１３ａ）と境界データを含む境界セル（１３ｂ）とに区分するセル区分ステップ（Ｂ）と、前記境界データによる境界セル（１３ｂ）の稜線の切断点を求める切断点決定ステップ（Ｃ）と、求めた切断点を結ぶ多角形を境界面のセル内部データとする境界面決定ステップ（Ｄ）と、流れ場の境界に、ＶＯＦ法を併用したカットセル有限体積法を適用して解析する解析ステップ（Ｅ）と、を備えたことを特徴とするＶ−ＣＡＤデータを直接用いた非圧縮性粘性流体の流れ場の数値解析方法が提供される。
また、本発明によれば、非圧縮性粘性流体と接する対象物（１）の境界データからなる外部データ（１２）を入力する入力装置（２）と、形状と物理量を統合した実体データとその記憶演算プログラムを記憶する外部記憶装置（３）と、前記記憶プログラムを実行するための内部記憶装置（４）及び中央処理装置（５）と、実行結果を出力する出力装置（６）とを備え、
前記外部データを境界が直交する複数のセル（１３）に分割し、分割された各セルを対象物の内側に位置する内部セル（１３ａ）と境界データを含む境界セル（１３ｂ）とに区分し、前記境界データによる境界セル（１３ｂ）の稜線の切断点を求め、求めた切断点を結ぶ多角形を境界面のセル内部データとし、流れ場の境界に、ＶＯＦ法を併用したカットセル有限体積法を適用して解析する、ことを特徴とするＶ−ＣＡＤデータを直接用いた非圧縮性粘性流体の流れ場の数値解析装置が提供される。
さらに、本発明によれば、コンピュータに、非圧縮性粘性流体と接する対象物の境界データからなる外部データ（１２）を境界が直交する複数のセル（１３）に分割する分割ステップ（Ａ）と、分割された各セルを対象物の内側に位置する内部セル（１３ａ）と境界データを含む境界セル（１３ｂ）とに区分するセル区分ステップ（Ｂ）と、前記境界データによる境界セル（１３ｂ）の稜線の切断点を求める切断点決定ステップ（Ｃ）と、求めた切断点を結ぶ多角形を境界面のセル内部データとする境界面決定ステップ（Ｄ）と、流れ場の境界に、ＶＯＦ法を併用したカットセル有限体積法を適用して解析する解析ステップ（Ｅ）と、を実行させるための非圧縮性粘性流体の流れ場の数値解析プログラムが提供される。
上記本発明の方法及び装置により、精度と安定性を有し、計算コストがあまりかからず、かつ、格子生成を完全に自動化することができ、物体境界の表現が容易であり、比較的簡単な計算処理により、精度の高いシミュレーションを短時間で行うことができることが分かった。
本発明の好ましい実施形態によれば、前記解析ステップ（Ｅ）において、空間積分について、対流項に二次元のＱＵＩＣＫ補間スキームを使用し、拡散項に二次精度の中心差分を用い、時間進行について、対流項と拡散項を合わせて二次精度のＡｄａｍｓ−Ｂａｓｈｆｏｒｔｈ法を適用し、圧力勾配項に一次精度のＥｕｌｅｒ陰解法を用いる。
この方法により、格子生成を完全に自動化することができ、かつ安定性と離散化精度を確保することができる。
二次元境界セルにおいて、有限体積法における支配方程式を式（７）であらわす。
有限体積法における支配方程式における対流項、圧力勾配項および拡散項を式（８）、（９）、（１０）であらわす。
固体境界の積分において、固体境界でノースリップ境界条件の場合に、対流項は零とし、圧力勾配項と拡散項に対しては、切断線の中間点Ｐの値を平均値として用いて積分し、空間積分に対してはすべての項に開口率を適用する。
この方法により、直交格子の完全自動化が保たれ、かつ流体計算におけるコントロールボリュームで保存則を厳密に満たすことができる。
カットセルにおいて算出される物理変数の定義点を、ＶＯＦ＝０．０１を閾値としてそれより小さい境界セルは完全な固体とみなし、それより大きい境界セルにおいて算出される変数をセルの中心に置く、また、稜線の変数定義点をセル稜線の中心に定義し、更に、線分４の中心点の変数値を線形補間により求める。
この方法によっても、直交格子の完全自動化が保たれ、かつ流体計算におけるコントロールボリュームで保存則を厳密に満たすことができる。
物体に働く抗力（流れ方向の力）と揚力（流れの垂直方向の力）を、式（１２）（１３）であらわす。
移動境界を伴う流体構造連成解析において、所定の時間刻み毎に、流体系と構造系を別々に解析し、それぞれの境界条件を陽的に受け渡す、ことが好ましい。
この解析方法により、流体系と構造系を連立して統一的に扱う強連成に比較して計算処理を簡略化でき、計算時間を短縮化できる。また境界の移動距離が一定の大きさのメッシュの倍数に制約されないので、計算不安定に落ちいることがほとんどない。
円柱の強制振動の流体構造連成解析において、円柱が弾性支持された剛体構造物として、流れ方向に対して直角方向に振動する１質点１自由度系モデルと仮定し、円柱中心のＹ方向の変位を式（１７）で与え、円柱表面のＹ方向の速度境界条件を、式（１８）で与える。
この解析方法により、（１）振動する円柱に作用する抗力および揚力が静止円柱のものより大きい、（２）Ｓｔｒｏｕｈａｌ数０．２付近の０．１５＜ｆｃ＜０．２５において、強制振動による渦放出周波数のロックイン現象がよく再現される、（３）ロックイン領域での抗力や揚力が増大する、等、実験データと整合性の高い解析結果が確認された。
また円柱の渦放出による自励振動の流体構造連成解析において、、有次元の振動方程式を１質点１自由度のダンパ・ばねモデルとて、式（１９）又は式（２０）であらわす。
また円柱の初期変位と初期速度を０とし、かつ式（２０）の揚力は現在の値を用いて陽的に与え、式（２０）の振動方程式をニューマークのβ法により積分し、円柱の振動変位と振動速度を求める。
この解析方法により、（１）無次元流速Ｖｒ＝４（円柱の固有振動数ｆｏ＝１／４＝０．２５）の場合に円柱の振幅がもっとも大きくなる、（２）最大変位ｙが、Ｖｒ＝２、３、４ではｙの最大変位は次第に大きくなり、Ｖｒ＝４で最大になる、（３）その後、Ｖｒが大きくなるにつれて、ｙの最大変位は波状に小さくなっていく、等、ＡＬＥ有限要素法による計算結果と定性的に一致していることが確認された。
本発明のその他の目的及び有利な特徴は、添付図面を参照した以下の説明から明らかになろう。
図１は、先行出願の実体データの記憶方法のフロー図である。
図２は、本発明の数値解析方法を実行するための数値解析装置の構成図である。
図３は、本発明の数値解析方法とそのプログラムのフロー図である。
図４Ａと図４Ｂは、二次元における開口率の定義を示す図である。
図５は、境界セルにおける有限体積法による空間離散化の説明図である。
図６Ａ，Ｂ，Ｃは、解析した３ケースのＶＯＦ分布を示す画像である。
図７Ａ，Ｂ，Ｃは、解析した３ケースの流れ方向速度の分布と理論解析解との比較図である。
図８は、解析格子とＶＯＦ分布を示す画像である。
図９は、速度ベクトルを示す画像である。
図１０は、圧力等値線を示す画像である。
図１１は、本発明による平均抗力係数を示す図である。
図１２は、本発明による平均Ｓｔｒｏｕｈａｌ数を示す図である。
図１３Ａは、直交適合格子法による静止円柱の抗力係数を示す図であり、図１３Ｂは、その平均Ｓｔｒｏｕｈａｌ数を示す図である。
図１４Ａは、一般座標系境界適合格子による静止円柱の抗力係数を示す図であり、図１４Ｂはその平均Ｓｔｒｏｕｈａｌ数を示す図である。
図１５Ａは、Ｒｅ＝３００，Ｔ＝２．５の場合の実験による流脈線を示す画像であり、図１５Ｂはその流線を示す画像である。
図１６Ａは、Ｒｅ＝５５０，Ｔ＝２．５の場合の実験による流脈線を示す画像であり、図１６Ｂは、その流線を示す画像である。
図１７Ａは、抗力と揚力の履歴変化を示す本発明の結果であり、図１７Ｂは、その部分境界適合直交格子法の結果である。
図１８Ａ〜Ｆは、実用問題に向けての応用例を示す画像である。
図１９Ａ〜Ｉは、静止円柱と各振動数での円柱に作用する流体力の抗力係数ＣＤおよび揚力係数ＣＬの時間履歴の変動波形を示す図である。
図２０は、流れの渦放出状態を表すＳｔｒｏｕｈａｌ数Ｓｔと強制振動周波数ｆｃの関係図である。
図２１Ａ〜Ｃは、無次元流速Ｖｒが２，３，４における変位ｙ（左図）と抗力係数ＣＤ及び揚力係数ＣＬ（右図）の時間履歴を示す図である。
図２２Ｄ〜Ｇは、無次元流速Ｖｒが５，６，７，８における変位ｙ（左図）と抗力係数ＣＤ及び揚力係数ＣＬ（右図）の時間履歴を示す図である。
図２３は、各々の無次元流速Ｖｒに対するｙの最大変位を示す図である。
好ましい実施例の説明
以下、本発明の好ましい実施形態を図面を参照して説明する。
図２は、本発明の数値解析方法を実行するための数値解析装置の構成図である。
この図に示すように、本発明の数値解析装置１０は、入力装置２、外部記憶装置３、内部記憶装置４、中央処理装置５および出力装置６を備える。
入力装置２は、例えばキーボードであり、対象物１の形状データからなる外部データ１２を入力する。外部記憶装置３は、ハードディスク、フロッピィーディスク、磁気テープ、コンパクトディスク等であり、形状と物理量を統合した実体データとその記憶演算プログラムを記憶する。内部記憶装置４は、例えばＲＡＭ，ＲＯＭ等であり、演算情報を保管する。中央処理装置５（ＣＰＵ）は、演算や入出力等を集中的に処理し、内部記憶装置４と共に、記憶プログラムを実行する。出力装置６は、例えば表示装置とプリンタであり、記憶した実体データと記憶プログラムの実行結果を出力するようになっている。
本発明の記憶演算装置１０は、上述した外部記憶装置３、内部記憶装置４、及び中央処理装置５により、外部データを境界が直交する複数のセル１３に分割し、分割された各セルを対象物の内側に位置する内部セル１３ａと境界データを含む境界セル１３ｂとに区分し、境界データによる境界セル１３ｂの稜線の切断点を求め、求めた切断点を結ぶ多角形を境界面のセル内部データとし、流れ場の境界に、ＶＯＦ法を併用したカットセル有限体積法を適用して解析する。
図３は、本発明の数値解析方法とそのプログラムのフロー図である。この図に示すように、本発明の方法及び変換プログラムは、分割ステップ（Ａ）、セル区分ステップ（Ｂ）、切断点決定ステップ（Ｃ）、境界面決定ステップ（Ｄ）、及び解析ステップ（Ｅ）からなる。
外部から入力する外部データ１２は、多面体を表すポリゴンデータ、有限要素法に用いる四面体又は六面体要素、３次元ＣＡＤ又はＣＧツールに用いる曲面データ、或いはその他の立体の表面を部分的な平面や曲面で構成された情報で表現するデータである。
外部データ１２は、このようなデータ（Ｓ−ＣＡＤデータと呼ぶ）のほかに、（１）Ｖ−ＣＡＤ独自のインターフェース（Ｖ−ｉｎｔｅｒｆａｃｅ）により人間の入力により直接作成されたデータと、（２）測定機やセンサ、デジタイザなどの表面のデジタイズデータや、（３）ＣＴスキャンやＭＲＩ、および一般的にＶｏｌｕｍｅレンダリングに用いられているボクセルデータなどの内部情報ももつＶｏｌｕｍｅデータであってもよい。
分割ステップ（Ａ）では、外部データ取得ステップ（図示せず）で取得した非圧縮性粘性流体と接する対象物の境界データからなる外部データ１２を境界平面が直交する複数のセル１３に分割する。この分割は、３次元の場合は、八分木分割であり、２次元の場合には４分割である。
すなわちこの分割ステップ（Ａ）における分割とは、目的とする非圧縮性粘性流体と接する対象物を含む、基準となる直方体（または矩形）を分割（８分割又は４分割）し、それぞれの領域の中に立体が完全に含まれるか、含まれなくなるまで再帰的に分割処理を繰り返す。この分割によりボクセル表現よりも大幅にデータ量を減らすことができる。
空間分割により分割された一つの空間領域をセル１３とよぶ。セルは境界が直交する直方体または矩形である。セルによる階層構造、分割数もしくは分解能によって空間中に占める領域を表現する。これにより空間全体の中で対象は大きさの異なるセルを積み重ねたものとして表現される。
セル区分ステップ（Ｂ）では、分割された各セルを対象物の内側に位置する内部セル１３ａと境界データを含む境界セル１３ｂとに区分する。
すなわち本発明では非圧縮性粘性流体と接する対象物の内側または外側に完全に内部に含まれるものはその最大の大きさをもつ内部セル１３ａ（立方体）とし、外部データ１２からの境界情報を含むセルは境界セル１３ｂとする。
切断点決定ステップ（Ｃ）では、境界データによる境界セル１３ｂの稜線の切断点を求める。
境界面決定ステップ（Ｄ）では、求めた切断点を結ぶ多角形を境界面のセル内部データとする。以下、かかる切断点を結ぶ多角形を含むセルを「カットセル」と呼ぶ。
解析ステップ（Ｅ）では、上述した内部セル１３ａと境界セル１３ｂに対して、流れ場の境界に、ＶＯＦ法を併用したカットセル有限体積法を適用して非圧縮性粘性流体の流れ場の数値解析を行う。
以下、本発明を更に詳細に説明する。
１．本発明ではＶ−ＣＡＤプログラムにおける実用性のある流体解析技術を目指し、カットセル（ＫＴＣ）直交格子による任意形状の非圧縮粘性流れの解析方法を開発した。本発明において、流れ場の境界での扱いについては、Ｈｉｒｔらが提案しているＶＯＦ（Ｖｏｌｕｍｅ Ｏｆ Ｆｌｕｉｄ）法［非特許文献１２］を併用したカットセル有限体積法を用いた。また本発明の方法による内部流であるチャンネル内流れと、外部流の円柱まわり流れを解析し、実験データ、理論解、及び既存の方法による解析結果との比較を行った。更に、本解析法の応用例として、いくつかの計算結果を合わせて示した。
２．基礎方程式と計算方法
２．１ 基礎方程式
本発明で用いる基礎支配方程式は［数４］の式（１）（２）に示す非圧縮粘性流体のＮａｖｉｅｒ・Ｓｔｏｋｅｓ方程式と連続の式である。
ここでＲｅは流れ場の代表長さと代表速度で定義されるレイノルズ数と呼ばれる無次元数であり、物理的には流れ場における慣性力と粘性力の比を表している。ｕは速度、ｐは圧力である。また、ｉ＝１，２，３、ｊ＝１，２，３は直交座標系での各方向を表している。なお、ｉ，ｊに関しては縮約を取ることとして、これ以降もそのように表すものとする。式（１）より流れ場における外力ｆｉを考慮しない場合に非圧縮性粘性流体の運動はただ一つのパラメータＲｅで支配されていることが分かる。
２．２ 計算精度・計算コスト及び計算方法
流体のシミュレーションは一種の数値実験であり、つねに一定の誤差を伴う。流体解析を精度よく行うには、以下の四つの条件が必要となる。
▲１▼流れの最小長さスケール（境界層・渦・衝撃波・火炎面など）を捕られるほど細かい空間解像度。
▲２▼流れの最大長さスケールを十分に捕獲するほど大きい計算領域と、人為的な流入・流出境界条件・壁によるブロッキングなど影響を無視できるほどの計算領域。
▲３▼打ち切り誤差や数値拡散などを無視できるほど十分な空間的・時間的な離散化精度。
▲４▼特定問題に即したモデル（壁面モデル・乱流ＳＧＳモデル・燃焼モデルなど）やスキーム（Ｋ−Ｋ・ＱＵＩＣＫ・ＭＵＳＣＬ・ＴＶＤ・ＥＮＯなど）、差分法を援用すれば、できるだけ直交等方性のある格子。
流体解析の計算コスト（計算時間・必要なメモリ）は、その計算方法に求められる解析精度によって決められるといえる。実用的な流体解析では、ある程度の精度（例えば、普通の物理実験の精度範囲や、ユーザーに求められる精度）があれば、十分と考えられる。
一方、流れ場は内部流、外部流及びそのほか（例えば噴流）に分類される。一般に外部流の解析領域は内部流のものより大きくとるのが普通である。本発明では実用性の観点から、外部流の場合には、流れ場の代表長さスケールＤ（＝１）として、解析領域の大きさはすべて１０Ｄ×１０Ｄとした。
本発明では、差分法を併用した有限体積法を用いる。ある程度の安定性や離散化精度を確保するため、空間積分について、対流項に二次元のＱＵＩＣＫ補間スキームを使用し、拡散項に二次精度の中心差分を用いる。時間進行について、対流項と拡散項を合わせて二次精度のＡｄａｍｓ−Ｂａｓｈｆｏｒｔｈ法を適用し、圧力勾配項に一次精度のＥｕｌｅｒ陰解法を用いる。解析アルゴリズムとしての圧力と速度のカップリングには、Ｈｉｒｔらにより提案された圧力と速度を同時緩和するＳＯＬＡ−ＨＳＭＡＣ法［非特許文献１３］を使用し、緩和係数１．６５でＳＯＲ法を用いて反復計算を行う。収束判定は連続方程式（２）の残差が０．０００２とする。また圧力振動を防ぐために三方向の速度ｕ，ｖ，ｗ，及び圧力ｐの格子上の定義点を、半メッシュをずれるスタガード格子を配置する。
２．３ Ｖ−ＣＡＤデータによる固体・流体境界の取扱い
Ｖｏｌｕｍｅ−ＣＡＤシステムの中でもっとも重要なアルゴリズムは加瀬・手嶋［非特許文献１４］により提案されるセルの切断点からモノの表面形状を復元するＫＴＣアルゴリズムである。このＫＴＣアルゴリズムは流体解析分野における直交格子カットセル（Ｃｕｔ Ｃｅｌｌ）法と同じ概念であることを注意されたい。また、ＫＴＣは二次元においても一個のセルの四つの稜線上に２個以上の交点（切断点）があり得るが、本発明では簡単のために二次元で一個のセルの稜線上に２個だけの切断点に限るとする。この場合のＫＴＣは図４Ａと図４Ｂに示すように二種類しかないことがわかる。
カットセル法における開口率は、セルの稜線に流体部分が占める比率として図４Ａと図４Ｂに示すように流体・固体境界とセル稜線との交点情報から［数５］の式（３）で定義される。
ここで、Δｘｉ，ｊ、Δｙｉ，ｊはそれぞれＸ方向とＹ方向の格子幅であり、Ａｉ， ｊ、Ｂｉ，ｊはＸ方向とＹ方向の開口率である。また、Ｖｘｉ，ｊ、Ｖｙｉ，ｊはそのセル稜線上に固体が 占める線分を表す。
各方向の開口率が分かれば、ＶＯＦ法における流体の体積占有率は［数５］の式（４）と式（５）により求められる。
２．４ ＶＯＦ有限体積法による離散化
有限体積法で保存型の支配方程式を離散化する場合、方程式（１）の積分形をグリーン定理（発散定理）に適用するため、微分型の方程式（１）をテンソル形の発散型に書き換えると［数６］の式（６）のようになる。
ここで、○に×はテンソル積であり、ベクトルＩはＫｒｏｎｅｃｋｅｒδｉｊに対応する単位テンソルである。
流体固体境界での離散化は、本発明の肝心なところになる。境界セル以外の格子における有限体積法による直交格子での離散化は従来の方法でもあまり問題が生じないのでここでは省略する。流体固体境界では図４のようなＫＴＣセルを想定して式（６）の空間離散化を行う。また式（６）の時間離散化について、流体占有率ＶＯＦを用いてコントロールボリューム（ＣＶ）内における速度の時間勾配を平均値とみなして、従来の方法通りで行えるため、境界セルでの時間離散化の説明も省く。
さて、図５の二次元境界セルの流体部分をコントロールボリューム（ＣＶ）Ｖｉ，ｊとして、式（６）を空間積分する（ここで簡単のため外力項ベクトルｆ＝０とした）と、有限体積法における支配方程式は［数６］の式（７）になる。
式（７）における対流項、圧力勾配項および拡散項をそれぞれグリーン定理（発散定理）に適用すると、コントロールボリュームにおける面積分（三次元の場合に体積積分）は、閉直線の線積分（三次元の場合に閉平面の面積分）に転換され、図５のコントロールボリュームの五つの稜線（１→２→３→４→５→１）に（反時計まわり方向をとる）沿って線積分することになる。本発明のデカルト座標では、上記の対流項、圧力勾配項および拡散項はそれぞれ［数７］のように離散化される。
ここで、速度変数の下添え字はスタガード格子ステンシルにより変化することを注意されたい。またはベクトルｉ，ｊはそれぞれＸ方向とＹ方向の方程式を表す。対流項の離散式（８）の中に上添え字（ｘ）（ｙ）が付いた速度は、それぞれＸ方向とＹ方向に二次元ＱＵＩＣＫ補間を適用する。拡散項の離散式（１０）の中に上添え字ｘ、ｙが付いた速度勾配は、それぞれＸ方向とＹ方向の速度勾配を表し、そこで各々二次中心差分を適用する。
カットセルによる離散化には、二つのキーポイントがあることを注意されたい。
キーポイントの一つは固体境界（図５の切断線４）の積分の取り扱いである。有限体積法では、一階偏微分の対流項と圧力勾配項の境界条件はディリクレ条件であり、二階偏微分の拡散項の境界条件はノイマン問題となる。固体境界でノースリップ境界条件の場合に、図５の切断線４のところに対流項が零となるから対流項を積分しなくてもいいが、圧力勾配項と拡散項に対して切断線４のところに積分しなければいけない。これは式（９）と式（１０）の各々最後の項に切断線４の中間点Ｐの値を（平均値として）用いて積分するわけである。また空間積分に対してすべての項に開口率を適用する必要もあることを注意されたい。
もう一つのキーポイントは、カットセルにおいて算出される物理変数の定義点である。有限体積法では流体占有率ＶＯＦで体積平均を行うため、物理変数の定義点は流体のコントロールボリュームの中心（例えば図５の○に・点で示すＣ点）に置くのが自然だが、そうすると直交等間隔格子の完全自動化メリットがなくなる。その故、本発明ではＶＯＦ＝０．０１を閾値としてそれより小さい境界セルは完全の固体とみなし、それより大きい境界セルにおいて算出される変数はコントロールボリュームの中心ではなく、セルの中心（例えば図５の○に＋印で示すＤ点）に置くこととする。同じく図５の稜線３と稜線５の変数定義点も線分の中心○ではなく、セル稜線の中心◎に定義される。また線分４の中心点Ｐに変数が定義されていないため、本発明でその点の変数値を線形補間により求める。こう処理すると、固体境界の再現性を損するものの、直交格子の完全自動化が保たれ、実用化に向ける。また流体の計算はコントロールボリュームで保存則が厳密に満たされ、上記の処理は流れ場の計算精度に及ぼす影響がそれほど大きくないと考えられる。
２．５ 物体に働く抗力と揚力について
物体まわりの流れを解析するとき、その物体に働く抗力と揚力の計算について、一般座標系の境界適合格子の場合に物体の第１格子点に沿って積分すれば簡単だが、直交格子で取り扱う場合にはやや煩雑となる。そこで本発明では、以下の方法で抗力と揚力を求める。
物体に働く力は、その物体表面に作用する流体応力テンソルを閉積分して、グリーン定理を適用すれば、二次元の場合に［数８］の式（１１）ように得られる。
ここでσｎは法線方向の平面に働く応力テンソル、σは物体表面の点での応力テンソルである。また二次元の二階応力テンソルの発散はそのテンソルの四つの方向からなることを注意されたい。それでＸ方向を流れ方向、Ｙ方向を流れの垂直方向とすると、物体に働く抗力（流れ方向の力）と揚力（流れの垂直方向の力）は、グリーン定理を再び適用すると［数８］の式（１２）（１３）のようになる。
抗力と揚力が分かれば、無次元量としての抗力係数ＣＤと揚力係数ＣＬを、主流速度Ｕ、流体密度ρ及び物体の基準面積Ａ（二次元の場合は基準長さ）を用いて［数５］の式（１４）のように無次元化される。
また非定常現象を表わす無次元量のストローハル数Ｓｔは、ＣＬの周波数ｆＬから［数８］の式（１５）で求める。
３．計算結果および考察
流体の数値解析において解析スキームや数値モデルなどを検証することにあたって、内部流のベンチマークテストとしてチャネル流れと、外部流のベンチマークテストとしての円柱まわりの流れはしばしば検証の対象となる。本発明でも、上記のＶＯＦ法を併用したカットセル有限体積法を検証するために、解析理論解のある二次元Ｐｏｉｓｅｕｉｌｌｅ流と、多くの実験や計算データのある静止円柱まわりの流れについて数値解析を行い、他の研究者のデータとの比較を行った。なお、実用問題に向けての応用例として、バックステップを過ぎる流れ、狭窄管内の流れ、チャネル内障害物周りの流れ、分岐管内の流れ、複数Ｂｌｕｆｆ Ｂｏｄｙまわりの流れや予混合燃焼器内の流れなどについて解析を行い、その計算結果の可視化例を示す。
３．１ チャネルＰｏｉｓｅｕｉｌｌｅ流による検証
セルのカットされたＶＯＦの効果を検証するため、１０Ｄ×１０Ｄ正方形の解析領域におけるダクトを、０度、１０度及び４５度を斜めにして解析を行った。全体の格子数は６４×６４、ダクトの半幅は２．５Ｄ、この半幅で定義されるレイノルズ数Ｒｅ＝１、時間刻み△ｔ＝０．０００１、解析時間ｓｔｅｐｓ＝１００００とする。ダクトの流入条件と流出条件はともに放物型の分布を仮定した。流れ方向の圧力勾配は１．０とした。またＨＳＭＡＣ法を用い、入口・出口および壁での圧力境界条件は直接与えていない。三ケースのＶＯＦ分布を図６Ａ，Ｂ，Ｃに示す。図６Ａ，Ｂ，Ｃにおける破線は、境界を示している。その三ケースに対応する流れ方向速度分布を理論解析解と比較して図７Ａ，Ｂ，Ｃに示す。
図７Ａ，Ｂ，Ｃの中に線Ａのａｎａｌｙｓｉｓは理論解析解であり、線Ｂのｉ＝３２は第３２格子点断面にある流れ軸中心での流れ垂直断面の速度分布である。ダクトを傾けない（角度０）計算結果の速度分布は解析解とよく一致し、傾斜角度を１０度、４５度にすると、軸中心の最大速度がわずかに過小評価されるものの、全体として解析結果は理論解と一致していると言える。
３．２ 静止円柱まわり流れによる検証
解析スキームやモデルなどの検証対象となる円柱まわり流れの解析は、直交格子にとって最も困難な例題であろう。本発明では、３．１の解析領域（１０Ｄ×１０Ｄ）と格子数（６４×６４）のままで、円柱直径Ｄで定義されるレイノルズ数のＲｅ＝１からＲｅ＝２００００までの２４ケースについて、一様流中に置かれている静止円柱まわり流れの解析を行った。計算の無次元時間は１００とし、円柱の円心座標の図６のようにＸ＝２．５Ｄ，Ｙ＝５．０Ｄとした。速度の境界条件として、円柱表面で滑り無し条件、流入側の流入条件は一様流（流れ方向速度＝１．０）条件、流出側の流出条件はＸ方向に自由流出（速度勾配０），Ｙ方向にフリースリップ条件とした。また、ＨＳＭＡＣ法を用いて速度・圧力を同時に緩和するため、流入・流出および壁での圧力境界条件を直接与える必要がない。
図８は解析格子と解析全領域のＶＯＦ分布を示す。図９、図１０はレイノルズ数Ｒｅ＝３００の場合の、速度ベクトルと圧力等値線であり、円柱後方に周期的なカルマン渦列の形成が観察されており、非定常流の流れパタンーの特徴を捉えている。これはレイノルズ数Ｒｅ＞８０の場合に共通なので、他のレイノルズ数における計算結果は省略した。
図１１、図１２にそれぞれ実験データなどと比較した無次元時間５０〜１００にわたる抗力係数の平均値と、平均Ｓｔｒｏｕｈａｌ数を示す。抗力係数について、低レイノルズ数（Ｒｅ＜６０）の本発明による計算結果は今井の式による結果（Ｂ）とよく一致している。レイノルズ数Ｒｅ＝１００付近あたりにおいて、本発明による計算結果は実験データ（Ｃ）と比較してわずかに小さくなっている。レイノルズ数Ｒｅ＞４００では、本発明による計算結果は実験データ（Ｃ）とほとんど相違がない結果となっている。
これは主に本発明で用いる空間解像度（１０／６４＝０．１５）が低いので、Ｒｅ＝１００（Ｏ）オーダー以上になると、円柱壁における境界層は十分に解像されていないと考えられる。図１１から分かるようにレイノルズ数Ｒｅ＞５０では、カットセル無しの階段近似による計算結果（Ｄ）は実験より明らかに過大評価される。これに対して、本発明で提案したカットセルＶＯＦ法による計算は、計算結果（Ｄ）に比べて大幅に改善されていることが分かる。
図１２の平均Ｓｔｒｏｕｈａｌ数について、本発明での計算結果は、実験データより大きめであるが、レイノルズ数Ｒｅ＞２００の適正値とされる０．２〜０．２５の範囲に収まっていることがわかる。
図１３Ａ，Ｂに他の研究者による直交適合格子法［非特許文献１５］を使用した円柱の抗力係数とＳｔｒｏｕｈａｌ数の解析結果を示す。また図１４Ａ，Ｂに一般座標系適合格子法［非特許文献１６］による抗力係数とＳｔｒｏｕｈａｌ数の解析結果を示す。これら従来の方法による結果では、レイノルズ数Ｒｅ＞２００の適正値とされる０．２〜０．２５の範囲から部分的に外れていることがわかる。従って本発明による平均Ｓｔｒｏｕｈａｌ数の計算結果（図１１、図１２）は、これら従来の境界適合格子による計算結果と比べて計算精度が優れていることが分かる。
レイノルズ数Ｒｅ＝３００の場合に無次元時間Ｔ＝２．５の時の可視化実験［非特許文献１７］の結果と、本発明での同じ条件による計算結果の円柱近傍拡大図を図１５Ａ，Ｂに示す。またレイノルズ数Ｒｅ＝５５０の場合（無次元時間Ｔ＝２．５）に実験［非特許文献１７］と対照する円柱付近の計算結果を図１６Ａ，Ｂに示す。
上記の図から分かるようにいずれのケースで円柱付近の双子渦の形状や大きさ、円柱直後から渦の循環再合流までの距離など、本発明での計算結果は可視化実験とよく一致している。
図１７Ａ，Ｂにレイノルズ数Ｒｅ＝２００場合の、本発明の方法による円柱抗力と揚力の時間履歴変化を、朴ら［非特許文献６］により提案された部分境界適合直交格子法による計算結果と比較して示す。この図から、本発明の等間隔直交格子による抗力およびＳｔｒｏｕｈａｌ数は、部分境界適合直交格子法による結果とよく一致していることが確認された。なお、本発明の揚力の変動幅は部分境界適合直交格子法による結果より小さく評価されているが、これは主に［非特許文献６］の許算領域は３０Ｄ×１５Ｄであり、流出境界として放射条件を使用しているためと考えられる。
ここで図に示されていないが、本発明でレイノルズ数Ｒｅ＝１からＲｅ＝２００００までの２４ケース計算を行ったが、レイノルズ数８０以上でカルマン渦が出るようになった。実験や理論の臨界レイノルズ数Ｒｅ＝５０位に比べてやや大きめであるが、実用上十分であると思われる。
３．３ 実用問題に向けての応用例
今回開発した解析方法をいくつかの例に応用してみた。それらは図１８Ａ〜Ｆのようにバックステップを過ぎる流れ（Ａ）、狭窄管内の流れ（Ｂ）、チャネル内障害物周りの流れ（Ｃ）、分岐管内の流れ（Ｄ）、複数Ｂｌｕｆｆ Ｂｏｄｙまわりの流れ（Ｅ）や予混合燃焼器内の流れ（Ｆ）などである。
ここに示しているのはＸ方向速度ｕ分布のみであり、これらの応用例の計算は３．１のダクト計算や３．２の円柱まわりの計算と同じように計算領域は１０Ｄ×１０Ｄで、格子数は６４×６４とした。またレイノルズ数はＲｅ＝１００、流入条件や流出条件等は３．２の計算と同じにした。つまり、開発した計算プログラムを変更することなく、Ｖ−ＣＡＤデータを直接に入力するだけで、任意形状の流れ場を計算できるようになった。
４． 次に本発明の解析方法を移動境界に伴う流れの数値計算（流体構造連成解析）に適用する。流体構造連成解析は、風工学や熱交換器など工学的見地において特に重要である。
流体構造連成解析の具体例として、（１）一様流の中に置かれた二次元円柱の強制振動と、（２）円柱の渦自励振動（渦励振）の数値計算を行い、構造物に致命的な損傷を与える代表的な流体構造連成問題としてのロックイン現象や渦励振による共振状態を検証し、移動境界の流体構造連成解析における本発明の解析方法の有効性を示す。
４．１ 基礎方程式
円柱の振動計算に関しては、弾性支持された剛体構造物として、流れ方向に対して直角方向に振動する１質点１自由度系モデルを仮定している。
円柱が強制振動の場合には、円柱を一様流と直角な方向（Ｙ方向）にｓｉｎ波で強制的に振動させる。即ち円柱中心のＹ方向の変位を［数９］の式（１７）で与える。ここで、Ａは振幅、ｆｃは強制振動周波数、ｔは時間である。
また、円柱表面のＹ方向の速度境界条件は、［数９］の式（１８）で与えられる。
円柱が渦放出による自励振動をする場合には、有次元の振動方程式は１質点１自由度のダンパ・ばねモデルによって［数９］の式（１９）に示す。
式（１９）の中でそれぞれｍは円柱の質量、ｃはダンパの粘性減衰係数、ｋはバネの弾性回復係数、ρは流体の密度、Ｕ０は一様流の流速、Ｄは円柱の直径、ＣＬは揚力係数である。式（１９）を無次元化して円柱の固有振動数ｆ０＝（ｋ／ｍ）０．５／（２π）、無次元減衰係数ｈ＝ｃ／（２（ｍｋ）０．５）、スクルートン数Ｓｃ＝２πｈ（２ρｅ／ρ）などの振動応答を支配するパラメータで整理すると、［数９］の式（２０）になる。
４．２ 数値解法
円柱が振動する場合には、式（１８）または式（２０）から求められた振動する円柱の運動速度を、流れ場の速度境界条件として計算時間ステップごとに変化して与えている。
円柱の渦自励振動の解析については、円柱の初期変位と初期速度を０として、式（２０）の振動方程式をニューマークのβ法により積分し、円柱の振動変位や振動速度を求められる。
本発明の流体構造連成解析は、流体系と構造系を連立して統一的に扱う強連成ではなく、流体と構造の方程式を別々に離散化し、それぞれの境界条件を陽的に受け渡す弱連成となる。複雑な流体構造系の連成解析における時間刻みは、特に強連成の場合、流体または構造計算の時間刻みの小さい方に制約されることは原則であるが、本発明での構造系の計算には単純なダンパ・バネモデルを用いるため、計算の時間刻みは流体解析のＣＦＬ条件のみから決められる。また、ニューマークのβ法により将来の加速度を積分する際に、式（２０）の揚力は現在の値を用いて陽的に与える。
４．３ 強制振動円柱まわりの流れ
流れ場のレイノルズ数Ｒｅ＝１０００の場合には、３．２の解析条件のままで、式（１７）と式（１８）による強制振動円柱まわりの流れ解析を行った。レイノルズ数Ｒｅ＝１０００の際に、静止円柱の無次元渦放出周波数のＳｔｒｏｕｈａｌ数はＳｔ＝０．２位であることが実験値で知られている。強制振動の計算条件として振幅Ａを０．１、強制振動周波数を０．２に中心としてｆｃ＝０．０５、０．１、０．１５、０．２、０．２５、０．３、０．４、０．５の８ケースで計算を行なった。静止円柱を含めて各振動数での円柱に作用する流体力の抗力係数ＣＤおよび揚力係数ＣＬの時間履歴の変動波形をそれぞれ図１９Ａ〜Ｉに示す。
図１９Ａ〜Ｉの抗力係数と揚力係数の時間履歴波形図から、振動周波数ｆｃ＝０．５を除いて、振動する円柱に作用する抗力および揚力が静止円柱のものより大きいことが分かった。
揚力の変動波により求められた流れの渦放出状態を表わすＳｔｒｏｕｈａｌ数Ｓｔと強制振動周波数ｆｃの関係を図２０に示す。Ｓｔｒｏｕｈａｌ数０．２付近の０．１５＜ｆｃ＜０．２５におけるｆｃ／Ｓｔ＝１の計算結果から、強制振動による渦放出周波数のロックイン現象がよく再現されていることが分かった。また、図１９Ａ〜Ｉからロックイン領域での抗力や揚力が増大することも捉えている。
本発明とＯｋａｍｏｔｏらの実験［非特許文献１８］とはレイノルズ数（Ｒｅ＝１９０００）や円柱振動振幅（Ａ＝０．０５）とも異なるが、Ｏｋａｍｏｔｏらの実験による０．１８６＜ｆｃ＜０．２１６のロックイン範囲は本発明のロックイン領域に含まれることが確認された。
４．４ 円柱の渦放出による自励振動
４．３の解析条件で、ニューマークのβ法により式（２０）の振動方程式を積分し、円柱の渦自励振動の振動変位や振動速度を求める。円柱の振動速度をさらに流れ場の速度境界条件として課する。流れ場のレイノルズ数Ｒｅを１００００、スクルートン数Ｓｃを５、減衰係数ｈを１％、無次元流速Ｖｒ＝Ｕｏ／（ｆｏＤ）を２，３，４，５，６，７，８の七ケースで計算を行い、各々の無次元流速に対する抗力係数ＣＤ、揚力係数ＣＬ、そして変位ｙの時間履歴を図２１Ａ〜Ｃと図２２Ｄ〜Ｇに示す。
図２１Ａ〜Ｃと図２２Ｄ〜Ｇに示された円柱の渦励振解析について、流れ方向に対して直角方向に振動する弾性支持された剛体構造物の１質点１自由度系モデルを仮定しているが、図２１Ａ〜Ｃと図２２Ｄ〜Ｇの計算結果から分かるように、式（２０）の振動方程式の右辺の揚力ＣＬ項は、非定常な渦放出により不規則的な変化を呈している。そして流体構造の連成運動としてさらに振動円柱が流れ場に非定常な擾乱をも与えているため、４．３のｓｉｎ強制振動に比べて、渦励振による揚力と抗力の平均値は、静止円柱のものより増大することが一致しているものの、その時刻歴変化曲線はもっと不規則的な波形となっていることが見られる。また無次元流速Ｖｒ＝４（円柱の固有振動数ｆｏ＝１／４＝０．２５）の場合に円柱の振幅がもっとも大きくなることが分かった。これは静止円柱の流れ場のＳｔｒｏｕｈａｌ数Ｓｔ＝０．２５にも関連すると考えられる。
各々の無次元流速Ｖｒに対するｙの最大変位を整理して図２３に示す。最大変位ｙについてみると、Ｖｒ＝２、３、４ではｙの最大変位は次第に大きくなり、Ｖｒ＝４で最大になる。その後、Ｖｒが大きくなるにつれて、ｙの最大３変位は波状に小さくなっていく。本計算での渦励振応答振幅の変化傾向は、近藤［非特許文献１９］のＡＬＥ有限要素法による計算結果と定性的に一致していると確認された。
５．結論
（１）Ｖ−ＣＡＤデータを直接利用する流体解析方法を開発し、内部流・外部流ともにベンチマークテストを用いて検証した。その結果、本発明で提案した流体解析方法は十分な安定性を有し、二次元静止円柱まわりの流れの計算では、従来のボクセル法より良く、境界適合格子と比べて遜色がない計算精度を持つことが検証された。
（２）流れ場の形状に関しては、どんな複雑な形状があっても、Ｖ−ＣＡＤデータさえを読み込めば、短時間で直接二次元の非圧縮性粘性流れの数値シミュレーションができるようになった。
（３）二次元円柱の流体構造連成振動の計算では、円柱まわり流れのレイノルズ数Ｒｅ＝１０００の場合に、円柱の強制振動によるロックイン（ｆｃ／Ｓｔ＝１）領域は０．１５＜ｆｃ＜０．２５であることが予測された。
また、レイノルズ数Ｒｅ＝１００００ かつ無次元流速Ｖｒ＝４の時、つまり、円柱固有振動数（ｆｏ＝１／４）が静止円柱のＳｔｒｏｕｈａｌ数（Ｓｔ＝０．２５）と等しくなる時に、渦自励振動による共振状態が起こることが示された。
これらの計算結果は、実験データやＡＬＥ有限要素法による結果と定性的に一致している。したがって、移動境界の流体構造連成問題においても、提案した流体解析方法の有効性が示されたと言える。
図２は、本発明の数値解析方法を実行するための数値解析装置の構成図である。
図３は、本発明の数値解析方法とそのプログラムのフロー図である。
図４Ａと図４Ｂは、二次元における開口率の定義を示す図である。
図５は、境界セルにおける有限体積法による空間離散化の説明図である。
図６Ａ，Ｂ，Ｃは、解析した３ケースのＶＯＦ分布を示す画像である。
図７Ａ，Ｂ，Ｃは、解析した３ケースの流れ方向速度の分布と理論解析解との比較図である。
図８は、解析格子とＶＯＦ分布を示す画像である。
図９は、速度ベクトルを示す画像である。
図１０は、圧力等値線を示す画像である。
図１１は、本発明による平均抗力係数を示す図である。
図１２は、本発明による平均Ｓｔｒｏｕｈａｌ数を示す図である。
図１３Ａは、直交適合格子法による静止円柱の抗力係数を示す図であり、図１３Ｂは、その平均Ｓｔｒｏｕｈａｌ数を示す図である。
図１４Ａは、一般座標系境界適合格子による静止円柱の抗力係数を示す図であり、図１４Ｂはその平均Ｓｔｒｏｕｈａｌ数を示す図である。
図１５Ａは、Ｒｅ＝３００，Ｔ＝２．５の場合の実験による流脈線を示す画像であり、図１５Ｂはその流線を示す画像である。
図１６Ａは、Ｒｅ＝５５０，Ｔ＝２．５の場合の実験による流脈線を示す画像であり、図１６Ｂは、その流線を示す画像である。
図１７Ａは、抗力と揚力の履歴変化を示す本発明の結果であり、図１７Ｂは、その部分境界適合直交格子法の結果である。
図１８Ａ〜Ｆは、実用問題に向けての応用例を示す画像である。
図１９Ａ〜Ｉは、静止円柱と各振動数での円柱に作用する流体力の抗力係数ＣＤおよび揚力係数ＣＬの時間履歴の変動波形を示す図である。
図２０は、流れの渦放出状態を表すＳｔｒｏｕｈａｌ数Ｓｔと強制振動周波数ｆｃの関係図である。
図２１Ａ〜Ｃは、無次元流速Ｖｒが２，３，４における変位ｙ（左図）と抗力係数ＣＤ及び揚力係数ＣＬ（右図）の時間履歴を示す図である。
図２２Ｄ〜Ｇは、無次元流速Ｖｒが５，６，７，８における変位ｙ（左図）と抗力係数ＣＤ及び揚力係数ＣＬ（右図）の時間履歴を示す図である。
図２３は、各々の無次元流速Ｖｒに対するｙの最大変位を示す図である。
好ましい実施例の説明
以下、本発明の好ましい実施形態を図面を参照して説明する。
図２は、本発明の数値解析方法を実行するための数値解析装置の構成図である。
この図に示すように、本発明の数値解析装置１０は、入力装置２、外部記憶装置３、内部記憶装置４、中央処理装置５および出力装置６を備える。
入力装置２は、例えばキーボードであり、対象物１の形状データからなる外部データ１２を入力する。外部記憶装置３は、ハードディスク、フロッピィーディスク、磁気テープ、コンパクトディスク等であり、形状と物理量を統合した実体データとその記憶演算プログラムを記憶する。内部記憶装置４は、例えばＲＡＭ，ＲＯＭ等であり、演算情報を保管する。中央処理装置５（ＣＰＵ）は、演算や入出力等を集中的に処理し、内部記憶装置４と共に、記憶プログラムを実行する。出力装置６は、例えば表示装置とプリンタであり、記憶した実体データと記憶プログラムの実行結果を出力するようになっている。
本発明の記憶演算装置１０は、上述した外部記憶装置３、内部記憶装置４、及び中央処理装置５により、外部データを境界が直交する複数のセル１３に分割し、分割された各セルを対象物の内側に位置する内部セル１３ａと境界データを含む境界セル１３ｂとに区分し、境界データによる境界セル１３ｂの稜線の切断点を求め、求めた切断点を結ぶ多角形を境界面のセル内部データとし、流れ場の境界に、ＶＯＦ法を併用したカットセル有限体積法を適用して解析する。
図３は、本発明の数値解析方法とそのプログラムのフロー図である。この図に示すように、本発明の方法及び変換プログラムは、分割ステップ（Ａ）、セル区分ステップ（Ｂ）、切断点決定ステップ（Ｃ）、境界面決定ステップ（Ｄ）、及び解析ステップ（Ｅ）からなる。
外部から入力する外部データ１２は、多面体を表すポリゴンデータ、有限要素法に用いる四面体又は六面体要素、３次元ＣＡＤ又はＣＧツールに用いる曲面データ、或いはその他の立体の表面を部分的な平面や曲面で構成された情報で表現するデータである。
外部データ１２は、このようなデータ（Ｓ−ＣＡＤデータと呼ぶ）のほかに、（１）Ｖ−ＣＡＤ独自のインターフェース（Ｖ−ｉｎｔｅｒｆａｃｅ）により人間の入力により直接作成されたデータと、（２）測定機やセンサ、デジタイザなどの表面のデジタイズデータや、（３）ＣＴスキャンやＭＲＩ、および一般的にＶｏｌｕｍｅレンダリングに用いられているボクセルデータなどの内部情報ももつＶｏｌｕｍｅデータであってもよい。
分割ステップ（Ａ）では、外部データ取得ステップ（図示せず）で取得した非圧縮性粘性流体と接する対象物の境界データからなる外部データ１２を境界平面が直交する複数のセル１３に分割する。この分割は、３次元の場合は、八分木分割であり、２次元の場合には４分割である。
すなわちこの分割ステップ（Ａ）における分割とは、目的とする非圧縮性粘性流体と接する対象物を含む、基準となる直方体（または矩形）を分割（８分割又は４分割）し、それぞれの領域の中に立体が完全に含まれるか、含まれなくなるまで再帰的に分割処理を繰り返す。この分割によりボクセル表現よりも大幅にデータ量を減らすことができる。
空間分割により分割された一つの空間領域をセル１３とよぶ。セルは境界が直交する直方体または矩形である。セルによる階層構造、分割数もしくは分解能によって空間中に占める領域を表現する。これにより空間全体の中で対象は大きさの異なるセルを積み重ねたものとして表現される。
セル区分ステップ（Ｂ）では、分割された各セルを対象物の内側に位置する内部セル１３ａと境界データを含む境界セル１３ｂとに区分する。
すなわち本発明では非圧縮性粘性流体と接する対象物の内側または外側に完全に内部に含まれるものはその最大の大きさをもつ内部セル１３ａ（立方体）とし、外部データ１２からの境界情報を含むセルは境界セル１３ｂとする。
切断点決定ステップ（Ｃ）では、境界データによる境界セル１３ｂの稜線の切断点を求める。
境界面決定ステップ（Ｄ）では、求めた切断点を結ぶ多角形を境界面のセル内部データとする。以下、かかる切断点を結ぶ多角形を含むセルを「カットセル」と呼ぶ。
解析ステップ（Ｅ）では、上述した内部セル１３ａと境界セル１３ｂに対して、流れ場の境界に、ＶＯＦ法を併用したカットセル有限体積法を適用して非圧縮性粘性流体の流れ場の数値解析を行う。
以下、本発明を更に詳細に説明する。
１．本発明ではＶ−ＣＡＤプログラムにおける実用性のある流体解析技術を目指し、カットセル（ＫＴＣ）直交格子による任意形状の非圧縮粘性流れの解析方法を開発した。本発明において、流れ場の境界での扱いについては、Ｈｉｒｔらが提案しているＶＯＦ（Ｖｏｌｕｍｅ Ｏｆ Ｆｌｕｉｄ）法［非特許文献１２］を併用したカットセル有限体積法を用いた。また本発明の方法による内部流であるチャンネル内流れと、外部流の円柱まわり流れを解析し、実験データ、理論解、及び既存の方法による解析結果との比較を行った。更に、本解析法の応用例として、いくつかの計算結果を合わせて示した。
２．基礎方程式と計算方法
２．１ 基礎方程式
本発明で用いる基礎支配方程式は［数４］の式（１）（２）に示す非圧縮粘性流体のＮａｖｉｅｒ・Ｓｔｏｋｅｓ方程式と連続の式である。
２．２ 計算精度・計算コスト及び計算方法
流体のシミュレーションは一種の数値実験であり、つねに一定の誤差を伴う。流体解析を精度よく行うには、以下の四つの条件が必要となる。
▲１▼流れの最小長さスケール（境界層・渦・衝撃波・火炎面など）を捕られるほど細かい空間解像度。
▲２▼流れの最大長さスケールを十分に捕獲するほど大きい計算領域と、人為的な流入・流出境界条件・壁によるブロッキングなど影響を無視できるほどの計算領域。
▲３▼打ち切り誤差や数値拡散などを無視できるほど十分な空間的・時間的な離散化精度。
▲４▼特定問題に即したモデル（壁面モデル・乱流ＳＧＳモデル・燃焼モデルなど）やスキーム（Ｋ−Ｋ・ＱＵＩＣＫ・ＭＵＳＣＬ・ＴＶＤ・ＥＮＯなど）、差分法を援用すれば、できるだけ直交等方性のある格子。
流体解析の計算コスト（計算時間・必要なメモリ）は、その計算方法に求められる解析精度によって決められるといえる。実用的な流体解析では、ある程度の精度（例えば、普通の物理実験の精度範囲や、ユーザーに求められる精度）があれば、十分と考えられる。
一方、流れ場は内部流、外部流及びそのほか（例えば噴流）に分類される。一般に外部流の解析領域は内部流のものより大きくとるのが普通である。本発明では実用性の観点から、外部流の場合には、流れ場の代表長さスケールＤ（＝１）として、解析領域の大きさはすべて１０Ｄ×１０Ｄとした。
本発明では、差分法を併用した有限体積法を用いる。ある程度の安定性や離散化精度を確保するため、空間積分について、対流項に二次元のＱＵＩＣＫ補間スキームを使用し、拡散項に二次精度の中心差分を用いる。時間進行について、対流項と拡散項を合わせて二次精度のＡｄａｍｓ−Ｂａｓｈｆｏｒｔｈ法を適用し、圧力勾配項に一次精度のＥｕｌｅｒ陰解法を用いる。解析アルゴリズムとしての圧力と速度のカップリングには、Ｈｉｒｔらにより提案された圧力と速度を同時緩和するＳＯＬＡ−ＨＳＭＡＣ法［非特許文献１３］を使用し、緩和係数１．６５でＳＯＲ法を用いて反復計算を行う。収束判定は連続方程式（２）の残差が０．０００２とする。また圧力振動を防ぐために三方向の速度ｕ，ｖ，ｗ，及び圧力ｐの格子上の定義点を、半メッシュをずれるスタガード格子を配置する。
２．３ Ｖ−ＣＡＤデータによる固体・流体境界の取扱い
Ｖｏｌｕｍｅ−ＣＡＤシステムの中でもっとも重要なアルゴリズムは加瀬・手嶋［非特許文献１４］により提案されるセルの切断点からモノの表面形状を復元するＫＴＣアルゴリズムである。このＫＴＣアルゴリズムは流体解析分野における直交格子カットセル（Ｃｕｔ Ｃｅｌｌ）法と同じ概念であることを注意されたい。また、ＫＴＣは二次元においても一個のセルの四つの稜線上に２個以上の交点（切断点）があり得るが、本発明では簡単のために二次元で一個のセルの稜線上に２個だけの切断点に限るとする。この場合のＫＴＣは図４Ａと図４Ｂに示すように二種類しかないことがわかる。
カットセル法における開口率は、セルの稜線に流体部分が占める比率として図４Ａと図４Ｂに示すように流体・固体境界とセル稜線との交点情報から［数５］の式（３）で定義される。
ここで、Δｘｉ，ｊ、Δｙｉ，ｊはそれぞれＸ方向とＹ方向の格子幅であり、Ａｉ， ｊ、Ｂｉ，ｊはＸ方向とＹ方向の開口率である。また、Ｖｘｉ，ｊ、Ｖｙｉ，ｊはそのセル稜線上に固体が 占める線分を表す。
各方向の開口率が分かれば、ＶＯＦ法における流体の体積占有率は［数５］の式（４）と式（５）により求められる。
有限体積法で保存型の支配方程式を離散化する場合、方程式（１）の積分形をグリーン定理（発散定理）に適用するため、微分型の方程式（１）をテンソル形の発散型に書き換えると［数６］の式（６）のようになる。
ここで、○に×はテンソル積であり、ベクトルＩはＫｒｏｎｅｃｋｅｒδｉｊに対応する単位テンソルである。
流体固体境界での離散化は、本発明の肝心なところになる。境界セル以外の格子における有限体積法による直交格子での離散化は従来の方法でもあまり問題が生じないのでここでは省略する。流体固体境界では図４のようなＫＴＣセルを想定して式（６）の空間離散化を行う。また式（６）の時間離散化について、流体占有率ＶＯＦを用いてコントロールボリューム（ＣＶ）内における速度の時間勾配を平均値とみなして、従来の方法通りで行えるため、境界セルでの時間離散化の説明も省く。
さて、図５の二次元境界セルの流体部分をコントロールボリューム（ＣＶ）Ｖｉ，ｊとして、式（６）を空間積分する（ここで簡単のため外力項ベクトルｆ＝０とした）と、有限体積法における支配方程式は［数６］の式（７）になる。
カットセルによる離散化には、二つのキーポイントがあることを注意されたい。
キーポイントの一つは固体境界（図５の切断線４）の積分の取り扱いである。有限体積法では、一階偏微分の対流項と圧力勾配項の境界条件はディリクレ条件であり、二階偏微分の拡散項の境界条件はノイマン問題となる。固体境界でノースリップ境界条件の場合に、図５の切断線４のところに対流項が零となるから対流項を積分しなくてもいいが、圧力勾配項と拡散項に対して切断線４のところに積分しなければいけない。これは式（９）と式（１０）の各々最後の項に切断線４の中間点Ｐの値を（平均値として）用いて積分するわけである。また空間積分に対してすべての項に開口率を適用する必要もあることを注意されたい。
もう一つのキーポイントは、カットセルにおいて算出される物理変数の定義点である。有限体積法では流体占有率ＶＯＦで体積平均を行うため、物理変数の定義点は流体のコントロールボリュームの中心（例えば図５の○に・点で示すＣ点）に置くのが自然だが、そうすると直交等間隔格子の完全自動化メリットがなくなる。その故、本発明ではＶＯＦ＝０．０１を閾値としてそれより小さい境界セルは完全の固体とみなし、それより大きい境界セルにおいて算出される変数はコントロールボリュームの中心ではなく、セルの中心（例えば図５の○に＋印で示すＤ点）に置くこととする。同じく図５の稜線３と稜線５の変数定義点も線分の中心○ではなく、セル稜線の中心◎に定義される。また線分４の中心点Ｐに変数が定義されていないため、本発明でその点の変数値を線形補間により求める。こう処理すると、固体境界の再現性を損するものの、直交格子の完全自動化が保たれ、実用化に向ける。また流体の計算はコントロールボリュームで保存則が厳密に満たされ、上記の処理は流れ場の計算精度に及ぼす影響がそれほど大きくないと考えられる。
２．５ 物体に働く抗力と揚力について
物体まわりの流れを解析するとき、その物体に働く抗力と揚力の計算について、一般座標系の境界適合格子の場合に物体の第１格子点に沿って積分すれば簡単だが、直交格子で取り扱う場合にはやや煩雑となる。そこで本発明では、以下の方法で抗力と揚力を求める。
物体に働く力は、その物体表面に作用する流体応力テンソルを閉積分して、グリーン定理を適用すれば、二次元の場合に［数８］の式（１１）ように得られる。
ここでσｎは法線方向の平面に働く応力テンソル、σは物体表面の点での応力テンソルである。また二次元の二階応力テンソルの発散はそのテンソルの四つの方向からなることを注意されたい。それでＸ方向を流れ方向、Ｙ方向を流れの垂直方向とすると、物体に働く抗力（流れ方向の力）と揚力（流れの垂直方向の力）は、グリーン定理を再び適用すると［数８］の式（１２）（１３）のようになる。
抗力と揚力が分かれば、無次元量としての抗力係数ＣＤと揚力係数ＣＬを、主流速度Ｕ、流体密度ρ及び物体の基準面積Ａ（二次元の場合は基準長さ）を用いて［数５］の式（１４）のように無次元化される。
また非定常現象を表わす無次元量のストローハル数Ｓｔは、ＣＬの周波数ｆＬから［数８］の式（１５）で求める。
流体の数値解析において解析スキームや数値モデルなどを検証することにあたって、内部流のベンチマークテストとしてチャネル流れと、外部流のベンチマークテストとしての円柱まわりの流れはしばしば検証の対象となる。本発明でも、上記のＶＯＦ法を併用したカットセル有限体積法を検証するために、解析理論解のある二次元Ｐｏｉｓｅｕｉｌｌｅ流と、多くの実験や計算データのある静止円柱まわりの流れについて数値解析を行い、他の研究者のデータとの比較を行った。なお、実用問題に向けての応用例として、バックステップを過ぎる流れ、狭窄管内の流れ、チャネル内障害物周りの流れ、分岐管内の流れ、複数Ｂｌｕｆｆ Ｂｏｄｙまわりの流れや予混合燃焼器内の流れなどについて解析を行い、その計算結果の可視化例を示す。
３．１ チャネルＰｏｉｓｅｕｉｌｌｅ流による検証
セルのカットされたＶＯＦの効果を検証するため、１０Ｄ×１０Ｄ正方形の解析領域におけるダクトを、０度、１０度及び４５度を斜めにして解析を行った。全体の格子数は６４×６４、ダクトの半幅は２．５Ｄ、この半幅で定義されるレイノルズ数Ｒｅ＝１、時間刻み△ｔ＝０．０００１、解析時間ｓｔｅｐｓ＝１００００とする。ダクトの流入条件と流出条件はともに放物型の分布を仮定した。流れ方向の圧力勾配は１．０とした。またＨＳＭＡＣ法を用い、入口・出口および壁での圧力境界条件は直接与えていない。三ケースのＶＯＦ分布を図６Ａ，Ｂ，Ｃに示す。図６Ａ，Ｂ，Ｃにおける破線は、境界を示している。その三ケースに対応する流れ方向速度分布を理論解析解と比較して図７Ａ，Ｂ，Ｃに示す。
図７Ａ，Ｂ，Ｃの中に線Ａのａｎａｌｙｓｉｓは理論解析解であり、線Ｂのｉ＝３２は第３２格子点断面にある流れ軸中心での流れ垂直断面の速度分布である。ダクトを傾けない（角度０）計算結果の速度分布は解析解とよく一致し、傾斜角度を１０度、４５度にすると、軸中心の最大速度がわずかに過小評価されるものの、全体として解析結果は理論解と一致していると言える。
３．２ 静止円柱まわり流れによる検証
解析スキームやモデルなどの検証対象となる円柱まわり流れの解析は、直交格子にとって最も困難な例題であろう。本発明では、３．１の解析領域（１０Ｄ×１０Ｄ）と格子数（６４×６４）のままで、円柱直径Ｄで定義されるレイノルズ数のＲｅ＝１からＲｅ＝２００００までの２４ケースについて、一様流中に置かれている静止円柱まわり流れの解析を行った。計算の無次元時間は１００とし、円柱の円心座標の図６のようにＸ＝２．５Ｄ，Ｙ＝５．０Ｄとした。速度の境界条件として、円柱表面で滑り無し条件、流入側の流入条件は一様流（流れ方向速度＝１．０）条件、流出側の流出条件はＸ方向に自由流出（速度勾配０），Ｙ方向にフリースリップ条件とした。また、ＨＳＭＡＣ法を用いて速度・圧力を同時に緩和するため、流入・流出および壁での圧力境界条件を直接与える必要がない。
図８は解析格子と解析全領域のＶＯＦ分布を示す。図９、図１０はレイノルズ数Ｒｅ＝３００の場合の、速度ベクトルと圧力等値線であり、円柱後方に周期的なカルマン渦列の形成が観察されており、非定常流の流れパタンーの特徴を捉えている。これはレイノルズ数Ｒｅ＞８０の場合に共通なので、他のレイノルズ数における計算結果は省略した。
図１１、図１２にそれぞれ実験データなどと比較した無次元時間５０〜１００にわたる抗力係数の平均値と、平均Ｓｔｒｏｕｈａｌ数を示す。抗力係数について、低レイノルズ数（Ｒｅ＜６０）の本発明による計算結果は今井の式による結果（Ｂ）とよく一致している。レイノルズ数Ｒｅ＝１００付近あたりにおいて、本発明による計算結果は実験データ（Ｃ）と比較してわずかに小さくなっている。レイノルズ数Ｒｅ＞４００では、本発明による計算結果は実験データ（Ｃ）とほとんど相違がない結果となっている。
これは主に本発明で用いる空間解像度（１０／６４＝０．１５）が低いので、Ｒｅ＝１００（Ｏ）オーダー以上になると、円柱壁における境界層は十分に解像されていないと考えられる。図１１から分かるようにレイノルズ数Ｒｅ＞５０では、カットセル無しの階段近似による計算結果（Ｄ）は実験より明らかに過大評価される。これに対して、本発明で提案したカットセルＶＯＦ法による計算は、計算結果（Ｄ）に比べて大幅に改善されていることが分かる。
図１２の平均Ｓｔｒｏｕｈａｌ数について、本発明での計算結果は、実験データより大きめであるが、レイノルズ数Ｒｅ＞２００の適正値とされる０．２〜０．２５の範囲に収まっていることがわかる。
図１３Ａ，Ｂに他の研究者による直交適合格子法［非特許文献１５］を使用した円柱の抗力係数とＳｔｒｏｕｈａｌ数の解析結果を示す。また図１４Ａ，Ｂに一般座標系適合格子法［非特許文献１６］による抗力係数とＳｔｒｏｕｈａｌ数の解析結果を示す。これら従来の方法による結果では、レイノルズ数Ｒｅ＞２００の適正値とされる０．２〜０．２５の範囲から部分的に外れていることがわかる。従って本発明による平均Ｓｔｒｏｕｈａｌ数の計算結果（図１１、図１２）は、これら従来の境界適合格子による計算結果と比べて計算精度が優れていることが分かる。
レイノルズ数Ｒｅ＝３００の場合に無次元時間Ｔ＝２．５の時の可視化実験［非特許文献１７］の結果と、本発明での同じ条件による計算結果の円柱近傍拡大図を図１５Ａ，Ｂに示す。またレイノルズ数Ｒｅ＝５５０の場合（無次元時間Ｔ＝２．５）に実験［非特許文献１７］と対照する円柱付近の計算結果を図１６Ａ，Ｂに示す。
上記の図から分かるようにいずれのケースで円柱付近の双子渦の形状や大きさ、円柱直後から渦の循環再合流までの距離など、本発明での計算結果は可視化実験とよく一致している。
図１７Ａ，Ｂにレイノルズ数Ｒｅ＝２００場合の、本発明の方法による円柱抗力と揚力の時間履歴変化を、朴ら［非特許文献６］により提案された部分境界適合直交格子法による計算結果と比較して示す。この図から、本発明の等間隔直交格子による抗力およびＳｔｒｏｕｈａｌ数は、部分境界適合直交格子法による結果とよく一致していることが確認された。なお、本発明の揚力の変動幅は部分境界適合直交格子法による結果より小さく評価されているが、これは主に［非特許文献６］の許算領域は３０Ｄ×１５Ｄであり、流出境界として放射条件を使用しているためと考えられる。
ここで図に示されていないが、本発明でレイノルズ数Ｒｅ＝１からＲｅ＝２００００までの２４ケース計算を行ったが、レイノルズ数８０以上でカルマン渦が出るようになった。実験や理論の臨界レイノルズ数Ｒｅ＝５０位に比べてやや大きめであるが、実用上十分であると思われる。
３．３ 実用問題に向けての応用例
今回開発した解析方法をいくつかの例に応用してみた。それらは図１８Ａ〜Ｆのようにバックステップを過ぎる流れ（Ａ）、狭窄管内の流れ（Ｂ）、チャネル内障害物周りの流れ（Ｃ）、分岐管内の流れ（Ｄ）、複数Ｂｌｕｆｆ Ｂｏｄｙまわりの流れ（Ｅ）や予混合燃焼器内の流れ（Ｆ）などである。
ここに示しているのはＸ方向速度ｕ分布のみであり、これらの応用例の計算は３．１のダクト計算や３．２の円柱まわりの計算と同じように計算領域は１０Ｄ×１０Ｄで、格子数は６４×６４とした。またレイノルズ数はＲｅ＝１００、流入条件や流出条件等は３．２の計算と同じにした。つまり、開発した計算プログラムを変更することなく、Ｖ−ＣＡＤデータを直接に入力するだけで、任意形状の流れ場を計算できるようになった。
４． 次に本発明の解析方法を移動境界に伴う流れの数値計算（流体構造連成解析）に適用する。流体構造連成解析は、風工学や熱交換器など工学的見地において特に重要である。
流体構造連成解析の具体例として、（１）一様流の中に置かれた二次元円柱の強制振動と、（２）円柱の渦自励振動（渦励振）の数値計算を行い、構造物に致命的な損傷を与える代表的な流体構造連成問題としてのロックイン現象や渦励振による共振状態を検証し、移動境界の流体構造連成解析における本発明の解析方法の有効性を示す。
４．１ 基礎方程式
円柱の振動計算に関しては、弾性支持された剛体構造物として、流れ方向に対して直角方向に振動する１質点１自由度系モデルを仮定している。
円柱が強制振動の場合には、円柱を一様流と直角な方向（Ｙ方向）にｓｉｎ波で強制的に振動させる。即ち円柱中心のＹ方向の変位を［数９］の式（１７）で与える。ここで、Ａは振幅、ｆｃは強制振動周波数、ｔは時間である。
また、円柱表面のＹ方向の速度境界条件は、［数９］の式（１８）で与えられる。
円柱が渦放出による自励振動をする場合には、有次元の振動方程式は１質点１自由度のダンパ・ばねモデルによって［数９］の式（１９）に示す。
式（１９）の中でそれぞれｍは円柱の質量、ｃはダンパの粘性減衰係数、ｋはバネの弾性回復係数、ρは流体の密度、Ｕ０は一様流の流速、Ｄは円柱の直径、ＣＬは揚力係数である。式（１９）を無次元化して円柱の固有振動数ｆ０＝（ｋ／ｍ）０．５／（２π）、無次元減衰係数ｈ＝ｃ／（２（ｍｋ）０．５）、スクルートン数Ｓｃ＝２πｈ（２ρｅ／ρ）などの振動応答を支配するパラメータで整理すると、［数９］の式（２０）になる。
円柱が振動する場合には、式（１８）または式（２０）から求められた振動する円柱の運動速度を、流れ場の速度境界条件として計算時間ステップごとに変化して与えている。
円柱の渦自励振動の解析については、円柱の初期変位と初期速度を０として、式（２０）の振動方程式をニューマークのβ法により積分し、円柱の振動変位や振動速度を求められる。
本発明の流体構造連成解析は、流体系と構造系を連立して統一的に扱う強連成ではなく、流体と構造の方程式を別々に離散化し、それぞれの境界条件を陽的に受け渡す弱連成となる。複雑な流体構造系の連成解析における時間刻みは、特に強連成の場合、流体または構造計算の時間刻みの小さい方に制約されることは原則であるが、本発明での構造系の計算には単純なダンパ・バネモデルを用いるため、計算の時間刻みは流体解析のＣＦＬ条件のみから決められる。また、ニューマークのβ法により将来の加速度を積分する際に、式（２０）の揚力は現在の値を用いて陽的に与える。
４．３ 強制振動円柱まわりの流れ
流れ場のレイノルズ数Ｒｅ＝１０００の場合には、３．２の解析条件のままで、式（１７）と式（１８）による強制振動円柱まわりの流れ解析を行った。レイノルズ数Ｒｅ＝１０００の際に、静止円柱の無次元渦放出周波数のＳｔｒｏｕｈａｌ数はＳｔ＝０．２位であることが実験値で知られている。強制振動の計算条件として振幅Ａを０．１、強制振動周波数を０．２に中心としてｆｃ＝０．０５、０．１、０．１５、０．２、０．２５、０．３、０．４、０．５の８ケースで計算を行なった。静止円柱を含めて各振動数での円柱に作用する流体力の抗力係数ＣＤおよび揚力係数ＣＬの時間履歴の変動波形をそれぞれ図１９Ａ〜Ｉに示す。
図１９Ａ〜Ｉの抗力係数と揚力係数の時間履歴波形図から、振動周波数ｆｃ＝０．５を除いて、振動する円柱に作用する抗力および揚力が静止円柱のものより大きいことが分かった。
揚力の変動波により求められた流れの渦放出状態を表わすＳｔｒｏｕｈａｌ数Ｓｔと強制振動周波数ｆｃの関係を図２０に示す。Ｓｔｒｏｕｈａｌ数０．２付近の０．１５＜ｆｃ＜０．２５におけるｆｃ／Ｓｔ＝１の計算結果から、強制振動による渦放出周波数のロックイン現象がよく再現されていることが分かった。また、図１９Ａ〜Ｉからロックイン領域での抗力や揚力が増大することも捉えている。
本発明とＯｋａｍｏｔｏらの実験［非特許文献１８］とはレイノルズ数（Ｒｅ＝１９０００）や円柱振動振幅（Ａ＝０．０５）とも異なるが、Ｏｋａｍｏｔｏらの実験による０．１８６＜ｆｃ＜０．２１６のロックイン範囲は本発明のロックイン領域に含まれることが確認された。
４．４ 円柱の渦放出による自励振動
４．３の解析条件で、ニューマークのβ法により式（２０）の振動方程式を積分し、円柱の渦自励振動の振動変位や振動速度を求める。円柱の振動速度をさらに流れ場の速度境界条件として課する。流れ場のレイノルズ数Ｒｅを１００００、スクルートン数Ｓｃを５、減衰係数ｈを１％、無次元流速Ｖｒ＝Ｕｏ／（ｆｏＤ）を２，３，４，５，６，７，８の七ケースで計算を行い、各々の無次元流速に対する抗力係数ＣＤ、揚力係数ＣＬ、そして変位ｙの時間履歴を図２１Ａ〜Ｃと図２２Ｄ〜Ｇに示す。
図２１Ａ〜Ｃと図２２Ｄ〜Ｇに示された円柱の渦励振解析について、流れ方向に対して直角方向に振動する弾性支持された剛体構造物の１質点１自由度系モデルを仮定しているが、図２１Ａ〜Ｃと図２２Ｄ〜Ｇの計算結果から分かるように、式（２０）の振動方程式の右辺の揚力ＣＬ項は、非定常な渦放出により不規則的な変化を呈している。そして流体構造の連成運動としてさらに振動円柱が流れ場に非定常な擾乱をも与えているため、４．３のｓｉｎ強制振動に比べて、渦励振による揚力と抗力の平均値は、静止円柱のものより増大することが一致しているものの、その時刻歴変化曲線はもっと不規則的な波形となっていることが見られる。また無次元流速Ｖｒ＝４（円柱の固有振動数ｆｏ＝１／４＝０．２５）の場合に円柱の振幅がもっとも大きくなることが分かった。これは静止円柱の流れ場のＳｔｒｏｕｈａｌ数Ｓｔ＝０．２５にも関連すると考えられる。
各々の無次元流速Ｖｒに対するｙの最大変位を整理して図２３に示す。最大変位ｙについてみると、Ｖｒ＝２、３、４ではｙの最大変位は次第に大きくなり、Ｖｒ＝４で最大になる。その後、Ｖｒが大きくなるにつれて、ｙの最大３変位は波状に小さくなっていく。本計算での渦励振応答振幅の変化傾向は、近藤［非特許文献１９］のＡＬＥ有限要素法による計算結果と定性的に一致していると確認された。
５．結論
（１）Ｖ−ＣＡＤデータを直接利用する流体解析方法を開発し、内部流・外部流ともにベンチマークテストを用いて検証した。その結果、本発明で提案した流体解析方法は十分な安定性を有し、二次元静止円柱まわりの流れの計算では、従来のボクセル法より良く、境界適合格子と比べて遜色がない計算精度を持つことが検証された。
（２）流れ場の形状に関しては、どんな複雑な形状があっても、Ｖ−ＣＡＤデータさえを読み込めば、短時間で直接二次元の非圧縮性粘性流れの数値シミュレーションができるようになった。
（３）二次元円柱の流体構造連成振動の計算では、円柱まわり流れのレイノルズ数Ｒｅ＝１０００の場合に、円柱の強制振動によるロックイン（ｆｃ／Ｓｔ＝１）領域は０．１５＜ｆｃ＜０．２５であることが予測された。
また、レイノルズ数Ｒｅ＝１００００ かつ無次元流速Ｖｒ＝４の時、つまり、円柱固有振動数（ｆｏ＝１／４）が静止円柱のＳｔｒｏｕｈａｌ数（Ｓｔ＝０．２５）と等しくなる時に、渦自励振動による共振状態が起こることが示された。
これらの計算結果は、実験データやＡＬＥ有限要素法による結果と定性的に一致している。したがって、移動境界の流体構造連成問題においても、提案した流体解析方法の有効性が示されたと言える。
上述したように、次世代ＣＡＤとして期待されているＶ−ＣＡＤボクセルデータを直接利用する任意形状の非圧縮粘性流れ場の数値解析方法の開発を行った。流れ場の境界には、ＶＯＦ法を併用したカットセル有限体積法を用いた。二次元のチャンネル内流れと円柱まわり流れを用いて調べた結果、この方法はある程度の精度と安定性を有し、計算コストがあまりかからない実用計算方法であることが分かった。
従って、本発明のＶ−ＣＡＤデータを直接用いた非圧縮性粘性流体の流れ場の数値解析方法と装置は、格子生成を完全に自動化することができ、物体境界の表現が容易であり、比較的簡単な計算処理により、精度の高いシミュレーションを短時間で行うことができる等の優れた効果を有する。
従って、本発明のＶ−ＣＡＤデータを直接用いた非圧縮性粘性流体の流れ場の数値解析方法と装置は、格子生成を完全に自動化することができ、物体境界の表現が容易であり、比較的簡単な計算処理により、精度の高いシミュレーションを短時間で行うことができる等の優れた効果を有する。
Claims (15)
- 非圧縮性粘性流体と接する対象物の境界データからなる外部データ（１２）を境界が直交する複数のセル（１３）に分割する分割ステップ（Ａ）と、
分割された各セルを対象物の内側に位置する内部セル（１３ａ）と境界データを含む境界セル（１３ｂ）とに区分するセル区分ステップ（Ｂ）と、
前記境界データによる境界セル（１３ｂ）の稜線の切断点を求める切断点決定ステップ（Ｃ）と、
求めた切断点を結ぶ多角形を境界面のセル内部データとする境界面決定ステップ（Ｄ）と、
流れ場の境界に、ＶＯＦ法を併用したカットセル有限体積法を適用して解析する解析ステップ（Ｅ）と、を備えたことを特徴とするＶ−ＣＡＤデータを直接用いた非圧縮性粘性流体の流れ場の数値解析方法。 - 前記解析ステップ（Ｅ）において、空間積分について、対流項に二次元のＱＵＩＣＫ補間スキームを使用し、拡散項に二次精度の中心差分を用い、時間進行について、対流項と拡散項を合わせて二次精度のＡｄａｍｓ−Ｂａｓｈｆｏｒｔｈ法を適用し、圧力勾配項に一次精度のＥｕｌｅｒ陰解法を用いる、ことを特徴とする請求項１に記載の非圧縮性粘性流体の流れ場の数値解析方法。
- 二次元境界セルにおいて、有限体積法における支配方程式を［数１］の式（７）であらわす、
- 有限体積法における支配方程式における対流項、圧力勾配項および拡散項を［数２］の式（８）、（９）、（１０）であらわす、
- 固体境界の積分において、固体境界でノースリップ境界条件の場合に、対流項は零とし、圧力勾配項と拡散項に対しては、切断線の中間点Ｐの値を平均値として用いて積分し、空間積分に対してはすべての項に開口率を適用する、ことを特徴とする請求項３に記載の非圧縮性粘性流体の流れ場の数値解析方法。
- カットセルにおいて算出される物理変数の定義点を、ＶＯＦ＝０．０１を閾値としてそれより小さい境界セルは完全な固体とみなし、それより大きい境界セルにおいて算出される変数をセルの中心に置く、また、稜線の変数定義点をセル稜線の中心に定義し、更に、線分４の中心点の変数値を線形補間により求める、ことを特徴とする請求項３に記載の非圧縮性粘性流体の流れ場の数値解析方法。
- 物体に働く抗力（流れ方向の力）と揚力（流れの垂直方向の力）を、［数３］の式（１２）（１３）であらわす、
- 移動境界を伴う流体構造連成解析において、所定の時間刻み毎に、流体系と構造系を別々に解析し、それぞれの境界条件を陽的に受け渡す、ことを特徴とする請求項３に記載の非圧縮性粘性流体の流れ場の数値解析方法。
- 円柱の強制振動解析において、円柱が弾性支持された剛体構造物として、流れ方向に対して直角方向に振動する１質点１自由度系モデルと仮定し、円柱中心のＹ方向の変位を式（１７）で与え、円柱表面のＹ方向の速度境界条件を、式（１８）で与える、
ことを特徴とする請求項８に記載の非圧縮性粘性流体の流れ場の数値解析方法。 - 式（１８）から求められた振動する円柱の運動速度を、流れ場の速度境界条件として計算時間ステップごとに変化して与える、ことを特徴とする請求項９に記載の非圧縮性粘性流体の流れ場の数値解析方法。
- 円柱の渦放出による自励振動解析において、有次元の振動方程式を１質点１自由度のダンパ・ばねモデルとて、式（１９）又は式（２０）であらわす、
ことを特徴とする請求項８に記載の非圧縮性粘性流体の流れ場の数値解析方法。 - 式（２０）から求められた振動する円柱の運動速度を、流れ場の速度境界条件として計算時間ステップごとに変化して与える、ことを特徴とする請求項１１に記載の非圧縮性粘性流体の流れ場の数値解析方法。
- 円柱の初期変位と初期速度を０とし、かつ式（２０）の揚力は現在の値を用いて陽的に与え、式（２０）の振動方程式をニューマークのβ法により積分し、円柱の振動変位と振動速度を求める、ことを特徴とする請求項１１に記載の非圧縮性粘性流体の流れ場の数値解析方法。
- 非圧縮性粘性流体と接する対象物（１）の境界データからなる外部データ（１２）を入力する入力装置（２）と、形状と物理量を統合した実体データとその記憶演算プログラムを記憶する外部記憶装置（３）と、前記記憶プログラムを実行するための内部記憶装置（４）及び中央処理装置（５）と、実行結果を出力する出力装置（６）とを備え、
前記外部データを境界が直交する複数のセル（１３）に分割し、分割された各セルを対象物の内側に位置する内部セル（１３ａ）と境界データを含む境界セル（１３ｂ）とに区分し、前記境界データによる境界セル（１３ｂ）の稜線の切断点を求め、求めた切断点を結ぶ多角形を境界面のセル内部データとし、流れ場の境界に、ＶＯＦ法を併用したカットセル有限体積法を適用して解析する、ことを特徴とするＶ−ＣＡＤデータを直接用いた非圧縮性粘性流体の流れ場の数値解析装置。 - コンピュータに、非圧縮性粘性流体と接する対象物の境界データからなる外部データ（１２）を境界が直交する複数のセル（１３）に分割する分割ステップ（Ａ）と、
分割された各セルを対象物の内側に位置する内部セル（１３ａ）と境界データを含む境界セル（１３ｂ）とに区分するセル区分ステップ（Ｂ）と、
前記境界データによる境界セル（１３ｂ）の稜線の切断点を求める切断点決定ステップ（Ｃ）と、
求めた切断点を結ぶ多角形を境界面のセル内部データとする境界面決定ステップ（Ｄ）と、
流れ場の境界に、ＶＯＦ法を併用したカットセル有限体積法を適用して解析する解析ステップ（Ｅ）と、を実行させるための非圧縮性粘性流体の流れ場の数値解析プログラム。
Applications Claiming Priority (3)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| JP2002379214 | 2002-12-27 | ||
| JP2002379214 | 2002-12-27 | ||
| PCT/JP2003/015971 WO2004061723A1 (ja) | 2002-12-27 | 2003-12-12 | Ｖ−ｃａｄデータを直接用いた非圧縮性粘性流体の流れ場の数値解析方法と装置 |
Publications (2)
| Publication Number | Publication Date |
|---|---|
| JPWO2004061723A1 true JPWO2004061723A1 (ja) | 2006-05-18 |
| JP4442765B2 JP4442765B2 (ja) | 2010-03-31 |
Family
ID=32708381
Family Applications (1)
| Application Number | Title | Priority Date | Filing Date |
|---|---|---|---|
| JP2004564488A Expired - Fee Related JP4442765B2 (ja) | 2002-12-27 | 2003-12-12 | Ｖ−ｃａｄデータを直接用いた非圧縮性粘性流体の流れ場の数値解析方法と装置 |
Country Status (3)
| Country | Link |
|---|---|
| US (1) | US7430500B2 (ja) |
| JP (1) | JP4442765B2 (ja) |
| WO (1) | WO2004061723A1 (ja) |
Families Citing this family (39)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| US20040034304A1 (en) * | 2001-12-21 | 2004-02-19 | Chikayoshi Sumi | Displacement measurement method and apparatus, strain measurement method and apparatus elasticity and visco-elasticity constants measurement apparatus, and the elasticity and visco-elasticity constants measurement apparatus-based treatment apparatus |
| EP1574988B1 (en) * | 2004-03-08 | 2014-06-18 | Siemens Product Lifecycle Management Software Inc. | Determining and using geometric feature data |
| JP4783100B2 (ja) * | 2005-09-12 | 2011-09-28 | 独立行政法人理化学研究所 | 境界データのセル内形状データへの変換方法とその変換プログラム |
| US8010326B2 (en) | 2005-12-28 | 2011-08-30 | Caterpillar Inc. | Method and apparatus for automated grid formation in multi-cell system dynamics models |
| US20070150245A1 (en) * | 2005-12-28 | 2007-06-28 | Caterpillar Inc. | Method and apparatus for solving transport equations in multi-cell computer models of dynamic systems |
| US7542890B2 (en) | 2005-12-28 | 2009-06-02 | Convergent Thinking, Llc | Method and apparatus for implementing multi-grid computation for multi-cell computer models with embedded cells |
| US7590515B2 (en) * | 2005-12-28 | 2009-09-15 | Convergent Thinking, Llc | Method and apparatus for treating moving boundaries in multi-cell computer models of fluid dynamic systems |
| US7921002B2 (en) * | 2007-01-04 | 2011-04-05 | Honda Motor Co., Ltd. | Method and system for simulating flow of fluid around a body |
| US8326585B1 (en) * | 2007-06-18 | 2012-12-04 | The United States Of America As Represented By The Secretary Of The Army | Computational fluid dynamics shock wave identification |
| CN102160057B (zh) * | 2008-09-18 | 2014-07-16 | 国立大学法人京都大学 | 用于粒子法的界面粒子的判定方法及装置 |
| WO2010062710A1 (en) * | 2008-11-03 | 2010-06-03 | Saudi Arabian Oil Company | Three dimensional well block radius determiner machine and related computer implemented methods and program products |
| US8306798B2 (en) * | 2008-12-24 | 2012-11-06 | Intel Corporation | Fluid dynamics simulator |
| US8056424B2 (en) * | 2009-09-17 | 2011-11-15 | Sean P. Palacios | Multi-channel flow sensor with extended flow range and faster response |
| JP5277129B2 (ja) * | 2009-09-25 | 2013-08-28 | 株式会社日立製作所 | 流体解析装置 |
| ES2387170B1 (es) * | 2009-11-30 | 2013-08-20 | Airbus Operations S.L. | Metodos y sistemas para optimizar el diseño de superficies aerodinamicas |
| JP5718564B2 (ja) | 2009-12-04 | 2015-05-13 | トヨタ自動車株式会社 | 車両用駆動ユニットの数値解析方法 |
| JP5516949B2 (ja) * | 2009-12-25 | 2014-06-11 | 独立行政法人理化学研究所 | 数値流体計算方法、数値流体計算装置、プログラム、及び、記録媒体 |
| US8457939B2 (en) * | 2010-12-30 | 2013-06-04 | Aerion Corporation | Generating inviscid and viscous fluid-flow simulations over an aircraft surface using a fluid-flow mesh |
| JP5576309B2 (ja) | 2011-01-31 | 2014-08-20 | インターナショナル・ビジネス・マシーンズ・コーポレーション | 粒子法における自由表面に位置する粒子の正確な判定 |
| US8676552B2 (en) * | 2011-02-16 | 2014-03-18 | Adobe Systems Incorporated | Methods and apparatus for simulation of fluid motion using procedural shape growth |
| US8538738B2 (en) | 2011-03-22 | 2013-09-17 | Aerion Corporation | Predicting transition from laminar to turbulent flow over a surface |
| US8892408B2 (en) | 2011-03-23 | 2014-11-18 | Aerion Corporation | Generating inviscid and viscous fluid flow simulations over a surface using a quasi-simultaneous technique |
| US8917283B2 (en) | 2011-03-23 | 2014-12-23 | Adobe Systems Incorporated | Polygon processing techniques in procedural painting algorithms |
| US8744812B2 (en) * | 2011-05-27 | 2014-06-03 | International Business Machines Corporation | Computational fluid dynamics modeling of a bounded domain |
| WO2013063433A1 (en) * | 2011-10-26 | 2013-05-02 | Engine Simulation Partners | Method and apparatus for modeling interactions of the fluid with system boundaries in fluid dynamic systems |
| US9208269B2 (en) * | 2012-03-08 | 2015-12-08 | Engine Simulation Partners, Llc | Boundaries in fluid dynamic systems |
| JP6039527B2 (ja) * | 2013-10-18 | 2016-12-07 | 楽天株式会社 | 動画生成装置、動画生成方法及び動画生成プログラム |
| JP5464768B1 (ja) * | 2013-10-31 | 2014-04-09 | 株式会社ソフトウェアクレイドル | 情報解析装置、情報解析方法、及びプログラム |
| JP2017162269A (ja) * | 2016-03-10 | 2017-09-14 | ソニー株式会社 | 情報処理装置、電子機器、情報処理方法、及びプログラム |
| JP6576303B2 (ja) * | 2016-06-06 | 2019-09-18 | 住友化学株式会社 | シミュレーション装置、シミュレーション装置の制御方法、および制御プログラム |
| CN111651831B (zh) * | 2020-04-13 | 2022-04-08 | 北京航空航天大学 | 飞行器定常粘性可压缩绕流的分区扰动域更新计算方法 |
| CN111783277B (zh) * | 2020-06-04 | 2024-04-12 | 海仿(上海)科技有限公司 | 流体固体界面解耦算法、装置及设备 |
| CN111950174B (zh) * | 2020-07-08 | 2024-05-14 | 上海交通大学 | 流固耦合的计算方法、装置及电子设备 |
| CN112487610B (zh) * | 2020-11-09 | 2021-10-08 | 河北工业大学 | 具有复杂几何特征分析对象的形变确定方法及系统 |
| CN112800643B (zh) * | 2020-12-30 | 2024-03-08 | 新源动力股份有限公司 | 一种波纹流道燃料电池多物理场耦合计算简化方法 |
| AT525205A1 (de) * | 2021-07-07 | 2023-01-15 | Mehmet Bugra Akin | Verfahren zum Erstellen von Simulationszellen für kontinuumsmechanische Simulationen eines Objekts |
| CN114444348B (zh) * | 2021-12-31 | 2024-04-16 | 海油发展珠海管道工程有限公司 | 一种螺旋列板涡激振动抑制装置的动力学设计方法 |
| CN114638173B (zh) * | 2022-01-25 | 2023-06-02 | 中国空气动力研究与发展中心计算空气动力研究所 | 一种高阶非线性激波捕捉空间离散方法 |
| CN114925627B (zh) * | 2022-05-12 | 2024-03-15 | 南京航空航天大学 | 一种基于图形处理器的直升机流场数值模拟系统及方法 |
Family Cites Families (4)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| JP3308770B2 (ja) * | 1994-07-22 | 2002-07-29 | 三菱電機株式会社 | 情報処理装置および情報処理装置における計算方法 |
| US6161057A (en) * | 1995-07-28 | 2000-12-12 | Toray Industries, Inc. | Apparatus for analyzing a process of fluid flow, and a production method of an injection molded product |
| JP2000057127A (ja) * | 1998-08-11 | 2000-02-25 | Matsushita Electric Ind Co Ltd | 流体解析装置及び、プログラム記録媒体 |
| JP3468464B2 (ja) | 2001-02-01 | 2003-11-17 | 理化学研究所 | 形状と物性を統合したボリュームデータ生成方法 |
-
2003
- 2003-12-12 JP JP2004564488A patent/JP4442765B2/ja not_active Expired - Fee Related
- 2003-12-12 WO PCT/JP2003/015971 patent/WO2004061723A1/ja active Application Filing
- 2003-12-12 US US10/540,861 patent/US7430500B2/en not_active Expired - Fee Related
Also Published As
| Publication number | Publication date |
|---|---|
| US20060089803A1 (en) | 2006-04-27 |
| JP4442765B2 (ja) | 2010-03-31 |
| WO2004061723A1 (ja) | 2004-07-22 |
| US7430500B2 (en) | 2008-09-30 |
Similar Documents
| Publication | Publication Date | Title |
|---|---|---|
| JP4442765B2 (ja) | Ｖ−ｃａｄデータを直接用いた非圧縮性粘性流体の流れ場の数値解析方法と装置 | |
| Lomtev et al. | A discontinuous Galerkin ALE method for compressible viscous flows in moving domains | |
| Peller et al. | High‐order stable interpolations for immersed boundary methods | |
| Young Yoon et al. | A particle–gridless hybrid method for incompressible flows | |
| CA2853625C (en) | Computer simulation of physical processes | |
| Koshizuka et al. | Moving-particle semi-implicit method for fragmentation of incompressible fluid | |
| Burt et al. | Novel Cartesian implementation of the direct simulation Monte Carlo method | |
| Cole et al. | Validation of a commercial fluid-structure interaction solver with applications to air cushion vehicle flexible seals | |
| Zorrilla et al. | A modified finite element formulation for the imposition of the slip boundary condition over embedded volumeless geometries | |
| Ameur et al. | r-adaptive algorithms for supersonic flows with high-order Flux Reconstruction methods | |
| Launchbury | Unsteady turbulent flow modelling and applications | |
| Towara | Discrete adjoint optimization with OpenFOAM | |
| Nguyen | Material point method: basics and applications | |
| Ranocha et al. | On error-based step size control for discontinuous Galerkin methods for compressible fluid dynamics | |
| Frisani et al. | On the immersed boundary method: Finite element versus finite volume approach | |
| Müller | Multiresolution schemes for conservation laws | |
| Gornak et al. | A goal oriented survey on immersed boundary methods | |
| Chao et al. | A review on the applications of finite element method to heat transfer and fluid flow | |
| Apte et al. | A formulation for fully resolved simulation (FRS) of particle–turbulence interactions in two-phase flows | |
| Tai et al. | A finite volume unstructured multigrid method for efficient computation of unsteady incompressible viscous flows | |
| Donders et al. | Parameter uncertainty and variability in the structural dynamics modeling process | |
| Neusius et al. | A Cartesian cut-cell method for the isothermal compressible viscous Navier-Stokes Equations | |
| Venkatesan-Crome et al. | Discrete Adjoint for Unsteady Incompressible Flows Using a Density-based Formulation | |
| Montanari et al. | Isogeometric models for impact analysis with LS-DYNA | |
| Pattinson | A cut-cell, agglomerated-multigrid accelerated, Cartesian mesh method for compressible and incompressible flow |
Legal Events
| Date | Code | Title | Description |
|---|---|---|---|
| A621 | Written request for application examination |
Free format text: JAPANESE INTERMEDIATE CODE: A621 Effective date: 20061005 |
|
| A131 | Notification of reasons for refusal |
Free format text: JAPANESE INTERMEDIATE CODE: A131 Effective date: 20090821 |
|
| A521 | Request for written amendment filed |
Free format text: JAPANESE INTERMEDIATE CODE: A523 Effective date: 20091005 |
|
| TRDD | Decision of grant or rejection written | ||
| A01 | Written decision to grant a patent or to grant a registration (utility model) |
Free format text: JAPANESE INTERMEDIATE CODE: A01 Effective date: 20100106 |
|
| A01 | Written decision to grant a patent or to grant a registration (utility model) |
Free format text: JAPANESE INTERMEDIATE CODE: A01 |
|
| A61 | First payment of annual fees (during grant procedure) |
Free format text: JAPANESE INTERMEDIATE CODE: A61 Effective date: 20100106 |
|
| R150 | Certificate of patent or registration of utility model |
Free format text: JAPANESE INTERMEDIATE CODE: R150 |
|
| FPAY | Renewal fee payment (event date is renewal date of database) |
Free format text: PAYMENT UNTIL: 20130122 Year of fee payment: 3 |
|
| FPAY | Renewal fee payment (event date is renewal date of database) |
Free format text: PAYMENT UNTIL: 20130122 Year of fee payment: 3 |
|
| FPAY | Renewal fee payment (event date is renewal date of database) |
Free format text: PAYMENT UNTIL: 20160122 Year of fee payment: 6 |
|
| LAPS | Cancellation because of no payment of annual fees |