MyFirstOpenFOAM

Size: px

Start display at page:

Download "MyFirstOpenFOAM"

せいてらわ
4 years ago
Views:

1 はじめての OpenFOAM

2 メニュー OpenFOAMというパッケージの特徴メッシュの生成定常乱流の計算境界条件近似解法数値解法計算結果の確認

3 Open FOAM というパッケージの特徴流体をシミュレートするプログラム群メインはプログラムのライブラリ実行可能なファイルはライブラリをリンクしたものオープンソース必要な機能を持つライブラリを組み合わせて独自のシミュレートプログラムを作成できる Click! アプリではないリンク (Link) ソルバー (Solver) ライブラリ (Library)

4 主な機能前処理メッシュの作成操作計算内容定常 / 非定常非圧縮性 / 圧縮性乱流の近似モデル ( 時間平均 / 空間平均 ) 熱の輸送多相の乱流処理燃焼のモデル化

5 主な機能計算内容 ( つづき ) ニュートン流体 / 非ニュートン流体粒子追跡法メッシュ ( 境界 ) 移動結果処理可視化 (ParaViewを利用) 空間平均値の算出時間平均値の算出

6 サンプルモデルの概要サンプルモデル簡易スロットルボディ ( 内部流 + 乱流 ) 定常乱流場流速 100m/sec ( 非圧縮場を想定 ) 流体 20 の空気

7 サンプルモデルの境界条件サンプルモデル ~ 境界条件 U: 0m/sec P: zero 勾配 U: 0m/sec P: zero 勾配

8 計算スキームの設定 fvschemes とは基礎式を近似する数値解法を指定する例えば Navier-Stokes 方程式では以下の関係にある t p x 2 ( ρu) + div( ρuu) = X + µ u 非定常項 (ddtschemes) 対流 ( 移流 ) 項 (divschemes) 外力項圧力勾配項 (gradschemes) 粘性項 (laplacianschemes)

9 計算スキームの設定 fvschemes とは基礎式を近似する数値解法を指定する例えばレイノルズ方程式では以下の関係にある 2 t p x ( ρu ) + div ( ρuu) = X + µ u + div( τ ) R 非定常項 (ddtschemes) 対流 ( 移流 ) 項 (divschemes) 外力項圧力勾配項 (gradschemes) 粘性項 (laplacianschemes) レイノルズ応力項 (divschemes)

10 計算スキームの設定離散化スキームの選択風上差分スキーム upwind 上流側のセルから外挿補間する手法 1 次近似式中心差分スキーム linear TVD スキーム limited Linear 上流下流双方のセルを線形補間し自セルの値を内挿する手法 2 次近似式風上差分スキームと中心差分スキームの結果を合わせて自セルの値を決める重み付平均が可能

11 計算スキームの設定 fvschemes の内容 (1) FoamFile version 2.0; format ascii; class dictionary; object fvschemes; // * * * * * * * * * * * * * * ddtschemes default steadystate; gradschemes default Gauss linear; grad(u) celllimited Gauss linear 1; divschemes default none; div(phi,u) bounded Gauss linearupwindv grad(u); div(phi,k) bounded Gauss upwind; div(phi,omega) bounded Gauss upwind; div((nueff*dev(t(grad(u))))) Gauss linear;

12 計算スキームの設定 fvschemes の内容 (2) laplacianschemes default Gauss linear corrected; interpolationschemes default linear; sngradschemes default corrected; fluxrequired default no; p;

13 計算スキームの設定 fvschemes の内容 (1) FoamFile version 2.0; format ascii; class dictionary; object fvschemes; // * * * * * * * * * * * * * * ddtschemes default steadystate; gradschemes default Gauss linear; grad(u) celllimited Gauss linear 1; 2 t t Reynolds 方程式の例 p x ( ρu ) + div ( ρuu) = X + µ u + div( τ ) 非定常項 (ddtschemes) 乱流の運動エネルギーの例圧力勾配項 (gradschemes) µ k x τ k x u x t 2 ( ρk ) + div( ρuk) = + µ u ε 非定常項 (ddtschemes) t 速度勾配項 (gradschemes) R

14 計算スキームの設定 t fvschemes の内容 (1) 2 t Reynolds 方程式の例 p x ( ρu ) + div ( ρuu) = X + µ u + div( τ ) 速度の移流項 (divschemes) 乱流の運動エネルギーの例 µ k x τ k x u x t 2 ( ρk ) + div( ρuk) = + µ u ε 乱流強度の移流項 (divschemes) t R divschemes default none; div(phi,u) bounded Gauss linearupwindv grad(u); div(phi,k) bounded Gauss upwind; div(phi,omega) bounded Gauss upwind; div((nueff*dev(t(grad(u))))) Gauss linear;

15 計算スキームの設定 fvschemes の内容 (2) laplacianschemes default Gauss linear corrected; t Reynolds 方程式の例 2 p x ( ρu ) + div ( ρuu) = X + µ u + div( τ ) R interpolationschemes default linear;

16 計算スキームの設定 fvschemes の内容 (2) sngradschemes セルの直行性補正手法の選択 sngradschemes default corrected; fluxrequired default no; p;

17 数値解法の設定 fvsolution とは連立方程式の数値解法を指定 (linear solver) u u u u u u f e d c b a f e d c b a f e d c b a f e d c b a f e d c b a f e d c b a u 1 u 4 u 2 u 5 u 3 u 6 先のスキームを使った方程式がセル数分作られるこの連立方程式を解く数値解法を指定する連立方程式は場ごとに個別に存在する (u,p,k,ε, )

18 連立方程式の解法の設定 fvsolution の内容 (1) FoamFile version 2.0; format ascii; class dictionary; object fvsolution; // * * * * * * * * * * * * * * * * solvers p solver GAMG; tolerance 1e-7; reltol 0.01; smoother GaussSeidel; npresweeps 0; npostsweeps 2; cacheagglomeration on; agglomerator faceareapair; U k ncellsincoarsestlevel 10; mergelevels 1; solver smoothsolver; smoother GaussSeidel; tolerance 1e-8; reltol 0.1; nsweeps 1; solver smoothsolver; smoother GaussSeidel; tolerance 1e-8; reltol 0.1; nsweeps 1;

19 連立方程式の解法の設定 fvsolution の内容 (2) omega solver smoothsolver; smoother GaussSeidel; tolerance 1e-8; reltol 0.1; nsweeps 1; SIMPLE nnonorthogonalcorrectors 0; potentialflow nnonorthogonalcorrectors 10; relaxationfactors fields p 0.3; equations U 0.7; k 0.7; omega 0.7; cache grad(u);

20 連立方程式の解法の設定 fvsolution の内容 (1) FoamFile version 2.0; format ascii; class dictionary; object fvsolution; // * * * * * * * * * * * * * * * * solvers p solver GAMG; tolerance 1e-7; reltol 0.01; smoother GaussSeidel; npresweeps 0; npostsweeps 2; cacheagglomeration on; agglomerator faceareapair; 圧力場の連立方程式に関する指定解法 :GAMG ( 代数的多重格子法 ) 予め荒い格子で概算を算出しその結果を初期値に本来の計算を行う手法 1 スイープの許容残差 1 ステップ中の許容残差発散する場合下げると改善する場合がある

21 連立方程式の解法の設定 fvsolution の内容 (1) U k ncellsincoarsestlevel 10; mergelevels 1; solver smoothsolver; smoother GaussSeidel; tolerance 1e-8; reltol 0.1; nsweeps 1; solver smoothsolver; smoother GaussSeidel; tolerance 1e-8; reltol 0.1; nsweeps 1;

22 連立方程式の解法の設定 fvsolution の内容 (2) omega solver smoothsolver; smoother GaussSeidel; tolerance 1e-8; reltol 0.1; nsweeps 1; SIMPLE nnonorthogonalcorrectors 0; potentialflow nnonorthogonalcorrectors 10; メッシュの直行性を補正する処理の繰り返し回数低品質メッシュに有効だが計算時間が増加する

23 連立方程式の解法の設定 fvsolution の内容 (2) 緩和係数次ステップに渡す補正量を低減させ急激な値の変化による解の振動を抑える SIMPLE 法は特に緩和係数が解の安定化に大きな役割を果たしている relaxationfactors fields p 0.3; equations U 0.7; k 0.7; omega 0.7; cache grad(u);

24 ベースメッシュの作成 blockmeshdict の内容 FoamFile version 2.0; format ascii; class dictionary; location system; object blockmeshdict; // * * * * * * * * * * * * * * converttometers 1; vertices( ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) ); blocks ( hex ( ) ( ) simplegrading (1 1 1) ); edges(); Patches ( patch BOXInlet ( ( ) ) patch BOXOutlet ( ( ) ) wall BOXRightWall ( ( ) ) wall BOXLeftWall ( ( ) ) wall BOXTopWall ( ( ) ) wall BOXGround ( ( ) ) ); mergepatchpairs();

25 ベースメッシュの作成 blockmeshdict の内容 FoamFile version 2.0; format ascii; class dictionary; location system; object blockmeshdict; // * * * * * * * * * * * * * * converttometers 1; vertices( ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) ); ファイルヘッダ dictionary blockmeshdict という箇所だけ注意単位変換の係数この係数倍した結果が m になるよう指定形状を決める頂点の座標上位順に頂点 0, 頂点 1, と番号が振られる

26 ベースメッシュの作成 blockmeshdict の内容メッシュ空間の外形と分割数間隔の指定 hex メッシュ空間外観が 6 面体であると定義 ( ) 頂点番号で 6 面体を定義順番は右手系 ( ) 空間の分割数順に X Y Z 方向の分割数 simplegrading 単純な倍率でメッシュ間隔を定義する (1 1 1) X,Y,Z 方向のセルの間隔倍率 1 で等距離境界条件の種類と固有名定義面の位置 patch: 何かしら値が定義される面 wall: 壁面境界面は囲む頂点で指定同時に複数個指定が可能 blocks ( hex ( ) ( ) simplegrading (1 1 1) ); edges(); Patches ( patch BOXInlet ( ( ) ) patch BOXOutlet ( ( ) ) wall BOXRightWall ( ( ) ) wall BOXLeftWall ( ( ) ) wall BOXTopWall ( ( ) ) wall BOXGround ( ( ) ) ); mergepatchpairs();

27 ベースメッシュの作成 blockmesh の実行次のようにファイルを配置するケースディレクトリ constant polymesh blockmeshdict system controldict コンソールから起動する blockmeshdict を polymesh の中に配置した場合 $> blockmesh blockmeshdict を他の場所に ( 例えば system の中 ) 配置した場合 $> blockmesh dict system/blockmeshdict

28 乱流境界層の近似 Log-law の適用とメッシュサイズ U 0.99U 境界層厚さ壁面近傍速度は Log の近似関数 ( 壁関数 ) で近似される

29 乱流境界層の近似 Log-law の適用とメッシュサイズセル中心境界層厚さ壁面境界条件が定義されたセルの中心が境界層厚さに等しいとき Log-law による壁面近傍近似が正確に算出される

30 乱流境界層の近似 U+ と y+ 壁面せん断応力 τ 複雑な壁面形状では τ は一意に求めることが出来ませんそこで平板表面を沿う流れ場の τ で近似します壁面せん断応力 τ ρ 2=0.0359Re. τ =ρ Re. y 摩擦速度 = τ ρ = Re. + = + = ν + = Re. ν = + ν Re.

31 乱流境界層の近似 Log-law の適用とメッシュサイズ壁面せん断応力 τ U + 摩擦速度 = τ ρ + = y + + = ν 壁関数の適用可能範囲ターゲット y + は 100 ( 対数軸 )

32 計算メッシュの作成サンプルモデルベースメッシュの準備 CAD に合わせてメッシュを抜き出しメッシュの詳細化ベースメッシュから CAD 形状にあったメッシュを抜き出して必要に応じて詳細メッシュ化サーフェス上へのレイヤーセルの追加を行うレイヤーセルの追加

33 計算メッシュの作成 snappyhexmeshdict の内容 (1) FoamFile version 2.0; format ascii; class dictionary; location system; object snappyhexmeshdict; // * * * * * * * * * * * * * * * // running steps castellatedmesh true; snap true; addlayers true; Geometry Throttle.stl type trisurfacemesh; name Thrittle; pipewall.stl type trisurfacemesh; name pipewall; Inlet.stl type trisurfacemesh; name Inlet; Outlet.stl type trisurfacemesh; name Outlet;

34 計算メッシュの作成 snappyhexmeshdict の内容 (2) densitybox0 type searchablebox; min ( ); max ( ); castellatedmeshcontrols maxlocalcells ; maxglobalcells ; maxloadunbalance 0.10; minrefinementcells 0; ncellsbetweenlevels 3; resolvefeatureangle 30; features ( );

35 計算メッシュの作成 snappyhexmeshdict の内容 (3) refinementsurfaces Throttle level (9 9); pipewall level (9 9); Inlet level (9 9); Outlet level (9 9); refinementregions densitybox0 mode inside; levels ((1E15 6)); locationinmesh ( ); allowfreestandingzonefaces true; snapcontrols nsmoothpatch 3; tolerance 1.0; nsolveiter 300; nrelaxiter 5; nfeaturesnapiter 10; implicitfeaturesnap explicitfeaturesnap multiresionfeaturesnap true; false; true;

36 計算メッシュの作成 snappyhexmeshdict の内容 (4) addlayerscontrols relativesizes false; layers Throttle_surf nsurfacelayers 1; pipewall_surf nsurfacelayers 1; finallayerthickness ; expansionratio 1.0; minthickness ; ngrow 0; featureangle 90; nsmoothsurfacenormals 1; nsmoothnormals 3; nsmooththickness 10; maxfacethicknessratio 0.5; maxthicknesstomedialratio 0.3; // minmedianaxisangle 80; nbuffercellsnoextrude 0; nlayeriter 50; meshqualitycontrols::relaxed.; nrelaxediter 20; nrelaxiter 5;

37 計算メッシュの作成 snappyhexmeshdict の内容 (5) meshqualitycontrols maxnonortho 65; maxboundaryskewness 10; maxinternalskewness 4; maxconcave 80; minvol 1E-14; mintetquality 1e-11; minflatness 0.5; minarea ; mintwist 0.02; mindeterminant 0.01; minfaceweight 0.02; minvolratio 0.01; mintriangletwist -1; //minvolcollapseratio 0.01; nsmoothscale 4; errorreduction 0.75; relaxed maxnonortho 65; debug 0; mergetolerance 1E-10; //

38 計算メッシュの作成 snappyhexmeshdict の内容 (1) FoamFile version 2.0; format ascii; class dictionary; location system; object snappyhexmeshdict; // * * * * * * * * * * * * * * * // running steps castellatedmesh true; snap true; addlayers true;

39 計算メッシュの作成 snappyhexmeshdict の内容 (1) CAD 等で作った形状データをメッシャに渡す STL ファイルのデータに名前を関連付ける Geometry Throttle.stl type trisurfacemesh; name Thrittle; pipewall.stl type trisurfacemesh; name pipewall; Inlet.stl type trisurfacemesh; name Inlet; Outlet.stl type trisurfacemesh; name Outlet;

40 計算メッシュの作成 snappyhexmeshdict の内容 (2) densitybox0 type searchablebox; min ( ); max ( ); 簡単な幾何学形状であれば STL ファイルが無くても定義できるここでは searchablebox キーワードで矩形柱定義している

41 計算メッシュの作成 snappyhexmeshdict の内容 (2) 並列処理 (MPI) でメッシュを作成することを前提とした諸設定 1 コア当たりのメッシュ生成数の上限全体のメッシュ生成数の上限コア間で分担するメッシュ数の不釣り合いが発生したらメッシュ生成中に各コアに再分配するこの係数は再分配を始める上限この例では 10% 以上不釣り合いになると再分配される castellatedmeshcontrols maxlocalcells ; maxglobalcells ; maxloadunbalance 0.10; minrefinementcells 0; ncellsbetweenlevels 3; resolvefeatureangle 30; features ( );

42 計算メッシュの作成 snappyhexmeshdict の内容 (3) refinementsurfaces Throttle level (9 9); pipewall level (9 9); Inlet level (9 9); Outlet level (9 9); refinementregions densitybox0 mode inside; levels ((1E15 6)); locationinmesh ( ); allowfreestandingzonefaces true; ベースメッシュを分割する回数形状に与えた名前ごとに定義できる level ([ 通常の分割回数 ] [ 角部の分割回数 ] ) で定義される指定した形状 ( 名 ) の内側に対しメッシュの分割回数を一律で定義する計算空間の中の任意の一点を指定する閉空間で区切られている場合他方の空間にあるメッシュは計算空間外と認識され削除される

43 計算メッシュの作成 snappyhexmeshdict の内容 (3) メッシュの面を形状面に沿うように変形させる処理に用いられる定義 nsolveriter を増やすと形状再現性が改善する [implicitfeaturesnap] 形状のエッジを探しながらメッシュの面を変形させる機能 [explicitfeaturesnap] 予め *.emesh 形式のファイルでエッジ線を準備しておきこれにメッシュのエッジを一致させる機能事前に不要なエッジを削除することができる snapcontrols nsmoothpatch 3; tolerance 1.0; nsolveiter 300; nrelaxiter 5; nfeaturesnapiter 10; implicitfeaturesnap explicitfeaturesnap multiresionfeaturesnap true; false; true;

44 計算メッシュの作成 snappyhexmeshdict の内容 (4) addlayerscontrols relativesizes false; layers Throttle_surf nsurfacelayers 1; pipewall_surf nsurfacelayers 1; finallayerthickness ; expansionratio 1.0; minthickness ; ngrow 0; featureangle 90; nsmoothsurfacenormals 1; nsmoothnormals 3; nsmooththickness 10; maxfacethicknessratio 0.5; maxthicknesstomedialratio 0.3; // 壁面近傍の境界層を精度よく再現するために用いるレイヤーセルの追加作業に用いられる設定壁面にあるメッシュのサイズに対する比率で指定するか直接メートル単位で指定するかを選択するこの例では直接メートル単位で指定しているレイヤーセルを配置する箇所とセルの層数を指定する場所の指定は [ 形状名 ]_[ パート名 ] で指定するパート名は STL ファイルであれば solid 名と一致している壁面から一番離れたレイヤーセル一層の厚さ壁面に接したセルのことではないので注意空間が狭い場合にどの距離までならレイヤーを追加するかという距離を指定これ以下であればレイヤーが削除される

45 計算メッシュの作成 snappyhexmeshdict の内容 (4) レイヤーセルの追加に伴い歪んだセルの形状を緩和する処理の繰り返し回数増やせばある程度品質が改善する minmedianaxisangle 80; nbuffercellsnoextrude 0; nlayeriter 50; meshqualitycontrols::relaxed.; nrelaxediter 20; nrelaxiter 5;

46 計算メッシュの作成 snappyhexmeshdict の内容 (5) meshqualitycontrols maxnonortho 65; maxboundaryskewness 10; maxinternalskewness 4; maxconcave 80; minvol 1E-14; mintetquality 1e-11; minflatness 0.5; minarea ; mintwist 0.02; mindeterminant 0.01; minfaceweight 0.02; minvolratio 0.01; mintriangletwist -1; //minvolcollapseratio 0.01; nsmoothscale 4; errorreduction 0.75; セルの非直行性を許容する最大角これ以上歪むセルを形状修正の対象にするメッシュ体積の最低サイズ負の体積セルが生成されることを防ぐメッシュの面に対する最低サイズ

47 計算メッシュの作成 snappyhexmeshdict の内容 (5) relaxed maxnonortho 65; debug 0; mergetolerance 1E-10; //

48 計算メッシュの作成 snappyhexmesh の実行次のようにファイルを配置するケースディレクトリ constant polymesh blockmeshdict コンソールから起動する単体コア (CPU) で実行する場合 $> snappyhexmesh -overwrite 複数コアで実行する場合 (MPI) $>decomposepar -constant $>foamjob p snappyhexmesh overwrite $>reconstructparmesh -constant system controldict snappyhexmeshdict fvschemes fvsolution 並列処理でメッシュを作成する場合 decomposepardict hostfile

49 計算条件の設定境界条件境界条件と計算結果は同じファイルに保存されるケースディレクトリ 0 p: 圧力 U: 速度 k: 乱流エネルキ epsilon: 乱流散逸率 constant polymesh blockmeshdict system controldict snappyhexmeshdict fvschemes fvsolution 並列処理でメッシュを作成する場合 decomposepardict

50 サンプルモデルの境界条件サンプルモデル ~ 境界条件 U: 0m/sec P: zero 勾配 U: 0m/sec P: zero 勾配

51 計算条件の設定境界条件乱流のパラメータの決め方乱流エネルギー [k] 3 k = UI 2 ( ) 2 ここで U は流入流速 I は乱流強度乱流強度のめやす : 高速圧縮性内部流などのケース:5~20% 非圧縮性 ~ 弱圧縮性などのケース :1~5% 外部流などのケース :1% 程度乱流散逸率 [ε] ε = C µ k l 3 2 ここで Cµ は乱流モデル定数 l は乱流長さスケール乱流長さスケールのめやす : 通常は水力直径 : 等価面積の円の直径の 3.8% 程度流入面の境界層厚さの 22% という方法もある

52 計算条件の設定境界条件 U U: 0m/sec FoamFile version 2.0; format ascii; class volvectorfield; object U; dimensions [ ]; internalfield uniform ( ); boundaryfield Inlet_surf type fixedvalue; value $internalfield; Outlet_surf type zerogradient; Throttle_surf type fixedvalue; value uniform (0 0 0); pipewall_surf type fixedvalue; value uniform (0 0 0);

53 計算条件の設定境界条件 p p:zero 勾配 FoamFile version 2.0; format ascii; class volscalarfield; object p; dimensions [ ]; internalfield uniform 0; boundaryfield Inlet_surf type zerogradient; Outlet_surf type fixedvalue; value uniform 0; Throttle_surf type zerogradient; pipewall_surf type zerogradient;

54 計算条件の設定境界条件 k k: 壁関数乱流強度のめやす : 高速圧縮性内部流などのケース:5~20% 非圧縮性 ~ 弱圧縮性などのケース :1~5% 外部流などのケース :1% 程度 FoamFile version 2.0; format ascii; class volscalarfield; object k; dimensions [ ]; internalfield uniform 0.05; boundaryfield Inlet_surf type fixedvalue; value $internalfield; Outlet_surf type zerogradient; Throttle_surf type kqrwallfunction; value $internalfield; pipewall_surf type kqrwallfunction; value $internalfield;

計算条件の設定境界条件 ε ε: : 壁関数上式とモデル定数より ε = k C µ l 3 2 0.05 = 0.09 0.1 3 2 0.01 ε = C µ k l 3 2 乱流長さスケール l のめやす : 水力直径の 3.8% 程度 Cµ:0.09 ( 乱流モデル定数 ) FoamFile version 2.

55 計算条件の設定境界条件 ε ε: : 壁関数上式とモデル定数より ε = k C µ l = ε = C µ k l 3 2 乱流長さスケール l のめやす : 水力直径の 3.8% 程度 Cµ:0.09 ( 乱流モデル定数 ) FoamFile version 2.0; format ascii; class volscalarfield; object epsilon; dimensions [ ]; internalfield uniform 0.01; boundaryfield Inlet_surf type fixedvalue; value $internalfield; Outlet_surf type zerogradient; Throttle_surf type epsilonwallfunction; value $internalfield; pipewall_surf type epsilonwallfunction; value $internalfield;

56 計算の実行シミュレーション実行の準備次のようにファイルを配置するケースディレクトリ 0 p: 圧力 U: 速度 k: 乱流エネルキ epsilon: 乱流散逸率 constant polymesh blockmeshdict system controldict snappyhexmeshdict fvschemes fvsolution 並列処理でメッシュを作成する場合 decomposepardict hostfile

57 計算の実行計算プログラムの実行 snappyhexmesh SIMPLE 法の計算を実行する定常 / 乱流 / 非圧縮性流体という特徴単体コア (CPU) で計算 $> simplefoam > log.simplefoam 2>&1 & もしくは $> foamjob simplefoam 複数コアで計算 (MPI) < ケースディレクトリ >/system/ 内に hostfile を作成以下のコマンドを実行 $> decomposepar force && foamjob p simplefoam

58 計算の実行計算を走らせてみましょう

Microsoft PowerPoint OpenFOAMの使い方(柴田).ppt [互換モード]

Microsoft PowerPoint OpenFOAMの使い方(柴田).ppt [互換モード] OpenFOAM 勉強会 for beginner 2 期第 1 回 OpenFOAM の使い方柴田貴裕目標既に OpenFOAM はインストール済みさまざまなチュートリアルに取り組めるようにチュートリアルの実行の方法を scalartransportfoam を例に用いて見ていくインストールは OpenCFD の HP の方法に従えば比較的容易にできる ) OpenCFD 社の HP http://www.openfoam.com/