/9/5 FD( 計算流体力学 ) の基礎理論 性能 運動分野 夏の学校 神戸大学大学院海事科学研究科勝井辰博 流体の質量保存 流体要素内の質量の増加率 [ 単位時間当たりの増加量 ] 単位時間に流体要素に流入する質量 流体要素 Fl lm (orol olm) v ( ) ガウスの定理 v( )
/9/5 = =( ) b=b =(b b b ) b= b = b + b + b アインシュタイン表記 ( 総和規約 ) ) ( v ガウスの定理運動量保存流体要素内の質量の増加率 [ 単位時間当たりの増加量 ] 単位時間に流体要素に流入する運動量流体要素に働く力の総和 ( ) より外力を F として F ) ( 流体要素 Fl lm (orol olm)
/9/5 ガウスの発散定理及びアインシュタイン表記を適用すると ( ) ( ) ( ) ( ) v( ) v( ) v( ) テンソル 書き換えると F 外力 は大きく つに分類される 面積力: 圧力 粘性力 境界面を通して流体要素に作用 体積力: 重力 遠心力 コリオリ力 電磁力 流体要素内の各部分に作用 T() 応力テンソル 検査面 左図のように点 を通る微小部分 を通して の内側から外側に作用する面積力を T() とするとき T() を応力テンソルと呼ぶ 考える点が同じであっても対応する面が異なれば応力は異なる為 T() も異なる ベクトルとテンソル速度ベクトル(yz)= () 場所と時間によっている 応力テンソル T() 場所と法線ベクトルと時間によっている 作用と反作用応力テンソルは面の表方向 裏方向に同様に作用している T()=- T( )
/9/5 応力テンソルは法線ベクトル によっているが は無限に存在する そこで下図のようにデカルト座標の座表面に平行な つの微小面の法線ベクトル を考える T( ) T( ) T( )=τ =τ +τ +τ T( )=τ T( )=τ T( ) : 応力テンソルの成分 方向に垂直な面に作用する 方向の単位面積当たりの面積力 τ を用いて T() はどう表されるか? - T() B:Δ :Δ B:Δ B:Δ B:Δ B = + + (= ) Δ = Δ Δ = Δ Δ = Δ - - 四角形 B に作用する面積力による力の総和は T( ) Δ + T( ) Δ + T( ) Δ + T() Δ= T( )=τ Δ = Δ とすると (τ ) Δ+ T()Δ= T()=(τ ) T() は τ で表される 4
/9/5 5 運動量保存 Ⅱ 運動量保存は次式で表される F ) ( 面積力を F とし 図よりガウスの発散定理と対称テンソルの性質を利用して F F ゆえに運動量保存は = τ F または非圧縮性流体の運動量方程式非圧縮性の流体の支配方程式 非圧縮性の流体の場合 τ は次のように表される < クロネッカーのデルタ > ) ( : ) ( :
/9/5 6 よって支配方程式は次のように書き換えられる 体積力の影響を考慮して移流項拡散項体積力 N 方程式 Fb バーガーズ方程式と偏微分方程式の性質 次元の N 方程式において圧力勾配を無視すると 移流速度を c 拡散係数を ( c ) として上式に代入し これをバーガーズ方程式と呼ぶ c バーガーズ方程式 階の偏微分方程式はの 種類に分類出来る 双曲型 放物型 楕円型
/9/5 7 双曲型 c ()= ( c) と置き c=x とする ) ( ) ( c X X X c ) ( ) ( X X X c ゆえに ) ( ) ( c c c つまり ()= ( c) のとき 必ず上記の微分方程式を満足する () は () を 方向に c 平行移動したものとなる () は時間とともに c 方向に移動速度 c で移動する 放物型 上に凸 時間とともに減少下に凸 時間とともに増加拡散を表す 楕円型 c 平均化
/9/5 8 有限差分法連続関数の導関数を離散点の値を用いて近似する方法ある関数 () の離散点を ( ) = とする は Δ おきに は Δ おきに定義された離散点テイラー展開より次の 式が得られる!!!! B +B より 4 O O 中心差分 より 前進差分 B より 後退差分
/9/5 9 陽解法と陰解法放物型の方程式において拡散係数を とすると この式から両辺に有限差分を適用して =+Δ の値 r = の値以上より陽解法及び陰解法は次のように定義できる r r 陽解法陰解法 ) ( r r
/9/5 次元定常移流拡散方程式 φ: 解くべき対象 : 移流速度 ( ここでは定数 ) Γ: 拡散係数 N 方程式 φ=( 解くべき で移流する 非線形方程式 ) 時間と共に φ は変化しない ργ: 一定 = で φ=φ =L で φ=φ L のディリクレ問題 解析解 L L :( ペクレ数 ) 常微分方程式 L φ φ L < = φ > L
/9/5 有限体積法と離散化 φ の補間直線近似 ( 結果は中心差分に相当 ) の評価 : 中心差分 ρ φ を一定とすると 4 4 局所ペクレ数 : L より L L
/9/5.775( ).5.689( ).9( 4) 4 4 ( ) L 8 ( ) L 4 ( 4) L ( 8) L 4 局所ペクレ数 : 移流と拡散の比 ペクレ数 ( 大 ) 移流が支配的となる ( 移流 : 上流の値がそのまま下流に伝わる ) つまり下流の影響は少なく 下流の情報が悪影響 が大 も大 遠い下流の情報が取り込まれる 一般に局所ペクレ数 で解が振動する 上流の影響だけを取り入れれば良い
/9/5 φ =φ φ =φ (> >) : 風上化 L L 4 8 6 4 4 5 L L L L L 風上化打ち切り誤差の評価風上補間セル界面 の上流計算点の値で φ を近似 次導関数を求める為に ( 流れ方向に依存して ) 後退 or 前進差分を用いることと同義 風上差分スキーム ( Drc chm) D ) ( ) (
/9/5 4 ρ > ρ > の時 : 風上差分 (ρ の評価点は少し異なる ) D は決して振動解を伴わないが数値拡散を伴う 拡散項の評価値と同じオーダーの量が移流項の誤差に入る 数値拡散 Δ を小さくしないと誤差が大きくなる 一方拡散項は ではが最大誤差 だけ少ない 次風上 が打ち切り誤差となる H ρ > の時 φ の 周りの T.( テイラー展開 )
/9/5 運動量方程式の離散化 直交正方スタッガード格子の導入 -+ + ++ -+ + ++ -+ + ++ - + - + - + -- - +- -- - +- -- - +- Δ Δ 移流項 - - -+ + ++ -+ + ++ - -+ + ++ - + - + - - + - + -- - +- -- - +- -- - -- - +- 5
/9/5 6 次風上で ~ を補間 - + - -- - - 圧力項セルが圧力により 方向に押される力 - + - -- - -
/9/5 7 拡散項 J - - - + - -- - - + - 移流 圧力 拡散 = y
/9/5 8 移流項 - - -+ ++ - + +- - -- + - -- - +- + -+ + - -- - +- -+ + ++ ++ + ++ + + - - + + + 次風上で ~ を補間 ++ + + - - + + +
/9/5 9 圧力項 セルが圧力により 方向に押される力 ++ + + - - + + + 拡散項 J - - ++ + + - - + + +
/9/5 y 移流 圧力 拡散 = locy rr olg IML(m Imlc rr Lk qo) by kr rcor orrcor = + = + = + corrc vl () : b b () : b b : 予測値 ステップ前の値 が求められる = + = + = + 正しい値 : 質量保存則と運動量保存則を満たす
/9/5 b b b b より b b b b の影響を ( とりあえず ) 無視 -+ ++ - + +- - -- + - -- - +- + -+ + - -- - +- -+ + ++ ++ + - + - + - 質量保存則
/9/5 + - + - + - 質量保存則 + - + - + - 質量保存則 N N