WO2004061723A1

WO2004061723A1 - Ｖ−ｃａｄデータを直接用いた非圧縮性粘性流体の流れ場の数値解析方法と装置

Info

Publication number: WO2004061723A1
Application number: PCT/JP2003/015971
Authority: WO
Inventors: Kangbin Lei; Masako Iwata; Ryutaro Himeno; Kiwamu Kase
Original assignee: Riken
Priority date: 2002-12-27
Filing date: 2003-12-12
Publication date: 2004-07-22
Also published as: US7430500B2; JP4442765B2; JPWO2004061723A1; US20060089803A1

Abstract

非圧縮性粘性流体と接する対象物の境界データからなる外部データ１２を境界が直交する複数のセル１３に分割する分割ステップ（Ａ）と、分割された各セルを対象物の内側に位置する内部セル１３ａと境界データを含む境界セル１３ｂとに区分するセル区分ステップ（Ｂ）と、境界データによる境界セル１３ｂの稜線の切断点を求める切断点決定ステップ（Ｃ）と、求めた切断点を結ぶ多角形を境界面のセル内部データとする境界面決定ステップＤと、流れ場の境界に、ＶＯＦ法を併用したカットセル有限体積法を適用して解析する解析ステップ（Ｅ）とを備える。

Description

明細書

V-C ADデータを直接用いた非圧縮性粘 ffi流の^れ場の数値,解析方法と装置 - 発明の背景

発明の技術分野

本発明は、形状と物理量を統合した実体データを記憶する V-CADデ一夕を' 直接用いた非圧縮性粘性流体の流れ場の数値解析方法と装置に関する。関連技術の説明

形状と物性を統合した実体データを小さい記憶容量で記憶することができ、これにより、物体の形状，構造 ·物性情報 '履歴を一元的に管理し、設計から加工、組立、試験、評価など一連の工程に関わるデータを同じデ一夕で管理することができ、 CADとシミュレーションを一元化することできる実体データの記憶方法として、 [特許文献 1] が開示されている。

【特許文献 1】

特開 2002— 230054号公報

[特許文献 1] の「形状と物性を統合した実体データの記憶方法」は、図 1に示すように、外部デ一夕入力ステップ（A) 、八分木分割ステップ（B) 、及びセルデータ記憶ステップ（C) からなり、外部デ一夕入力ステップ（A) では、外部データ取得ステップ S 1で取得した対象物の境界データからなる外部データ 1 2を本発明の方法を記憶したコンピュータ等に入力し、八分木分割ステップ (B) では、外部データ 1 2を八分木分割により境界平面が直交する直方体のセル 1 3に分割し、セルデータ記憶ステップ（C) では、各セル毎に種々の物性値を記憶するものである。

上述した [特許文献 1] の発明は、対象物の形状データからなる外部データを、八分木分割により境界平面が直交する直方体のセルに分割し、各セル毎に種々の物性値を記憶するものである。分割された各セルは対象物の内側に位置する内部セルと、境界面を含む境界セルとからなる。また、内部セルは、属性として 1種の物性値を持ち、境界セルは、対象物の内側と外側の 2種の物性値をもつものである。

以下、この方法によるデータを「V〜CADデータ」と呼び、このデータを用いた設計ゃシミュレーションを「ボリューム CAD」又は「V—CAD」と呼ぶ。図 1において 14が V— C ADデータである。

CFD (C omp u t a t i o n a l F l u i d Dyn am i c s ) が実用化するに従って格子生成に手間や時間がかかり、複雑形状では計算時間よりも格子生成の時間の方が長くなつてきている。このため、近年、直交格子による流体解析が再び話題となっている。直交格子による流体解析に関しては、 [非特許文献 1] 〜 [非特許文献 1 7] が知られている。

また、「強制振動円柱まわりの流れ」に関しては、 [非特許文献 1 8]に実験結果が開示され、「円柱の渦放出による自励振動」に関しては、 [非特許文献 1 9] に ALE有限要素法による計算結果が開示されている。

【非特許文献 1】

S a i k i , E. M. , B i r i n g e n, S . , 1 9 96, Nume r i c a 1 S i mu l a t i o n o f a Cy l i n d e r i n U n i f o r m f l ow : Ap p l i c a t i o n o f a V i r t u a l B o u n d a r y Me t h o d, J . C omp u t . P h y s . 1 23, 45 0-4 65.

【非特許文献 2】

矢部孝 · 肖鋒， 1 9 97，固体 ·液体 ·気体の統一解法と C I P法（2) , 数値流体力学， 7, 1 0 3-1 14.

【非特許文献 3】

Ye， T. , M i t t a l , R. , Ud a y kuma r, Η. S . , & S h vy， W.， 1 9 99， A C a r t e s i a n G r i d Me t h o d f o r V i s c o u s I n c omp r e s s i b l e F l ows w i t h C omp l e x I mme r s e d B o u n d a r i e s , A I A A- 9 9- 33 1 2, 545-5 5 7.

【弗特許文献 4】中村明 ·下村信雄 .里深信行， 1 99 5，デカルト格子系による任意形状物体周りの圧縮性粘性流計算，日本機械学会論文集， 6 1 B- 5 9 2, 4 3 1 9-4 3 26.

【非特許文献 5】

巿川治 ·藤井孝蔵， 20 02，直交格子を使用した三次元の任意形状物体まわりの流体シミュレーション，日本機械学会論文集， 6 8 B- 6 6 9, 1 3 2 9- 1 3 36.

【非特許文献 6】

朴炳湖 ·黒田成昭， 20 00，非圧縮性粘性流れの直交格子解法，ながれ， 1 9， 37-46.

【非特許文献 7】

〇 n o , K. , Tom i t a ， N. , F u j i t a n i K & H i m e n o , R. , 1 9 98， A n A p p 1 i c a t i o n o f V o X e 1 M

0 d e l i n g Ap p r o a c h t o P r e d i c t i o n o f En g i n e C o o 1 i n g F l ow, S o c i e t y o f A u t o m o t i v e E n g i n e e r s o f . J a P a n S P r i n g C o n v e n t i o n, No. 984 ， 1 6 5-1 68.

【非特許文献 8】

t t p ： / / kuwa h a r a. i s a s . a c . j p / i n d e x . h t m 1

【非特許文献 9】

寺本進 ·藤井孝蔵， 1 9 98，直交格子法による三次元物体周りの流れ解析，第 1 2回数値流体力学シンポジウム講演論文集， 29 9-300.

【非特許文献 1 0】

Qu i r k, J . J . , 1994, An A l t e r n a t i v e t o U n s t r u c t u r e d G r i d s f o r C omp u t i n g G a s Dyn am i c F l ows A r o u n d A r b i t r a r i l y C om p l e x Two - D i me n s i o n a l B o d i e s , C omp u t e r s F l u i d s , 23, 1 25-142. 【非特許文献 1 1】

Ka rma n, S. L. J r. ， 1 99 5, S PL I TF LOW： A 3D U n s t r u c t u r e d C a r t e s i a n/P r i sma t i c G r i d (1 2) y n a m i c s o f C F D C o d e f o r C omp l e x Ge ome- t r i e s , Α Ι ΑΑ 9 5-0343.

【非特許文献 12】

H i r t， C. W. , & N i c h o l s , B. D. , 1 9.8 1， V o 1 um e o f F l u i d (V〇F) Me t h o d f o r t h e D F r e e B o u n d a r i e s , J . C omp u t . P h y s . 3 9， 20 1-2 2 5.

【非特許文献 1 3】

H i r t， C. W. ， & C o o k, J . L. , 1 9 7 2, C a 1 c u 1 a t i n g Th r e e - d i me n s i o n a l F l ow s A r o u n d S t r u c t u r e s a n d Ov e r Ro u gh T e r r a i n, J . C omp u t . P h y s . 1 0， 324-340.

【非特許文献 14】

加瀨究 '手嶋吉法， 20 0 1, ボリューム CADの開発，理研シンポジウム ' ものつくり情報技術統合化研究，第 1回， 6- 1 1.

【非特許文献 1 5】

豊田郁夫 ·荒川忠一， 1 99 9, 直交適食格子法による円柱周り流れの解析，第 1 3回数値流体力学シンポジウム， F 03-1, CD-ROM.

【非特許文献 1 6】

松宮輝 ·木枝香織 ·谷口伸行 ·小林敏雄， 1 993，三次精度風上差分法による二次元円柱周り流れの数値シミュレーション，日本機械学会論文集， 5 9 B- 5 66, 2937-2943.

【非特許文献 1 7】

B o u a r d， R. , & C o u t a n c e a u， M. ， 1 980， Th e e a r l y s t a g e o f d e v e l o pme n t o f t h e w a k e b e h i n d a n i mp u l s i v e l y s t a r t e d c y 1 i n d e r f o r 40 <R e < 1 0 " 4 , J . F l u i d Me c h. , 1 0 1-3, 583-607 ·

【非特許文献 1 8】

O k a m o t o , S . ， Uema t s u, R . ， a n d T a g u w a , Y. ， F l u i d f o r c e a c t i n g o n two-d i me n s i o n a 1 c i r c u l a r c y l i n d e r i n L o k- i n p h e n om e n o n, J SME I n t e r n a t i o n a l ' J o u r n a l , B 45 , No. 4, (2002) , 8 5 0-856.

【非特許文献 1 9】

近藤典夫，円柱の空力挙動に関する数値シミュレーション，第 1 5回数値流体力学シンポジウム， E 09- 2, (200 1) ， CD-ROM. 流体解析では現在、複雑な三次元形状の流れ場も重合格子や非構造格子技法を使って計算できるようになったが、メッシュ生成がシミュレーシヨン全体の大きな部分を占めるようになった。このため、完全自動化できるメッシュ生成法として、直交格子を用いることが有望視されている。

直交格子系による任意形状の数値解析では、物体境界の取り扱いの難しさがよく知られている。近年、流れ場の境界近くの離散化や境界条件の扱う方法によつて、いくつかの直交格子法が提案されている。

例えば、仮想境界法 [非特許文献 1] (V i r t u a l B o un d a r y) , C I P [非特許文献 2] (Cu b i c-I n t e r p o l a t e d P r o p a g a t i o n) 密度関数法、 I mm e r s e d B o u n d a r y m e t h o d法 [非特許文献 3] , NPL C (Ne i g h b o r i n g P o i n t L o c a 1 C o l l o c a t i o n) 法 [非特許文献 4] ，格子点間に位置する境界までの距離を差分スキームに取り込む法 [非特許文献 5] 、部分境界適合直交格子法 [非特許文献 6] などがある。

これらの方法は物体境界が厳密に扱われるが、その分計算処理が複雑であり、必ずしも任意形状の三次元問題に向いているとはいえない。

一方、実用化の観点から、直交格子法は基本的に直交格子だけで階段状の境界を生成し、物体形状を近似する方法（例えば、小野 [非特許文献 7 ] · 日産自動車、桑原 [非特許文献 8 ] ·計算流体研）と、カットセルを導入して境界形状を取り扱う近似度を高める方法（例えば、藤井 [非特許文献 9 ] ·宇宙研、 Q u i r k [非特許文献 1 0 ] · J . J . ， N A S A) の二種類が有望である。

しかし、カットセルによる方法では、直交格子内を境界が任意の場所を通るため、境界上で隣り合ったセルの大きさに大きな差が生じることがあり、カツトセル直交格子では粘性流れの解析が難しいという報告 [非特許文献 1 1 ] もあった。上述したように、従来の重合格子や非構造格子を用いて、非圧縮粘性流体の流れ場の数値解析を行う場合には、格子生成を完全には自動化できず、そのため格子生成がシミュレーション時間全体に占める割合が高く、シミュレーション時間の短縮が困難な問題点があった。

一方、直交格子を用いた流れ場の数値解析は、格子生成を自動化できるものの、直交格子により物体境界を表現することが難しく、結果としてシミュレーション精度がわるい問題点があった。特に、移動境界に伴う流れの数値計算では、境界の移動距離が一定の大きさのメッシュの倍数に制約されるため、計算不安定に落ちる場合も多かった。

また特に、カットセルによる方法では、直交格子内を境界が任意の場所を通るため、境界上で隣り合ったセルの大きさに大きな差が生じることがあり、カットセル直交格子では粘性流れの解析が難しかった。発明の要約本発明は、かかる問題点を解決するために創案されたものである。すなわち、本発明の目的は、格子生成を完全に自動化することができ、物体境界の表現が容易であり、比較的簡単な計算処理により、計算不安定に落ちいることなく、精度の高いシミュレーションを短時間で行うことができる V- C A Dデータを直接用いた非圧縮性粘性流体の流れ場の数値解析方法と装置を提供することにある。本発明によれば、非圧縮性粘性流体と接する対象物の境界データからなる外部データ（ 1 2 )を境界が直交する複数のセル（ 1 3 )に分割する分割ステップ（A) と、分割された各セルを対象物の内側に位置する内部セル（13 a) と境界デー夕を含む境界セル（13 b) とに区分するセル区分ステップ（B) と、前記境界データによる境界セル（13 b) の稜線の切断点を求める切断点決定ステップ (C) と、求めた切断点を結ぶ多角形を境界面のセル内部データとする境界面決定ステップ（D) と、流れ場の境界に、 VOF法を併用したカットセル有限体積法を適用して解析する解析ステップ（E) と、を備えたことを特徴とする V - C A Dデータを直接用いた非圧縮性粘性流体の流れ場の数値解析方法が提供される。また、本発明によれば、非圧縮性粘性流体と接する対象物（1) の境界データからなる外部データ（12) を入力する入力装置（2) と、形状と物理量を統合した実体データとその記憶演算プログラムを記憶する外部記憶装置（3) と、前記記憶プログラムを実行するための内部記憶装置（4) 及び中央処理装置（5) と、実行結果を出力する出力装置（6) とを備え、

前記外部データを境界が直交する複数のセル（13) に分割し、分割された各セルを対象物の内側に位置する内部セル（13 a) と境界データを含む境界セル (13 b) とに区分し、前記境界データによる境界セル（13 b) の稜線の切断点を求め、求めた切断点を結ぶ多角形を境界面のセル内部データとし、流れ場の境界に、 V〇F法を併用したカットセル有限体積法を適用して解析する、ことを特徴とする V-CADデータを直接用いた非圧縮性粘性流体の流れ場の数値解析装置が提供される。

さらに、本発明によれば、コンピュータに、非圧縮性粘性流体と接する対象物の境界データからなる外部データ（12) を境界が直交する複数のセル（13) に分割する分割ステップ（A) と、分割された各セルを対象物の内側に位置する内部セル（13 a) と境界データを含む境界セル（13 b) とに区分するセル区分ステップ（B) と、前記境界デ一夕による境界セル（13 b) の稜線の切断点を求める切断点決定ステップ（C) と、求めた切断点を結ぶ多角形を境界面のセル内部データとする境界面決定ステップ（D) と、流れ場の境界に、 VOF法を併用したカットセル有限体積法を適用して解析する解析ステップ（E) と、を実行させるための非圧縮性粘性流体の流れ場の数値解析プログラムが提供される。上記本発明の方法及び装置により、精度と安定性を有し、計算コストがあまりかからず、かつ、格子生成を完全に自動化することができ、物体境界の表現が容易であり、比較的簡単な計算処理により、精度の高いシミュレーションを短時間で行うことができることが分かった。

本発明の好ましい実施形態によれば、前記解析ステップ（E) において、空間積分について、対流項に二次元の QU I CK補間スキームを使用し、拡散項に二次精度の中心差分を用い、時間進行について、対流項と拡散項を合わせて二次精度の Ad ams- B a s h f o r t h法を適用し、圧力勾配項に一次精度の E u 1 e r陰解法を用いる。

この方法により、格子生成を完全に自動化することができ、かつ安定性と離散化精度を確保することができる。

二次元境界セルにおいて、有限体積法における支配方程式を式（7) であらわす。

この式（7) は、非圧縮粘性流体の基礎支配方程式（1) をテンソル形の発散型に書き換えた式（6) を二次元境界セルの流体部分をコントロールボリューム (C V) Vi, 』として、空間積分したものであり、非圧縮粘性流体の基礎支配方程式（1) を満たすことができる。

有限体積法における支配方程式における対流項、圧力勾配項および拡散項を式 (8) 、（9) 、 (1 0) であらわす。

対流項：

= [Δ ( "

(8)

+ΓΑ ^ΐ ,, ^ ,,

+Ax(A_iJv^v_i - v v — , )]7 \only no - slip on wall

圧力勾配項：

= Bi Piu -Bi- jPi-v -P_P ^Bij -^Bi- j)f (9)

+Αχ[Α_υΡυ+υ2 - A _xP - p_p(A_UJ -4·,,-ι)]7

拡散項： ,

^div(grad(u))dV = grad(u) . MS = ^ grad (u)_m . nSS_m

y,j s,-₅

-B_{i j})grad{u)^x _p)

一 4, 4, 4, (10)

+Ax(A_{l j}grad(v)l_J+U2 - 4 (^v) 一 ~(4, ~

式（8) 〜（10) には、稜線の切断点を結ぶ対象物の境界データが含まれるため、境界における流れ場の非圧縮性粘性流体の数値解析ができる。

固体境界の積分において、固体境界でノースリップ境界条件の場合に、対流項は零とし、圧力勾配項と拡散項に対しては、切断線の中間点 Pの値を平均値として用いて積分し、空間積分に対してはすべての項に開口率を適用する。

この方法により、直交格子の完全自動化が保たれ、かつ流体計算におけるコントロールポリユームで保存則を厳密に満たすことができる。

カットセルにおいて算出される物理変数の定義点を、 V〇F = 0. 01を閾値としてそれより小さい境界セルは完全な固体とみなし、それより大きい境界セルにおいて算出される変数をセルの中心に置く、また、稜線の変数定義点をセル稜線の中心に定義し、更に、線分 4の中心点の変数値を線形補間により求める。

この方法によっても、直交格子の完全自動化が保たれ、かつ流体計算におけるコントロールポリユームで保存則を厳密に満たすことができる。

物体に働く抗カ（流れ方向の力）と揚力（流れの垂直方向の力）を、式（1 2 ) ( 1 3 ) であらわす。

抗カ：

^ Ji dx + 3y )

(12)

0；ひ） ~^σχχ ^v y Cartesian

'力

この式により、直交格子において、抗力と揚力を容易かつ精度よく求めるもとができる。

移動境界を伴う流体構造連成解析において、所定の時間刻み毎に、流体系と構造系を別々に解析し、それぞれの境界条件を陽的に受け渡す、ことが好ましい。

この解析方法により、流体系と構造系を連立して統一的に扱う強連成に比較して計算処理を簡略化でき、計算時間を短縮化できる。また境界の移動距離が一定の大きさのメッシュの倍数に制約されないので、計算不安定に落ちいることがほとんどない。

円柱の強制振動の流体構造連成解析において、円柱が弾性支持された剛体構造物として、流れ方向に対して直角方向に振動する 1質点 1自由度系モデルと仮定し、円柱中心の Y方向の変位を式（1 7) で与え、円柱表面の Y方向の速度境界条件を、式（1 8) で与える。

v„. = ^ cos(2^ (18)

また、式（1 8) から求められた振動する円柱の運動速度を、流れ場の速度境界条件として計算時間ステップごとに変化して与える。

この解析方法により、（1) 振動する円柱に作用する抗力および揚力が静止円柱のものより大きい、（2) S t r o u h a l数 0. 2付近の 0. 1 5く： f c < 0. 2 5において、強制振動による渦放出周波数のロックイン現象がよく再現される、（3) ロックイン領域での抗カゃ揚力が増大する、等、実験データと整合性の高い解析結果が確認された。

また円柱の渦放出による自励振動の流体構造連成解析において、、有次元の振動方程式を 1質点 1自由度のダンバ ·ばねモデルとて、式（1 9) 又は式（20) であらわす。 m^+c^+ ゆ (19)

dt² dt r ^{0 1} '

さらに式（20) から求められた振動する円柱の運動速度を、流れ場の速度境界条件として計算時間ステップごとに変化して与える。

また円柱の初期変位と初期速度を 0とし、かつ式（20) の揚力は現在の値を用いて陽的に与え、式（20) の振動方程式をニューマークの 3法により積分し、円柱の振動変位と振動速度を求める。

この解析方法により、（ 1 ) 無次元流速 V r = 4 (円柱の固有振動数 f o = 1/4= 0. 25) の場合に円柱の振幅がもっとも大きくなる、（2) 最大変位 yが、 V r = 2、 3、 4では yの最大変位は次第に大きくなり、 V r = 4で最大 1

12

になる、（3) その後、 V rが大きくなるにつれて、 yの最大変位は波状に小さくなつていく、等、 ALE有限要素法による計算結果と定性的に一致していることが確認された。

本発明のその他の目的及び有利な特徴は、添付図面を参照した以下の説明から明らかになろう。

図面の簡単な説明

図 1は、先行出願の実体データの記憶方法のフロー図である。

図 2は、本発明の数値解析方法を実行するための数値解析装置の構成図である。図 3は、本発明の数値解析方法とそのプログラムのフロー図である。

図 4 Aと図 4 Bは、二次元における開口率の定義を示す図である。

図 5'は、境界セルにおける有限体積法による空間離散化の説明図である。

図 6A, B, Cは、解析した 3ケースの VOF分布を示す画像である。

図 7A, B， Cは、解析した 3ケースの流れ方向速度の分布と理論解析解との比較図である。

図 8は、解析格子と VOF分布を示す画像である。

図 9は、速度ベクトルを示す画像である。

図 1 0は、圧力等値線を示す画像である。

図 1 1は、本発明による平均抗カ係数を示す図である。

図 1 2は、本発明による平均 S t r 0 u h a 1数を示す図である。

図 1 3Aは、直交適合格子法による静止円柱の抗カ係数を示す図であり、図 1 3 Bは、その平均 S t r o u h a 1数を示す図である。

図 14Aは、一般座標系境界適合格子による静止円柱の抗カ係数を示す図であり、図 14 Bはその平均 S t r o u h a 1数を示す図である。

図 1 5 Aは、 R e = 30 0， T= 2. 5の場合の実験による流脈線を示す画像であり、図 1 5 Βはその流線を示す画像である。

図 1 6 Αは、 Re = 55 0， T= 2. 5の場合の実験による流脈線を示す画像であり、図 1 6 Bは、その流線を示す画像である。

図 1 7 Aは、抗力と揚力の履歴変化を示す本発明の結果であり、図 1 7 Bは、その部分境界適合直交格子法の結果である。

図 1 8A〜Fは、実用問題に向けての応用例を示す画像である。

図 1 9A〜I は、静止円柱と各振動数での円柱に作用する流体力の抗カ係数 C _Dおよび揚力係数の時間履歴の変動波形を示す図である。

図 20は、流れの渦放出状態を表す S t r o u h a 1数 S tと強制振動周波数 f cの関係図である。

図 2 1 A〜(：は、無次元流速 Vrが 2, 3, 4における変位 y (左図）と抗カ係数 C_D及び揚力係数 (右図）の時間履歴を示す図である。

図 22D〜G は、無次元流速 Vr が 5， 6, 7， 8における変位 y (左図）と抗カ係数 C_D及び揚力係数 (右図）の時間履歴を示す図である。

図 23は、各々の無次元流速 V rに対する yの最大変位を示す図である。

好ましい実施例の説明以下、本発明の好ましい実施形態を図面を参照して説明する。

図 2は、本発明の数値解析方法を実行するための数値解析装置の構成図である。この図に示すように、本発明の数値解析装置 1 0は、入力装置 2、外部記憶装置 3、内部記憶装置 4、中央処理装置 5および出力装置 6を備える。

入力装置 2は、例えばキーボードであり、対象物 1の形状デ一夕からなる外部データ 1 2を入力する。外部記憶装置 3は、ハードディスク、フロッピィ一ディスク、磁気テープ、コンパクトディスク等であり、形状と物理量を統合した実体データとその記憶演算プログラムを記憶する。内部記憶装置 4は、例えば RAM, ROM等であり、演算情報を保管する。中央処理装置 5 (CPU) は、演算や入出力等を集中的に処理し、内部記憶装置 4と共に、記憶プログラムを実行する。出力装置 6は、例えば表示装置とプリン夕であり、記憶した実体データと記憶プログラムの実行結果を出力するようになっている。本発明の記憶演算装置 1 0は、上述した外部記憶装置 3、内部記憶装置 4、及ぴ中央処理装置 5により、外部データを境界が直交する複数のセル 1 3に分割し、分割された各セルを対象物の内側に位置する内部セル 1 3 aと境界データを含む境界セル 1 3 bとに区分し、境界データによる境界セル 1 3 bの稜線の切断点を求め、求めた切断点を結ぶ多角形を境界面のセル内部デ一夕とし、流れ場の境界に、 VOF法を併用したカツトセル有限体積法を適用して解析する。

図 3は、本発明の数値解析方法とそのプログラムのフロー図である。この図に示すように、本発明の方法及び変換プログラムは、分割ステップ（A) 、セル区分ステップ（B) 、切断点決定ステップ（C) 、境界面決定ステップ（D) 、及び解析ステップ（E) からなる。

外部から入力する外部データ 1 2は、多面体を表すポリゴンデ一夕、有限要素法に用いる四面体又は六面体要素、 3次元 CAD又は CGツールに用いる曲面デ —夕、或いはその他の立体の表面を部分的な平面や曲面で構成された情報で表現するデータである。

外部データ 1 2は、このようなデータ（S— CADデータと呼ぶ）のほかに、 ( 1) V— CAD独自のインターフェース（V— i n t e r f a c e ) により人間の入力により直接作成されたデータと、（2) 測定機やセンサ、デジタイザなどの表面のデジタイズデ一夕や、（3) CTスキャンや MR I、および一般的に V o 1 umeレンダリングに用いられているポクセルデータなどの内部情報ももつ V ο 1 um eデータであってもよい。

分割ステップ（A) では、外部データ取得ステップ（図示せず）で取得した非圧縮性粘性流体と接する対象物の境界データからなる外部データ 1 2を境界平面が直交する複数のセル 1 3に分割する。この分割は、 3次元の場合は、八分木分割であり、 2次元の場合には 4分割である。

すなわちこの分割ステップ（A) における分割とは、目的とする非圧縮性粘性流体と接する対象物を含む、基準となる直方体（または矩形）を分割（8分割又は 4分割）し、それぞれの領域の中に立体が完全に含まれるか、含まれなくなるまで再帰的に分割処理を繰り返す。この分割によりポクセル表現よりも大幅にデ一夕量を減らすことができる。空間分割により分割された一つの空間領域をセル 1 3とよぶ。セルは境界が直交する直方体または矩形である。セルによる階層構造、分割数もしくは分解能によって空間中に占める領域を表現する。これにより空間全体の中で対象は大きさの異なるセルを積み重ねたものとして表現される。

セル区分ステップ（B) では、分割された各セルを対象物の内側に位置する内部セル 1 3 aと境界データを含む境界セル 1 3 bとに区分する。

すなわち本発明では非圧縮性粘性流体と接する対象物の内側または外側に完全に内部に含まれるものはその最大の大きさをもつ内部セル 1 3 a (立方体）とし、外部データ 1 2からの境界情報を含むセルは境界セル 1 3 bとする。

切断点決定ステップ（C) では、境界データによる境界セル 1 3 bの稜線の切断点を求める。

境界面決定ステップ（D) では、求めた切断点を結ぶ多角形を境界面のセル内部データとする。以下、かかる切断点を結ぶ多角形を含むセルを「カットセル」と呼ぶ。 '

解析ステップ（E) では、上述した内部セル 1 3 aと境界セル 1 3 bに対して、流れ場の境界に、 VOF法を併用したカツトセル有限体積法を適用して非圧縮性粘性流体の流れ場の数値解析を行う。以下、本発明を更に詳細に説明する。

1. 本発明では V-CADプログラムにおける実用性のある流体解析技術を目指し、カットセル（KTC) 直交格子による任意形状の非圧縮粘性流れの解析方法を開発した。本発明において、流れ場の境界での扱いについては、 H i r tらが提案している VOF (Vo l ume O f F l u i d) 法 [非特許文献 12] を併用したカツトセル有限体積法を用いた。また本発明の方法による内部流であるチャンネル内流れと、外部流の円柱まわり流れを解析し、実験データ、理論解、及び既存の方法による解析結果との比較を行った。更に、本解析法の応用例として、いくつかの計算結果を合わせて示した。

2. 基礎方程式と計算方法 2. 1 基礎方程式

本発明で用いる基礎支配方程式は [数 4] の式（1) (2) に示す非圧縮粘性流体の N a v i e r - S t o k e s方程式と連続の式である。

【数 4】 (1)

^L=0 (2)

ここで R eは流れ場の.代表長さと代表速度で定義されるレイノルズ数と呼ばれる無次元数であり、物理的には流れ場における慣性力と粘性力の比を表している。 uは速度、 pは圧力である。また、 i = l， 2, 3、 j = 1 , 2， 3は直交座標系での各方向を表している。なお、 i , j に関しては縮約を取ることとして、これ以降もそのように表すものとする。式（ 1) より流れ場における外力 f iを考慮しない場合に非圧縮性粘性流体の運動はただ一つのパラメータ. R eで支配されていることが分かる。

2. 2 計算精度 ·計算コスト及び計算方法

流体のシミュレーションは一種の数値実験であり、つねに一定の誤差を伴う。流体解析を精度よく行うには、以下の四つの条件が必要となる。

① 流れの最小長さスケール（境界層 ·渦 ·衝撃波 ·火炎面など）を捕られるほど細かい空間解像度。

② 流れの最大長さスケールを十分に捕獲するほど大きい計算領域と、人為的な流入 ·流出境界条件 ·壁によるブロッキングなど影響を無視できるほどの計算領域。

③ 打ち切り誤差や数値拡散などを無視できるほど十分な空間的 ·時間的な離散化精度。

④ 特定問題に即したモデル（壁面モデル ·乱流 S GSモデル ·燃焼モデルなど）やスキーム（K— K ' QU I CK ' MUS CL ' TVD ' ENOなど）、差分法を援用すれば、できるだけ直交等方性のある格子。

流体解析の計算コスト（計算時間 '必要なメモリ）は、その計算方法に求められる解析精度によって決められるといえる。実用的な流体解析では、ある程度の精度（例えば、普通の物理実験の精度範囲や、ユーザーに求められる精度）があれば、十分と考えられる。

一方、流れ場は内部流、外部流及びそのほか（例えば噴流）に分類され.る。一般に外部流の解析領域は内部流のものより大きくとるのが普通である。本発明では実用性の観点から、外部流の場合には、流れ場の代表長さスケール D (= 1) として、解析領域の大きさはすべて 1 0 D X 1 0 Dとした。

本発明では、差分法を併用した有限体積法を用いる。ある程度の安定性や離散化精度を確保するため、空間積分について、対流項に二次元の QU I CK補間スキ一ムを使用し、拡散項に二次精度の中心差分を用いる。時間進行について、対流項と拡散項を合わせて二次精度の A d ams-B a s h f o r t h法を適用し、圧力勾配項に一次精度の Eu 1 e r陰解法を用いる。解析アルゴリズムとしての圧力と速度のカップリングには、 H i r tらにより提案された圧力と速度を同時緩和する S〇L A- HSMAC法 [非特許文献 1 3] を使用し、緩和係数 1. 6 5で SOR法を用いて反復計算を行う。収束判定は連続方程式（2) の残差が 0. 0002とする。また圧力振動を防ぐために三方向の速度 u.， V, w, 及び圧力 pの格子上の定義点を、半メッシュをずれるスタガード格子を配置する。

2. 3 V- CADデータによる固体 ·流体境界の取扱い

V o 1 ume-CADシステムの中でもっとも重要なアルゴリズムは加瀬 ·手嶋 [非特許文献 14] により提案されるセルの切断点からモノの表面形状を復元する KTCアルゴリズムである。この KTCアルゴリズムは流体解析分野における直交格子カツトセル（Cu t C e l l ) 法と同じ概念であることを注意されたい。また、 KT Cは二次元においても一個のセルの四つの稜線上に 2個以上の交点（切断点）があり得るが、本発明では簡単のために二次元で一個のセルの稜線上に 2個だけの切断点に限るとする。この場合の K T Cは図 4 Aと図 4 Bに示すように二種類しかないことがわかる。

カツトセル法における開口率は、セルの稜線に流体部分が占める比率として図 4 Aと図 4 Bに示すように流体 ·固体境界とセル稜線との交点情報から [数 5] の式（3) で定義される。

ここで、 AX i, j、し jはそれぞれ X方向と Y方向の格子幅であり、 Ai, j、 B jは X方向と Y方向の開口率である。また、 Vx ^ j Vy _{i ;} jはそのセル稜線上に固体が占める線分を表す。

各方向の開口率が分かれば、 VOF法における流体の体積占有率は [数 5] の式（4) と式（5) により求められる。

【数 5】ん . 、厂 _β ^lZ^l o<A_{i Bij}<X (3)

VF_U -1 - 0.5(2- -A )(l-B -B.__j) ( )

when{A_u + A_i > ^ nd^ + > 1)

ゾ=0.5(ん十 — 'X^+U (5)

when(A.j + A < l)or(5_{i y}.十 Β'.― ≤ 1)

2. 4 VOF有限体積法による離散化

有限体積法で保存型の支配方程式を離散化する場合、方程式（1) の積分形をグリーン定理（発散定理）に適用するため、微分型の方程式（1) をテンソル形の発散型に書き換えると [数 6 ] の式（6 ) のようになる。

ここで、〇に Xはテンソル積であり、ベクトル Iは K r o n e c k e r <5 _u に対応する単位テンソルである。

流体固体境界での離散化は、本発明の肝心なところになる。境界セル以外の格子における有限体積法による直交格子での離散化は従来の方法でもあまり問題が生じないのでここでは省略する。流体固体境界では図 4のような KTCセルを想定して式（6) の空間離散化を行う。また式（6) の時間離散化について、流体占有率 VOFを用いてコントロールポリューム（CV) 内における速度の時間勾配を平均値とみなして、従来の方法通りで行えるため、境界セルでの時間離散化の説明も省く。さて、図 5の二次元境界セルの流体部分をコントロールボリューム（CV) _{i t} jとして、式（6) を空間積分する（ここで簡単のため外力項べクトル f = とした）と、有限体積法における支配方程式は [数 6] の式（7) になる。

【数 6】

― - - ν(ΰ Θΰ)- ν(ρΐ) +— div(grad(u))+ f ( 6 )

dt Re

式（7) における対流項、圧力勾配項および拡散項をそれぞれグリーン定理（発散定理）に適用すると、コントロールポリュームにおける面積分（三次元の場合に体積積分）は、閉直線の線積分（三次元の場合に閉平面の面積分）に転換され、図 5のコントロールボリュームの五つの稜線（1→2→3→4→5→1) に（反時計まわり方向をとる）沿って線積分することになる。本発明のデカルト座標では、上記の対流項、圧力勾配項および拡散項はそれぞれ [数 7] のように離散化される。

【数 7】

対流項：

=[Ay¾ ^>«_/J-5_i__1J^^-w)

+Δ(4. ^;),_+1/2V" - . -i«M/2.y-_I/2V,,_y_₁)]? ( 8 )

•HA A. - ,；： , ,ゾ）

+Δ(4,ゾ vg)v"一

on wall

圧力勾配項：

= Ay[B_{i Pi},_V2

( 9 )

+Ax[A_i p_{i +l/2} - A_u__lPi _ 一 p_p(A 一 -i )V 拡散項：

grad{u)_m · n$S_m

= [Ay(B grad(u) _{U2 j}

— ;) );）

+Ax(A_{i i}gmd( _j+V2 - A,_Hg^md(^ll)ij- 2 - (4,ゾ― Aト、 )g^md(^ll)^} _P )f (10)

ここで、速度変数の下添え字はスタガード格子ステンシルにより変化することを注意されたい。またはベクトル i, jはそれぞれ X方向と Y方向の方程式を表す。対流項の離散式（8) の中に上添え字（X) (y) が付いた速度は、それぞれ X方向と Y方向に二次元 QU I CK補間を適用する。拡散項の離散式（10) の中に上添え字 x、 yが付いた速度勾配は、それぞれ X方向と Y方向の速度勾配を表し、そこで各々二次中心差分を適用する。

カツトセルによる離散化には、二つのキ一ボイントがあることを注意されたい。キーポイントの一つは固体境界（図 5の切断線 4) の積分の取り扱いである。有限体積法では、一階偏微分の対流項と圧力勾配項の境界条件はディリクレ条件であり、二階偏微分の拡散項の境界条件はノイマン問題となる。固体境界でノースリップ境界条件の場合に、図 5の切断線 4のところに対流項が零となるから対流項を積分しなくてもいいが、圧力勾配項と拡散項に対して切断線 4のところに積分しなければいけない。これは式（9 ) と式（1 0 ) の各々最後の項に切断線 4の中間点 Pの値を（平均値として）用いて積分するわけである。また空間積分に対してすべての項に開口率を適用する必要もあることを注意されたい。

もう一つのキ一ボイントは、カツトセルにおいて算出される物理変数の定義点である。有限体積法では流体占有率 V O Fで体積平均を行うため、物理変数の定義点は流体のコントロールポリュームの中心（例えば図 5の〇に ·点で示す C点）に置くのが自然だが、そうすると直交等間隔格子の完全自動化メリットがなくなる。その故、本発明では V O F = 0 . 0 1を閾値としてそれより小さい境界セルは完全の固体とみなし、それより大きい境界セルにおいて算出される変数はコン卜ロールポリュームの中心ではなく、セルの中心（例えば図 5の〇に +印で示す D点）に置くこととする。同じく図 5の稜線 3と稜線 5の変数定義点も線分の中心〇ではなく、セル稜線の中心 ©に定義される。また線分 4の中心点 Pに変数が定義されていないため、本発明でその点の変数値を線形補間により求める。こう処理すると、固体境界の再現性を損するものの、直交格子の完全自動化が保たれ、実用化に向ける。また流体の計算はコント口一ルポリュームで保存則が厳密に満たされ、上記の処理は流れ場の計算精度に及ぼす影響がそれほど大きくないと考えられる。

2 . 5 物体に働く坊力と揚力について

物体まわりの流れを解析するとき、その物体に働く抗力と揚力の計算について、一般座標系の境界適合格子の場合に物体の第 1格子点に沿って積分すれば簡単だが、直交格子で取り扱う場合にはやや煩雑となる。そこで本発明では、以下の方法で抗力と揚力を求める。

物体に働く力は、その物体表面に作用する流体応力テンソルを閉積分して、グリーン定理を適用すれば、二次元の場合に [数 8 ] の式（1 1 ) ように得られる。ここで σ _nは法線方向の平面に働く応力テンソル、 σは物体表面の点での応力テンソルである。また二次元の二階応力テンソルの発散はそのテンソルの四つの方向からなることを注意されたい。それで X方向を流れ方向、 Y方向を流れの垂直方向とすると、物体に働く抗カ（流れ方向の力）と揚力（流れの垂直方向の力）は、グリーン定理を再び適用すると [数 8] の式（1 2) (1 3) のようになる。抗力と揚力が分かれば、無次元量としての抗カ係数 C_Dと揚力係数 C_Lを、主流速度 U、流体密度 ιθ及び物体の基準面積 A (二次元の場合は基準長さ）を用いて [数 5] の式（14) のように無次元化される。

また非定常現象を表わす無次元量のストローハル数 S tは、 C_Lの周波数： f _L から [数 8] の式（1 5) で求める。

301S971

23

【数 8】 f 卡び

(11)

3σ„ び„

+ ^-) dxdy + ( +―^) dxdy =F +F_yj

dx dy dx dy 抗カ

du

cr 二び = ^=-ρ + 2μ— (14)

dx 、¾ 3

C, = (15)

^D pU^zD/2 ^L~ pU²D/2

St- f_LDlU (16)

3. 計算結果および考察

流体の数値解析において解析スキームや数値モデルなどを検証することにあたつて、内部流のベンチマークテストとしてチャネル流れと、外部流のベンチマークテストとしての円柱まわりの流れはしばしば検証の対象となる。本発明でも、上記の V OF法を併用したカツトセル有限体積法を検証するために、解析理論解のある二次元 P o i s e u i 1 1 e流と、多くの実験や計算データのある静止円柱まわりの流れについて数値解析を行い、他の研究者のデータとの比較を行った。なお、実用問題に向けての応用例として、パックステップを過ぎる流れ、狭窄管内の流れ、チャネル内障害物周りの流れ、分岐管内の流れ、複数 B l u f f B o d yまわりの流れゃ予混合燃焼器内の流れなどについて解析を行い、その計算結果の可視化例を示す。

3. 1 チャネル P 0 i s e u i 1 1 e流による検証

セルのカットされた V〇Fの効果を検証するため、 1 0DX 1 0D正方形の解析領域におけるダクトを、 0度、 1 0度及び 45度を斜めにして解析を行った。全体の格子数は 64 X 64、ダクトの半幅は 2. 5D、この半幅で定義されるレィノルズ数 R e = 1、時間刻み△ t = 0. 000 1、解析時間 s t e p s = 1 0 0 0 0とする。ダクトの流入条件と流出条件はともに放物型の分布を仮定した。流れ方向の圧力勾配は 1. 0とした。また HSMAC法を用い、入口 ' 出口および壁での圧力境界条件は直接与えていない。三ケースの V〇F分布を図 6 A， B, Cに示す。図 6A， B， Cにおける破線は、境界を示している。その三ケースに対応する流れ方向速度分布を理論解析解と比較して図 7 A， B, Cに示す。

図 7A, B, Cの中に線 Aの a n a 1 y s i sは理論解析解であり、線 Bの i = 32は第 32格子点断面にある流れ軸中心での流れ垂直断面の速度分布である。ダクトを傾けない（角度 0) 計算結果の速度分布は解析解とよく一致し、傾斜角度を 1 0度、 45度にすると、軸中心の最大速度がわずかに過小評価されるものの、全体として解析結果は理論解と一致していると言える。

3. 2 静止円柱まわり流れによる検証

解析スキームやモデルなどの検証対象となる円柱まわり流れの解析は、直交格子にとって最も困難な例題であろう。本発明では、 3. 1の解析領域（ 1 0 DX 1 0 D) と格子数（64 X 64) のままで、円柱直径 Dで定義されるレイノルズ数の R e = lから R e = 20000までの 24ケースについて、一様流中に置かれている静止円柱まわり流れの解析を行った。計算の無次元時間は 1 0 0とし、円柱の円心座標の図 6のように X=2. 5D， Y= 5. 0Dとした。速度の境界条件として、円柱表面で滑り無し条件、流入側の流入条件は一様流（流れ方向速度 =1. 0) 条件、流出側の流出条件は X方向に自由流出（速度勾配 0) ， Y方向にフリースリップ条件とした。また、 HSMAC法を用いて速度 ·圧力を同時に緩和するため、流入 ·流出および壁での圧力境界条件を直接与える必要がない。図 8は解析格子と解析全領域の VOF分布を示す。図 9、図 10はレイノルズ数 Re = 300の場合の、速度ベクトルと圧力等値線であり、円柱後方に周期的なカルマン渦列の形成が観察されており、非定常流の流れパタン一の特徴を捉えている。これはレイノルズ数 R e>80の場合に共通なの.で、他のレイノルズ数における計算結果は省略した。

図 1 1、図 12にそれぞれ実験データなどと比較した無次元時間 50〜100 にわた'る抗カ係数の平均値と、 '平均 S t r ouh a 1数を示す。抗カ係数について、低レイノルズ数（Re<60) の本発明による計算結果は今井の式による結果（B) とよく一致している。レイノルズ数 Re = l 00付近あたりにおいて、本発明による計算結果は実験データ（C)と比較してわずかに小さくなつている。レイノルズ数 Re>400では、本発明による計算結果は実験データ（C) とほとんど相違がない結果となっている。

これは主に本発明で用いる空間解像度（10/64=0. 15) が低いので、 Re = l 00 (0) オーダー以上になると、円柱壁における境界層は十分に解像されていないと考えられる。図 11から分かるようにレイノルズ数 R e>50では、カットセル無しの階段近似による計算'結果（D) は実験より明らかに過大評価される.。これに対して、本発明で提案したカットセル VOF法による計算は、計算結果（D) に比べて大幅に改善されていることが分かる。

図 12の平均 S t r o u h a 1数について、本発明での計算結果は、実験デー夕より大きめであるが、レイノルズ数 Re> 200の適正値とされる.0. 2〜0. 25の範囲に収まっていることがわかる。

図 13 A, Bに他の研究者による直交適合格子法 [非特許文献 15] を使用した円柱の抗カ係数と S t r o u h a 1数の解析結果を示す。また図 14A, Bに一般座標系適合格子法 [非特許文献 16] による抗カ係数と S t r ouh a 1数の解析結果を示す。これら従来の方法による結果では、レイノルズ数 Re：> 20 0の適正値とされる 0. 2〜0. 25の範囲から部分的に外れていることがわかる。従って本発明による平均 S t r o u h a 1数の計算結果（図 11、図 12) は、これら従来の境界適合格子による計算結果と比べて計算精度が優れていることが分かる。

レイノルズ数 Re = 300の場合に無次元時間 T= 2.5の時の可視化実験 [非特許文献 17] の結果と、本発明での同じ条件による計算結果の円柱近傍拡大図を図 15A., Βに示す。またレイノルズ数 Re = 550の場合（無次元時間 T =

2. 5) に実験 [非特許文献 17] と照する円柱付近の計算結果を図 16 Α， Βに示す。

上記の図からかるようにいずれのケースで円柱付近の双子渦の形状や大きさ、円柱直後から渦の循環再合流までの距離など、本発明での計算結果は可視化実験とよく一致している。

図 1.7 Α， Βにレイノルズ数 Re = 200場合の、本発明の方法による円柱抗力と揚力の時間履歴変化を、朴ら [非特許文献 6] により提案された部分境界適合直交格子法による計算結果と比較して示す。この図から、本発明の等間隔直交格子による抗カおよび S t r ouha 1数は、部分境界適合直交格子法による結果とよく一致して.いることが確認された。なお、本発明の揚力の変動幅は部分境界適合直交格子法による結果より小さく評価されているが、これは主に [非特許文献 6] の計算領域は 30DX 15Dであり、流出境界として放射条件を使用しているためと考えられる。

ここで図に示されていないが、本発明でレイノルズ数 R e.= 1から R e = 20 000までの 24ケース計算を行ったが、レイノルズ数 80以上でカルマン渦が出るようになった。実験や理論の臨界レイノルズ数 Re = 50位に比べてやや大きめであるが、実用上十分であると思われる。

3. 3 実用問題に向けての応用例

今回開発した解析方法をいくつかの例に応用してみた。それらは図 18A〜F のようにパックステップを'過ぎる流れ（A) 、狭窄管内の流れ（B) 、チャネル内障害物周りの流れ（C) 、分岐管内の流れ（D) 、複数 B l u f f Body まわりの流れ（E) ゃ予混合燃焼器内の流れ（F) などである。

ここに示しているのは X方向速度 u分布のみであり、これらの応用例の計算は 3. 1のダクト計算や 3. 2の円柱まわりの計算と同じように計算領域は 1 0D X 1 0 Dで、格子数は 64 X 64とした。またレイノルズ数は R e = 1 0 0、流入条件や流出条件等は 3. 2の計算と同じにした。つまり、開発した計算プログラムを変更することなく、 V-C ADデータを直接に入力するだけで、任意形状の流れ場を計算できるようになった。

4. 次に本発明の解析方法を移動境界に伴う流れの数値計算（流体構造連成解析）に適用する。流体構造連成解析は、風工学や熱交換器など工学的見地において特に重要である。

流体構造連成解析の具体例として、（1) 一様流の中に置かれた二次元円柱の強制振動と、（2) 円柱の渦自励振動（渦励振）の数値計算を行い、構造物に致命的な損傷を与える代表的な流体構造連成問題としてのロックイン現象や渦励振による共振状態を検証し、移動境界の流体構造連成解析における本発明の解析方法の有効性を示す。

4. 1 基礎方程式

円柱の振動計算に関しては、弾性支持された剛体構造物として、流れ方向に対して直角方向に振動する 1質点 1自由度系モデルを仮定している。

円柱が強制振動の場合には、円柱を一様流と直角な方向（Y方向）に s i n波で強制的に振動させる。即ち円柱中心の Y方向の変位を [数 9]の式（1 7) で与える。ここで、 Aは振幅、 fcは強制振動周波数、 tは時間である。

また、円柱表面の Y方向の速度境界条件は、 [数 9]の式（1 8) で与えられる。円柱が渦放出による自励振動をする場合には、有次元の振動方程式は 1質点 1 自由度のダンパ · ばねモデルによって [数 9]の式（ 1 9) に示す。

式（1 9) の中でそれぞれ mは円柱の質量、 cはダンパの粘性減衰係数、 kはパネの弾性回復係数、 pは流体の密度、 U。は一様流の流速、 Dは円柱の直径、 C ま揚力係数である。式（1 9) を無次元化して円柱の固有振動数 f 。= (k/ m) °' ⁵/ (2 κ) 、無次元減衰係数 h = cZ (2 (mk) 。· ⁵) 、スクルートン数 S c = 27t h (2 p _e/p) などの振動応答を支配するパラメータで整理すると、 [数 9]の式（2 0) になる。

【数 9】 y = Asi„ ) (17) v„, = ^ _ccos(2^ (¹⁸)

dt² dt 2 " ^L ^ ^d ,hf_n) H^f_oy_y^C_L (20)

dt at Sc

4. 2 数値解法

円柱が振動する場合には、式（1 8) または式（2 0) から求められた振動する円柱の運動速度を、流れ場の速度境界条件として計算時間ステップごとに変化して与えている。

円柱の渦自励振動の解析については、円柱の初期変位と初期速度を 0として、式（2 0) の振動方程式をニューマークの /3法により積分し、円柱の振動変位や振動速度を求められる。

本発明の流体構造連成解析は、流体系と構造系を連立して統一的に扱う強連成ではなく、流体と構造の方程式を別々に離散化し、それぞれの境界条件を陽的に受け渡す弱連成となる。複雑な流体構造系の連成解析における時間刻みは、特に強連成の場合、流体または構造計算の時間刻みの小さい方に制約されることは原則であるが、本発明での構造系の計算には単純なダンパ ·パネモデルを用いるため、計算の時間刻みは流体解析の CFL条件のみから決められる。また、ニューマークの i3法により将来の加速度を積分する際に、式（2 0) の揚力は現在の値を用いて陽的に与える。

4. 3 強制振動円柱まわりの流れ

流れ場のレイノルズ数 R e = 1 0 0 0の場合には、 3. 2の解析条件のままで、 1

29

式（1 7) と式（1 8) による強制振動円柱まわりの流れ解析を行った。レイノ. ルズ数 R e = 1 000の際に、静止円柱の無次元渦放出周波数の S t r o u h a 1数は S t = 0. 2位であることが実験値で知られている。強制振動の計算条件として振幅 Aを 0. 1、強制振動周波数を 0. 2に中心として f c = 0. 05、 0. 1、 0. 1 5、 0. 2、 0. 25、 0. 3、 0. 4、 0. 5の 8ケースで計算を行なった。静止円柱を含めて各振動数での円柱に作用する流体力の抗カ係数 C_Dおよび揚力係数の時間履歴の変動波形をそれぞれ図 1 9 A〜Iに示す。図 1 9 A〜Iの抗カ係数と揚力係数の時間履歴波形図から、振動周波数 f c = 0. 5を除いて、振動する円柱に作用する抗力および揚力が静止円柱のものより大きいことが分かった。

揚力の変動波により求められた流れの渦放出状態を表わす S t r o u h a 1数 S tと強制振動周波数 f cの関係を図 20に示す。 S t r o u h a l数 0. 2付近の 0. 1 5< f c<0. 25における f cノ S t = 1の計算結果から、強制振動による渦放出周波数のロックイン現象がよく再現されていることが分かった。また、図 1 9 A〜Iからロックイン領域での抗カゃ揚力が増大することも捉えている。

本発明と 0 k amo t oらの実験 [非特許文献 1 8]とはレイノルズ数（R e = 1 9000) や円柱振動振幅（A = 0. 05) とも異なるが、〇 k amo t oらの実験による 0. 1 86< f c<0. 2 1 6のロックイン範囲は本発明のロックィン領域に含まれることが確認された。

4. 4 円柱の渦放出による自励振動

4. 3の解析条件で、ニューマークの] 3法により式（20) の振動方程式を積分し、円柱の渦自励振動の振動変位や振動速度を求める。円柱の振動速度をさらに流れ場の速度境界条件として課する。流れ場のレイノルズ数 R eを 1 0000、スクル一トン数 S cを 5、減衰係数 hを 1 %、無次元流速 Vr=Uo/(foD)を 2 , 3， 4， 5, 6， 7, 8の七ケースで計算を行い、各々の無次元流速に対する抗力係数 C_D、揚力係数 Cいそして変位 yの時間履歴を図 2 1 A〜Cと図 22D〜 Gに示す。

図 2 1 A〜Cと図 22D〜Gに示された円柱の渦励振解析について、流れ方向に対して直角方向に振動する弾性支持された剛体構造物の 1質点 1自由度系モデルを仮定しているが、図 2 1 A〜Cと図 22D〜Gの計算結果から分かるように、式（20) の振動方程式の右辺の揚力項は、非定常な渦放出により不規則的な変化を呈している。そして流体構造の連成運動としてさらに振動円柱が流れ場に非定常な擾乱をも与えているため、 4. 3の s i n強制振動に比べて、渦励振による揚力と抗力の平均値は、静止円柱のものより増大することが一致しているものの、その時刻歴変化曲線はもつと不規則的な波形となっていることが見られる。また無次元流速 V r = 4 (円柱の固有振動数 f 0= 1 4 = 0. 2 5) の場合に円柱の振幅がもっとも大きくなることが分かった。これは静止円柱の流れ場の S t r o u h a 1数 S t = 0. 25にも関連すると考えられる。

各々の無次元流速 V rに対する yの最大変位を整理して図 23に示す。最大変位 yについてみると、 V r = 2、 3、 4では yの最大変位は次第に大きくなり、 V r = 4で最大になる。その後、 V rが大きくなるにつれて、 yの最大.変位は波状に小さくなつていく。本計算での渦励振応答振幅の変化傾向は、近藤 [非特許文献 1 9]の ALE有限要素法による計算結果と定性的に一致していると確認された。 5. 結論

(1 ) V-CADデータを直接利用する流体解析方法を開発し、内部流，外部流ともにベンチマークテストを用いて検証した。その結果、本発明で提案した流体解析方法は十分な安定性を有し、二次元静止円柱まわりの流れの計算では、従来のポクセル法より良く、境界適合格子と比べて遜色がない計算精度を持つことが検証された。

(2) 流れ場の形状に関しては、どんな複雑な形状があっても、 V- CADデー夕さえを読み込めば、短時間で直接二次元の非圧縮性粘性流れの数値シミュレ一ションができるようになった。

(3) 二次元円柱の流体構造連成振動の計算では、円柱まわり流れのレイノルズ数 R e = 1 0 00の場合に、円柱の強制振動によるロックイン（ f cZS t = 1) 領域は 0. 1 5< f c<0. 25であることが予測された。

また、レイノルズ数 R e = 1 0 000 かつ無次元流速 V r = 4の時、つまり、円柱固有振動数（ f 0 = 1/4) が静止円柱の S t r o u h a l数（S t = 0. 2 5 ) と等しくなる時に、渦自励振動による共振状態が起こることが示された。これらの計算結果は、実験データや ALE有限要素法による結果と定性的に一致している。したがって、移動境界の流体構造連成問題においても、提案した流体解析方法の有効性が示されたと言える。

【発明の効果】

上述したように、次世代 CADとして期待されている V-CADポクセルデー夕を直接利用する任意形状の非圧縮粘性流れ場の数値解析方法の開発を行った。流れ場の境界には、 VOF法を併用したカットセル有限体積法を用いた。二次元のチャンネル内流れと円柱まわり流れを用いて調べた結果、この方法はある程度の精度と安定性を有し、計算コストがあまりかからない実用計算方法であることが分かった。 - 従って、本発明の V-CADデータを直接用いた非圧縮性粘性流体の流れ場の数値解析方法と装置は、格子生成を完全に自動化することができ、物体境界の表現が容易であり、比較的簡単な計算処理により、精度の高いシミュレーションを短時間で行うことができる等の優れた効果を有する。

Claims

請求の範囲

1. 非圧縮性粘性流体と接する対象物の境界データからなる外部データ（1 2) を境界が直交する複数のセル（1 3) に分割する分割ステップ（A) と、分割された各セルを対象物の内側に位置する内部セル（1 3 a) と境界デ一タを含む境界セル（1 3 b) とに区分するセル区分ステップ（B) と、

前記境界デ一夕による境界セル（1 3 b) の稜線の切断点を求める切断点決定ステップ（C) と、

求めた切断点を結ぶ多角形を境界面のセル内部データとする境界面決定ステツプ（D) と、

流れ場の境界に、 VO F法を併用した力ットセル有限体積法を適用して解析する解析ステップ（E) と、を備えたことを特徴とする V-CADデータを直接用いた非圧縮性粘性流体の流れ場の数値解析方法。

2. 前記解析ステップ（E) において、空間積分について、対流項に二次元の QU I.CK補間スキームを使用し、拡散項に二次精度の中心差分を用い、時間進行について、対流項と拡散項を合わせて二次精度の Ad ams-B a s h f o r t h法を適用し、圧力勾配項に一次精度の Eu 1 e r陰解法を用いる、ことを特徴とする請求項 1に記載の非圧縮性粘性流体の流れ場の数値解析方法。

3. 二次元境界セルにおいて、有限体積法における支配方程式を [数 1] の式（7) であらわす、

【数 1】

ことを特徴とする請求項 2に記載の非圧縮性粘性流体の流れ場の数値解析方法。

4. 有限体積法における支配方程式における対流項、圧力勾配項および拡散項を [数 2] の式（8) 、 (9) 、 (1 0) であらわす、

【数 2】対流項： ·

=[△ ( O厂 -,.

+Ax(A_t u_i ^{^ v_{I J} - _J-i« /2_J-./2V_i,_y-₁ )] ( ⁸ )

- ；： ,

圧力勾配項：

= Ay[B_{i PiW2} -B^_jPi_ _j -ρ_ρ{Β^. -Β^.)ΐ ( 9 )

+Ax[A_iJp_u+ -A_i __]p_i __U2-p_p(A_i -A,.^)]] 拡散項： '.

(10)

+Ax(A grad(v) _j+V2 - 4,ゾー, ( — - (4, - 4, -】)

ことを特徴とする請求項 3に記載の非圧縮性粘性流体の流れ場の数値解析方法。

5. 固体境界の積分において、体境界でノースリップ境界条件の場合に、対流項は零とし、圧力勾配項と拡散項に対しては、切断線の中間点 Pの値を平均値として用いて積分し、空間積分に対してはすべての項に開口率を適用する、ことを特徴とする請求項 3に記載の非圧縮性粘性流体の流れ場の数値解析方法。

6. カットセルにおいて算出される物理変数の定義点を、 VOF= 0. 0 1を閾値としてそれより小さい境界セルは完全な固体とみなし、それより大きい境界セルにおいて算出される変数をセルの中心に置く、また、稜線の変数定義点をセル稜線の中心に定義し、更に、線分 4の中心点の変数値を線形補間により求める、ことを特徴とする請求項 3に記載の非圧縮性粘性流体の流れ場の数値解析方法。

7. 物体に働く抗カ（流れ方向の力）と揚力（流れの垂直方向の力）を、 [数 3 ] の式（1 2) ( 1 3) であらわす、

【数 3】

抗カ：

Wi ( +

揚力：

_ΓΓ σび d(T

y dx dy

= if(^ ) hdv+ + i ffif(- v ^ = び ('び _wし、一びし ,、) し、一び,"

ことを特徴とする請求項 3に記載の非圧縮性粘性体の流れ場の数値解析方法。

8. 移動境界を伴う流体構造連成解析において、所定の時間刻み毎に、流体系と構造系を別々に解析し、それぞれの境界条件を陽的に受け渡す、ことを特徴とする請求項 3に記載の非圧縮性粘性流体の流れ場の数値解析方法。

9. 円柱の強制振動解析において、円柱が弹性支持された剛体構造物として、流れ方向に対して直角方向に振動する 1質点 1自由度系モデルと仮定し、円柱中心の Y方向の変位を式（1 7) で与え、円柱表面の Y方向の速度境界条件を、式（1 8) で与える、 y = Asin(2irf_ct) (17)

ことを特徴とする請求項 8に記載の非圧縮性粘性流体の流れ場の数値解析方法。

1 0. 式（1 8) から求められた振動する円柱の運動速度を、流れ場の速度境界条件として計算時間ステップごとに変化して与える、ことを特徴とする請求項 9に記載の非圧縮性粘性流体の流れ場の数値解析方法。

1 1. 円柱の渦放出による自励振動解析において、有次元の振動方程式を 1質点 1自由度のダンバ ·ばねモデルとて、式（1 9) 又は式（20) であらわす、

1 2. 式（20) から求められた振動する円柱の運動速度を、流れ場の速度境界条件として計算時間ステップごとに変化して与える、ことを特徴とする請求項 1 1に記載の非圧縮性粘性流体の流れ場の数値解析方法。

1 3. 円柱の初期変位と初期速度を 0とし、かつ式（20) の揚力は現在の値を用いて陽的に与え、式（20) の振動方程式をニューマークの) 3法により積分し，円柱の振動変位と振動速度を求める、ことを.特徴とする請求項 1 1に記載の非圧縮性粘性流体の流れ場の数値解析方法。

14. 非圧縮性粘性流体と接する対象物（1) の境界データからなる外部データ（1 2〉を入力する入力装置（2) と、形状と物理量を統合した実体データとその記憶演算プログラムを記憶する外部記憶装置（3) と、前記記憶プログラムを実行するための内部記憶装置（4) 及び中央処理装置（5) と、実行結果を出力する出力装置（6) とを備え、

前記外部デ一夕を境界が直交する複数のセル（ 1 3.) に分割し、分割された各セルを対象物の内側に位置する内部セル（1 3 a) と境界データを含む境界セル (1 3 b) とに区分し、前記境界データによる境界セル（1 3 b) の稜線の切断点を求め、求めた切断点を結ぶ多角形を境界面のセル内部デ一夕とし、流れ場の境界に、 VOF法を併用したカットセル有限体積法を適用して解析する、ことを特徴とする V-CADデータを直接用いた非圧縮性粘性流体の流れ場の数値解析装置。

1 5. コンピュータに、非圧縮性粘性流体と接する対象物の境界データからなる外部データ（1 2) を境界が直交する複数のセル（1 3) に分割する分割ステップ（A) と、

分割された各セルを対象物の内側に位置する内部セル（1 3 a) と境界データを含む境界セル（1 3 b) とに区分するセル区分ステップ（B) と、

前記境界データによる境界セル（1 3 b) の稜線の切断点を求める切断点決定ステップ（C) と、

流れ場の境界に、 VOF法を併用したカツトセル有限体積法を適用して解析する解析ステップ（E) と、を実行させるための非圧縮性粘性流体の流れ場の数値解析プログラム。