名称未設定

  → DNS データーの理解が重要     高精度解析スキームの開発が重要     汎用性を考慮すると有限要素法の利用価値の向上     連成問題の必要性     LES の実用面への適用     RANS の地位の低下 第２章 流体力学基礎 2.1 基礎数学（テンソル代数）  ○スカラー表記（スカラーとは 0 階テンソルに対応している） !   代表的なスカラー量： 質量、密度、濃度、温度、エネルギー、圧力等  ○ベクトル表記（ベクトルとは 1 階テンソルに対応している） u = u

(

₁ u₂ u₃

)

= u v w

(

)

= u_i ｛２次元表現  u = u

(

1 u2

)

= u v

(

)

= ui｝   代表的なベクトル量： 運動量、速度、渦度、各種フラックス等  ○マトリクス表記（マトリクスとは 2 階テンソルに対応している） T = !_{1 1} !_{1 2} !_{1 3} !_{2 1} !_{2 2} !_{2 3} !_{3 1} !_{2 3} !_{3 3} " # $ $ % & ' ' =!i j   代表的な 2 階テンソル量： 応力、レイノルズ応力等 !_{i j} = 1 0 0 0 1 0 0 0 1 " # $ $ % & ' '   クロネッカーのデルタ  ○高階テンソル表記 !_ijm= 1 "1 0 # $ % & % (i, j,m) = (1,2,3), (2, 3,1), (3,1,2) (i, j,m) = (1,3, 2), (2,1,3), (3,2,1) i = j or j = m or m = i _{交代テンソル}  ○積 内積    A ! B = AiBi = A1B1+ A2B2+ A3B3    （ベクトル ベクトル → スカラー） 外積    A ! B ="ijm AjBm= A

(

2B3# A3B2 A3B1# A1B3 A1B2 # A2B1

)

   （ベクトル ベクトル → ベクトル） テンソル積  A ! B = A_iB_j = A₁B₁ A₁B₂ A₁B₃ A₂B₁ A₂B₂ A₂B₃ A₃B₁ A₃B₂ A₃B₃ " # $ $ % & ' '

   （ベクトル ベクトル → 2 階テンソル）   （注）Einstein の和の規約によりテンソルの足（添え字）がそろったもの      は各成分の和を取る。  ○ナブラ（勾配演算子） ! = " "x₁ " "x₂ " "x₃ # $ % & ' ( = " "x " "y " "z # $ % & ' = " "x_i ｛２次元表現 ! = " "x₁ " "x₂ # $ % & ' ( = " "x " "y # $ % & ' = " "x_i｝  ○ラプラシアン ! = " 2 "x2 + "2 "y2 + "2 "z2 = "2 "x_j"x_j  ○ベクトル解析表現との対応 gradF = !F = "F "x_i = "F "x "F "y "F "z # $ % & '    0 階テンソル → 1 階テンソル（テンソル積） divG = ! " G = #Gi #x_i = #G₁ #x + #G₂ #y + #G₃ #z    1 階テンソル → 0 階テンソル（内積） rotG = ! " G =#_ijm$Gm $x_j = $G₃ $y % $G₂ $z $G₁ $z % $G₃ $x $G₂ $x % $G₁ $y & ' ( ) *    1 階テンソル → 1 階テンソル（外積） 2.2 流体力学の基礎式  基礎方程式は保存則から導出する必要がある。   質量保存則    → 密度方程式   運動量保存則   → 運動方程式   エネルギー保存則 → エネルギー方程式  フラックス形式   検査体積内での物理量の時間的変化を検査面からの物理量の注入と検査体積内  の物理量の発生で表現する。 !F

!t dxdydz = (Mx" " Mx+)dydz + (My" " My+)dzdx + (Mz" " Mz+)dxdy + Sdxdydz    時間変化  yz 面の出入   zx 面の出入   xy 面の出入  発生 !F !t = " M_x+" M_x" dx " M_y+ " M_y" dy " M_z+ " M_z" dz + S = ! "M_x "x ! "M_y "y ! "M_z "z + S

 最終的なフラックス形式の方程式 !F !t = " !M_j !x_j + S   F が 0 階テンソル（スカラー）のとき、M は 1 階テンソル（ベクトル）   F が 1 階テンソル（ベクトル）のとき、M は 2 階テンソル   S はつねに F と同階テンソル

x

y

z

dx

dy

dz

フラックス形式の概念図  ○質量保存則（圧縮性流れ） F = ! M_j =!u_{j（移流）} S = 0 →  !" !t = # !"u_j !x_j   非圧縮性流れ（ ! = const ）の連続方程式 →  !u_j !x_j = 0  ○運動量保存則（圧縮性流れ） F_i =!u_i

M_{i j}=!u_iu_j"#_{i j（移流と応力）} S_i= f_{i（外力）} →  !"u_i !t = # !"u_iu_j !x_j + !$_{i j} !x_j + fi  Newton の応力の法則   !i j= " p#i j+ 2µ si j" 1 3 $u_m $x_m#i j % & ' ( ) *    s_{i j}= 1 2 !u_i !x_j + !u_j !x_i " # $ % & '  歪テンソル →  !"u_i !t = # !"u_iu_j !x_j # !p !x_i + 2µ ! !x_j si j# 1 3 !u_m !x_m$i j % & ' ( ) * + f_i !u_i !t = " !u_iu_j !x_j " ! !x_i p #+ µ # !2u_i !x_j!x_j + f_i # {変数変換 p ! p ", # = µ ", fi ! f_i " _}   非圧縮性流れにおける運動方程式 →  !u_i !t = " !u_iu_j !x_j " !p !x_i +# !2u_i !x_j!x_j + fi  ○エネルギー保存則（圧縮性流れ） F = !e +1 2!umum_{（内部エネルギーと運動エネルギーの和）} M_j =!eu_j+1 2!umumuj"#j mum + qj_{（移流と応力による仕事と熱供給）} S = 0 →  !"e !t + ! !t 1 2"umum # $ % & = ' !"eu_j !x_j ' ! !x_j 1 2"umumuj # $ % & + !(_{j m}u_m !x_j ' !q_j !x_j  Fourier の法則   qj = !" #T #x_j →  !"e !t + ! !t 1 2"umum # $ % & = ' !"eu_j !x_j ' ! !x_j 1 2"umumuj # $ % & ' !pu_j !x_j        +2µ ! !x_j sj mum" 1 3 !u_n !x_nuj # $ % & ' ( +) !2T !x_j!x_j

 ○パッシブスカラー量(温度や濃度)保存則（非圧縮性流れ） F = T M_j = Tu_j+ q_j S = 0 →  !T !t = " !Tu_j !x_j " !q_j !x_j   非圧縮条件と Fourier の法則により →  !T !t = "uj !T !x_j +# !2T !x_j!x_j 2.3 非圧縮性 Navier-Stokes 方程式周辺  この講義では非圧縮性流れを中心に話を進める。支配方程式である Navier-Stokes 方程式は以下のように表される。  ○連続方程式 !u_j !x_j = 0  ○運動方程式 !u_i !t + !u_iu_j !x_j = " !p !x_i +# !2u_i !x_j!x_j    左辺第１項 非定常項    左辺第２項 移流項＝非線形項    右辺第１項 圧力勾配項    右辺第２項 粘性拡散項   非定常項以外が全て空間微分でまとめられるこの表現を保存型とよぶ。   ガウスの定理 dV V !!! " #F =_!!SdSF # n   より、境界面からの流入・流出がなければ、保存するという意味。   連続方程式を用いて運動方程式の移流項を書き換えると !u_i !t + uj !u_i !x_j = " !p !x_i +# !2u_i !x_j!x_j   この表現を非保存型とよぶ。   （注）数学的には保存型も非保存型も同一のものであるが、数値計算では明瞭な違     いが現れる。   ラグランジェ微分 D Dt = ! !t+ uj ! !x_j

   ガリレイ不変性を満足する微分   移流項に現れた速度勾配テンソル !u_i !x_j = si j"#i j    歪テンソル   s_{i j}= 1 2 !u_i !x_j + !u_j !x_i " # $ % & '    渦度テンソル  !i j = 1 2 "u_j "x_i # "u_i "x_j $ % & ' ( )   歪テンソルは対称テンソルで、伸縮運動を意味している。   渦度テンソルは反対称テンソルで、回転運動を意味している。   渦度ベクトルの定義 !_i ="_ijm #um #x_j = #w #y $ #v #z #u #z $ #w #x #v #x $ #u #y % & ' ( )    流体要素の自転角速度の 2 倍を意味する。   渦度ベクトルと渦度テンソルの関係式 !_i ="_ijm!_{j m} !_{i j} = 1 2"ijm!m   渦度ベクトルを用いた Navier-Stokes 方程式 !u_i !t + "ijm#jum = $ ! !x_i p + 1 2ujuj % & ' ( +) !2u_i !x_j!x_j  ○圧力方程式   運動方程式の発散(div)をとると !2p !x_i!x_i = " !u_j !x_i !u_i !x_j = si jsi j"#i j#i j    ポアソン方程式(!p = A )   圧力は速度の 2 次の積で表現される。（圧力項＝非線形項）  ○渦度方程式   運動方程式の回転(rot)をとると D!_i Dt =!j "u_i "x_j +# "2!_i "x_j"x_j   となり、圧力項は消える。   歪テンソルを用いる表現

D!_i Dt =!jsi j+" #2!_i #x_j#x_j   右辺第１項は歪による渦の伸縮作用を表している。   速度と渦度の関係 ! = rotur   類似性： 磁場の強さと定常電流の関係（Maxwell 方程式）       ie = rotH → Biot-Savart の法則      定常電流がわかれば磁場の強さがわかる。    H(x) = 1 4! d 3 " x # ie( " x ) $ (x % " x ) x % " x 3 + grad&(x) 任意調和関数 !"(x) = 0 Biot-Savart の法則    渦度がわかれば速度がわかる。   u(x, t) = 1 4! d 3 " x # r $ ( " x ,t) % (x & " x ) x & " x 3 + grad'(x, t) 任意調和関数 !"(x, t) = 0  ○その他の重要な物理量   運動エネルギー   1 2uiui 非粘性流れでは保存量   エンストロフィー  1 2!i!i 散逸構造と関連 非粘性流れでは保存量   ヘリシティー   ui!i らせん運動と関連        3 次元特有の量 非粘性流れでは保存量  ○レイノルズ数   運動方程式の移流項（非線形項）と粘性拡散項の比率 Re = !u_iu_j !x_j " ! 2_u i !x_j!x_j =UL "    Re が小さい場合（粘性拡散項が強い）、流れは層流状態になる。    Re が大きい場合（非線形項が強い）、流れは乱流状態になる。   非線形項は流れを不安定にする効果がある。   粘性拡散項は流れを平滑化する。安定化効果がある。

   Case1 凸型    !2f !x2 < 0 → 粘性拡散項は減少するように働く。    Case2 凹型    !2f !x2 > 0 → 粘性拡散項は増加するように働く。 Case1 Case2 拡散項の働き   Navier-Stokes 方程式を特徴づける無次元パラメーター    代表長さ、速度、粘性率でスケーリングを行うと、Navier-Stokes 方程式は相    似則が成立している。 第３章 乱流基礎 3.1 乱流はなぜ生じるのか？ 3.1.1 乱流発生の研究  Reynolds(1883)は円管内流れにおいて水の流速や管の太さを変えると、    色素が筋を引いて流れる状態（層流）    色素が途中から急に乱れ出し管の断面全体に広がって流れる状態（乱流） が現れ、この状態の選択がレイノルズ数と呼ばれる無次元パラメーターで決められて いることを発見した。特に、彼の研究ではレイノルズ数 2300 が臨界値であった。  Ekman(1911)は追試によりレイノルズ数が 40000 を越えても層流であることを見 いだし、臨界レイノルズ数は 50000 以上となった。  → 層流から乱流への遷移に関する研究    微小擾乱に対する流れの安定性に関する研究 3.1.2 線形安定性理論  Navier-Stokes 方程式に対して、層流解とそれに負荷される２乗量を無視した擾乱

という立場から導出される線形偏微分方程式を解析する理論。  平行平板間の流れを考える。一定の力がx 方向に働いている。 F_i = F!_i1  この外力下で定常性一様性を仮定した Navier-Stokes 方程式は 0 = F_i+! " 2_u i "y2 となる。この条件下での層流解は以下である。（平面ポワズイユ流） u_i= F 2! L 2 " y2

(

)

#_i1 この流れを  代表長さ  L  代表速度 Uc = FL2 2! _{（チャネル中央速度）} で無次元化した微小擾乱を考慮した Navier-Stokes 方程式 !u_i !t + !u_iu_j !x_j = " !p !x_i + 1 Re !2u_i !x_j!x_j + 2 Re#i1 !u_j !x_j = 0 無次元化後の平面ポワズイユ流解はu i= 1! y

(

2

)

"i1_{と書ける。}

x

z

y

L

-L

平行平板間流と擾乱 層流解をu i、 p で表し、微小擾乱を ! u i、 ! p として微小擾乱の支配方程式を導くと

! " u _i !t + !u _iu " _j !x_j + ! " u _iu _j !x_j + ! " u _iu " _j !x_j = # ! " p !x_i + 1 Re !2u " _i !x_j!x_j ! " u _j !x_j = 0 となる。左辺第 4 項は非線形項である。微小擾乱の 2 乗項を無視すると ! " u _i !t + !u _iu " _j !x_j + ! " u _iu _j !x_j = # ! " p !x_i + 1 Re !2u " _i !x_j!x_j ! " u _j !x_j = 0 という線形辺微分方程式が得られる。この方程式において時間に関してラプラス変 換 Lf (x, ! ) =_#0"dtf (x,t)e$!t f (x, t) = 1 2!i const#i$d" const+ i$ % Lf (x, " )e"t を施す。結果として、同次方程式 !L " u _i+#u iL " u j #x_j + #L " u _iu _j #x_j + #L " p #x_i $ 1 Re #2L " u _i #x_j#x_j = 0 が求まる。この方程式を固有値問題として解析すればよい。その固有値! の中で実 部最大の固有値!Mを見いだせば、無限時間極限では u_i(x, t) ! f_i(x)e"Mt と振る舞う。 Re(!_M) < 0 _安定 擾乱は時間とともに減衰して、層流解に漸近していく。 Re(!M) = 0 中立 線形化の段階で近似を利用しているため、現実的にはほとんど 意味がない。 Re(!_M) > 0 _不安定 擾乱は時間とともに成長して、層流解から離れていく。  空間（x,z 方向）を波動で展開（フーリエ変換）する。 L ! v x,y,z,"

(

)

= ei kxx+ i kzz_{# k} x, y, kz,"

(

)

連続式を用いて圧力を消去すると! の方程式 u ! i " k_x # $ % & ' ( d2 dy2 ! kx 2 ! k_z2 # $ % & ' ( ) ! du dy = 1 ik_xRe d2 dy2 ! kx 2 ! k_z2 # $ % & ' ( 2 )

が得られる。この方程式は Orr-Sommerfeld 方程式と呼ばれる。この方程式から、固 有値! は波数kx、kzとレイノルズ数Re に依存する関数であることわかる。  臨界レイノルズ数の評価ためには! = 0 でのレイノルズ数の最小値を求めればよ く、方程式の対称性などを考慮すると k_x = 1.02 k_z = 0 Re_cr = 5772 が導かれる。この波動は主流方向のみに依存する２次元波で Tollmien-Schlichting 波と呼ばれるものである。  → 線形安定性理論で平面ポワズイユ流における微小擾乱は臨界レイノルズ数以    上で Tollmien-Schlichting 波を発生し、乱れが増幅することが見いだせた。  線形安定性理論での問題点  ○円管内の流れ（ポワズイユ流）では臨界レイノルズ数は見いだされておらず、    である可能性が高いが、現実に円管内乱流は存在している。  ○ Tollmien-Schlichting 波を発生以降の乱流状態への道筋が不明である。 → 非線形安定性理論？？の必要性 0 0.2 0.4 0.6 0.8 1 1.2 0 50000 100000 k x Re UNSTABLE STABLE 平面ポワズイユ流の線形安定特性（スケッチ） 3.1.3 非線形性の安定性への寄与  ある周波数のモード（Tollmien-Schlichting 波）が線形安定性理論にしたがって成 長すると、非線形相互作用が働き様々なモードの擾乱が発生する。  先にふれた平行平板間の流れではz 方向に依存した波動擾乱（３次元波）が生じて 最終的に乱流になることが実験的に観測されている。

  イメージ（拡散項：線形性）    新しいモードは生じず減衰していく。 0 0.5 1 1.5 2 0 1 2 3 4 5 6 7 8 9 10 !1 !2  →  0 0.5 1 1.5 2 0 1 2 3 4 5 6 7 8 9 10 !1 !2   イメージ（移流項：非線形性）    新しいモードは生じて、複雑化していく。 0 0.5 1 1.5 2 0 1 2 3 4 5 6 7 8 9 10 !2 !2"!1 !1 !2+!1  →  0 0.5 1 1.5 2 0 1 2 3 4 5 6 7 8 9 10 !2 !1 !2"!1 !2+!1 !2"2!1 2!1 3.1.4 乱流遷移の進展  以上のような層流から乱流への遷移は既に様々な研究により、唯一の遷移過程でな いことが知られている。  → 乱流への道筋の解明はできていない。  カオス研究として解の分岐をベルナール対流やテイラー・クエット流を中心として 多くの研究者によって研究が進められているのが現状である。

3.2 乱流の性質  ○ 非定常性   時々刻々激しく不規則に変動する。 -2 -1.5 -1 -0.5 0 0.5 1 1.5 2 0 1 2 3 4 5 6 7 8 9 10 u t   0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 0 1 2 3 4 5 6 7 8 9 10 K t ２次元一様等方性乱流の 1 点速度変動と平均運動エネルギーの時間変化  ○ 空間的変動性   空間的に激しく不規則に変動する -5 0 5 10 15 20 0 0.5 1 1.5 2 2.5 !"# $%&' ()#    !"#$ !$#% $#$ $#% "#$ "#$ $#% $#$ !$#% !"#$ 同心円環内乱流の軸方向速度      同心円環内乱流の瞬間渦度場

 ○ 強混合性（強拡散性）   層流に比べて、乱流は強い混合効果がある。  ○ 強抵抗性   層流に比べて、乱流は流れづらい。小さな流量。 0 10 20 30 40 50 60 70 80 90 0 0.5 1 1.5 2 2.5 !"# $"%&' $"()' 同心円環内乱流の軸方向速度  ○ 平均レベルにおける自己保存性  壁乱流（平板境界層、チャネル流）や自由乱流（噴流、混合層）には特定の速度ス ケール、長さスケールで無次元化すると平均速度分布等が一致する性質。  例 壁法則   摩擦速度（壁面摩擦の平方根）と動粘度（分子粘性）で無次元化    壁面からの距離の無次元量 y + ₌u!y "    平均速度の無次元量 U + = U u_!   粘性底層 U+ = y+   対数層 U+ = 2.45log y++ 5.0

0 2 4 6 8 10 12 14 16 18 20 0.2 1 10 100 200 U+ y+ チャネル乱流の主流方向速度 3.3 乱流の統計性 3.3.1 乱流の現象論(Kolmogorov 理論)  乱流の種類に依存しない普遍的な性質をもつ小さいスケールの運動に着目  速度のフーリエ変換 u_i(x, t) =!!!d3keik"xu ˜ _i(k,t)  運動エネルギー（揺動速度であれば乱流エネルギー） K =!!!d3x1 2ui(x, t)ui(x, t) = dk! E k

( )

 エネルギースペクトル関数 E k

( )

= 2!k2u ˜ _i(k,t) ˜ u _i("k, t)   ある波数に対応したスケールの渦がどれくらいエネルギーを持っているかを   意味している。高波数は小規模スケール、低波数は大規模スケールに対応。 1x10-13 1x10-12 1x10-11 1x10-10 1x10-9 1x10-8 1x10-7 1x10-6 1x10-5 1x10-4 1x10-3 1x10-2 1x10-1 1x100 1 10 100 300 E(k) k ２次元一様等方乱流のエネルギースペクトル

 ○ Kolmogorov の第１仮説   乱流のエネルギースペクトルの小スケール成分は平均のエネルギー散逸率! と   分子粘性! のみに依存する。    Kolmogorov 波数   kK =!1/ 4"#3/4    Kolmogorov 長さスケール    lK = 2!"#1 / 4$3 / 4    Kolmogorov 相似則    E k

( )

=!5/4"1/4f k / k

(

_K

)

    この相似則が成り立つ領域を普遍領域とよぶ。    → 平均散逸率と分子粘性との独立性を仮定  ○ Kolmogorov の第２仮説   レイノルズ数が極めて大きい乱流では普遍領域の低波数側に分子粘性には依存   しない領域が存在する。   Kolmogorov スペクトル（5/3 乗則） E k

( )

= C_K!2 / 3k"5/3 C_{K: Kolmogorov 定数 約 1.5}

E( k)

k

! k

"5/ 3 !"#$% &'() *+,() -.() 3.3.2 乱流場の正規性  乱流の揺動速度場の統計的性質は正規分布（ガウス分布）にほぼ従う。 Pdf (u) = 1 2!" exp # u2 2"2 $ % & ' ( )

 スキューネス < uuu >= 0   フラットネス < uuuu >= 3 < uu >2 3.3.3 間欠性理論  Kolmogorov 理論からの相違、乱流場の非正規性、非一様な散逸、フラクタル性 3.3.4 乱流構造  渦構造のチューブやシートタイプ、壁乱流における低速ストリーク構造、馬蹄渦 ２次元一様等方性乱流の瞬間渦度場 チャネル乱流のストリーク構造と渦構造

第４章 乱流計算 4.1 直接数値計算：DNS  ○流れの支配方程式（Navier-Stokes 方程式）をモデルを一切導入せずに、完全に計  算する。 !ui !t + !uiuj !xj = "!p !xi +# ! 2 ui !xj!xj !uj !xj = 0  ○散逸領域の渦まで解像する必要がある。   （乱流場のエネルギーが分子粘性により熱に変換される領域まで解像する    必要がある。）    → 乱流場のレイノルズ数に対する最低限の格子点数が決まる   エネルギー供給領域の渦のスケール（積分長） !I = K3/2 !   Kolmogorov スケール（粘性散逸領域の渦のスケール） !K = !3/ 4 "#1/ 4   両スケールを解像するのに必要な計算格子点数（１次元） N =!K !I = K 3/2 !3/ 4"3/ 4 = Re 3/ 4    ここで Re は乱流レイノルズ数である。   Re = K2 !"   両スケールを解像するのに必要な計算格子点数（3 次元） N3 = Re9 / 4    → 膨大な格子点数を用いた解析が必要。  ○３次元非定常解析    → 長い計算時間がかかる。  ○数値解を得ることができる。 → 実験に匹敵および凌駕   平均データや瞬間揺動場など全て得られる。  ○精度の高い計算スキームが必要。   擬スペクトル法が主流。高次中心差分法。   いかなる人工粘性も許可されない。→ 最良の計算方法が必須  ○航空関連分野では数値粘性などを用いたモデルを用いない解析を DNS と呼ぶこ   とがある。

E( k)

k ! k"5/ 3 !"#$% &'() *+,() -.() /0() DNS の計算領域 4.2 ラージ・エディ・シミュレーション：LES  ○流れの支配方程式（Navier-Stokes 方程式）にフィルター関数をかけて、以下のよ   うにグリッド成分とサブグリッド成分に分離する。 u_iGS

(

x, t

)

="""d3x ! G x, !

(

x

)

u_i

(

x , t!

)

u_iSGS

(

x, t

)

= u_i

(

x, t

)

! u_iGS

(

x, t

)

  フィルター関数としては    ガウシアンフィルター   G x

( )

= 6 !"2 exp # 6x2 "2 $ % & ' ( )    スペクトルカットオフフィルター   G x

( )

= 2 sin !x " # $ % & ' ( !x    トップハットフィルター    G x

( )

= 1 ! x " ! 2 0 x # ! 2 $ % & ' &   などが用いられている。ここで、Δ はフィルター代表長さである。   運動方程式 !u_iGS !t + !u_iGSuGS_j !x_j = " !pGS !x_i +# !2u_iGS !x_j!x_j " !$GS !x_j   連続の式

!u_jGS !x_j = 0    フィルター相関項  !GS = uiuj GS " u_iGSu_jGS = L_{i j}+ C_{i j}+ R_{i j}    Leonard 項   Li j= uiGSujGS GS ! uiGSuGSj    Cross 項  Ci j= uiGSujSGS GS+ uiSGSuGS GSj    SGS レイノルズ応力項   Ri j = uiSGSujSGS GS   フィルター相関各項は未知変数であるためモデル化が必要である。    → SGS モデリング  ○大きなスケールの運動（GS 成分）に関しては直接数値解析し、小さなスケールの   運動（SGS 成分）にのみモデルを用いて解析する手法  ○散逸領域にフィルターをかけると DNS と一致する。  ○普遍性が成り立っている慣性領域にフィルターをかけるべき。  ○３次元非定常解析    → それなりに計算時間がかかる。  ○平均データや瞬間揺動場などが得られる。  ○精度の高い計算スキームが必要。中心差分法。   高次の人工粘性も一般的には導入するべきではない。

E( k)

k ! k"5/ 3 !"#$% &'() *+,() -.() /0()

FILTER

LES の計算領域

4.2.1 スマゴリンスキーモデル  ０方程式型 SGS モデル L_{i j}= 0 C_{i j}= 0 R_{i j} =2 3K SGS !_{i j}"#SGS $ui GS $x_j + $u_jGS $x_i % & ' ( ) * !SGS = C

(

_S"

)

2 1 2 #u_iGS #x_j + #u_jGS #x_i $ % & ' ( ) #ui GS #x_j + #uGS_j #x_i $ % & ' ( )  ○高い安定性、計算実績   モデル定数の最適化か必要、壁減衰関数が必要 4.2.2 1 方程式型 SGS モデル  Schumann により発案されたモデル。

 Yoshizawa, Horiuti, Hamba, Okamoto により統計理論から構築されたモデル。

L_{i j}= 0 C_{i j}= 0 R_{i j} =2 3K SGS !_{i j}"#SGS $ui GS $x_j + $u_jGS $x_i % & ' ( ) * !SGS = C_!"KSGS1/ 2 DKSGS Dt = ! 1 2Ri j "u_iGS "x_j + "u_jGS "x_i # $ % & ' ( !)SGS+ DSGS !SGS= C_! K SGS3 / 2 " + 2# $ Ks $x_j $ Ks $x_j DSGS= ! !x_j Cd"K SGS1/2 +#

(

)

!KSGS !x_j $ % & ' ( )  ○高い安定性   モデル定数の最適化か必要、壁減衰関数が必要   → 近年ダイナミックタイプとして森西や梶島らにより再構築されている。 この系統のモデルにはこの他にも Bardina モデル、Clark モデル、非平衡 0 方程式 SGS モデル、非線形 SGS モデル、代数応力 SGS モデル、高次 1 方程式型 SGS モデル などがある。 4.2.3 ダイナミックモデル（ダイナミックスマゴリンスキーモデル）  Lilly の提案による新たな LES 計算法  フィルター操作を２回行う。

 １重フィルター操作の SGS 応力 !_{i j}= u_iu_j " u _iu _j  ２重フィルター操作の SGS 応力 Tij = uiuj ! !u"i " uj  １重フィルター操作の SGS 応力にフィルターをもう一度かけ、差をとった応力 lij = Tij! ! "ij= uiuj " !u!i ! u ! uiuj " + u" = uiuj " !iuj ! ui ! u  これは全て既知量でかけている。  仮にスマゴリンスキーモデルを仮定 Tij = C! 2 ! S S!ij ! !ij= C" 2 S Sij "  自動的にモデル定数が C = lijSij !2 S! S!ijSij" ! 2 S Sij "_S ij で決まる。  ○モデル定数の最適化か不必要、壁減衰関数が不必要、高性能   → 現在の LES はほぼダイナミックモデルで統一され、LES はモデリングレベル     の研究は理論的背景を有するならともかく、完成されているといえる。 4.3 アンサンブル平均モデル計算：RANS  ○流れの支配方程式（Navier-Stokes 方程式）を平均化する。   平均場と揺動場の分解 f = F + ! f f = F   平均操作の性質 F = F ! f = 0  平均操作 !u_i !t + !u_iu_j !x_j = " !p !x_i +# !2u_i !x_j!x_j ! u_i !t + ! u_iu_j !x_j = " ! p !x_i +# !2 u_i !x_j!x_j

!U_i !t + ! u_iu_j !x_j = " !P !x_i +# !2U_i !x_j!x_j  ２体相関項 u_iu_j = U_iU_j + U_iu ! _j + u _i! U_j + u _i! u ! _j    = UiUj + Ui u ! j + u i! Uj+ u i! u ! j = UiUj + u ! iu ! j  平均場方程式 !U_i !t + !U_iU_j !x_j + ! " u _iu " _j !x_j = # !P !x_i +$ !2U_i !x_j!x_j !U_j !x_j = 0  RANS ではこの方程式を解析する。   ! u iu ! j :レイノルズ応力で上式の未知変数   → モデル化の必要がある。  瞬間速度方程式 !U_i !t + ! " u _i !t + !U_iU_j !x_j + !U_iu " _j !x_j + ! " u _iU_j !x_j + ! " u _iu " _j !x_j = # !P !x_i # ! " p !x_i +$ !2U_i !x_j!x_j +$ !2u " _i !x_j!x_j !U_j !x_j + ! " u _j !x_j = 0  瞬間速度方程式から平均場方程式を引いて得られる揺動場方程式 ! " u _i !t + !U_iu " _j !x_j + ! " u _iU_j !x_j + ! " u _iu " _j !x_j # ! " u _iu " _j !x_j = # ! " p !x_i +$ !2u " _i !x_j!x_j ! " u _j !x_j = 0 乱流エネルギーの輸送方程式 D Dt 1 2 u i! u ! i " # $ % = & ! u iu ! j 'U_i 'x_j &( ' ! u _i 'x_j ' ! u _i 'x_j & ' 'x_j 1 2 u i! u i! u ! j +)i j p ! ! u i &( ' 'x_j 1 2 u ! iu i! " # * $ % + レイノルズ応力の輸送方程式 D ! u _iu ! _j Dt = " ! u ju m! #U_i #x_m " u i! u ! m #U_j #x_m + p ! # ! u _j #x_i + p ! # ! u _i #x_j " 2$ # ! u _i #x_m # ! u _j #x_m      ! " # u _iu # _ju _m# "x_m ! " p # # u _j "x_i ! " p # # u _i "x_j +$ "2 u # _iu # _j "x_m"x_m

 ○把握困難なほど非常に多くのモデルが提案されている。

E( k)

k

! k

"5/ 3 !"#$% &'() *+,() -.() /0() RANS の計算領域  ○レイノルズ応力を平均量を用いてモデル化する   → 渦粘性モデル     標準 K-ε モデルが代表的である。     渦粘性表現に関して線形モデルと非線形モデルに分けることができる。     予測能力は高くない。     計算安定性は高い。現在の実用面での主流。  ○レイノルズ応力の輸送方程式をモデル化する。   → 応力方程式モデル、代数応力モデル     渦粘性モデルに比べれば高い予測能力と有する。     計算安定性は乏しい。  ○壁関数で壁に接続するモデルを高レイノルズ数型モデルと呼ぶ。  ○壁まで適用できるモデルを低レイノルズ数型モデルと呼ぶ。  ○一般的には計算負荷は非常に軽い。  ○ RANS 自体は平均操作のため積分スケール程度までしか挙動をとらえられない   ので、非定常性（音）や３次元性（構造）では破綻する。 4.3.1 K-ε モデル  線形渦粘性表現 ! uiu!j = 2 3K"ij #$T %Ui %xj +%Uj %xi & ' ( ) * +

!T = 0.09 K2 "  乱流エネルギーと散逸率の輸送方程式 !K !t + Uj !K !xj = PK "# + ! !xj $ + $T %K & '( ) *+ !K !xj , -. /. 0 1 . 2. !" !t + Uj !" !xj = 1.44 " KPK # 1.92 "2 K + ! !xj $ +$T %_" & '( ) *+ !" !xj , -. /. 0 1 . 2. PK = ! ui"u"j #Ui #xj 4.3.2 SSG モデル  応力方程式 ! ui"u"j !t + Um ! u"iu"j !xm = Pij + #ij $%ij + ! !xm & + CDH K u"_mu_l" % ' () * +, ! ui"u"j !x/ -. / 0/ 1 2 / 3/  厳密な生成項 Pij = ! u"iu"m #Uj #xm ! u"ju"m #Ui #xm  再分配項のモデル表現 !ij "#ij = "# C1aij + $C1 aimamj " 1 3almaml%ij & '( ) *+ , -. / 0 1+ C2 P_Ka_ij + C3K 2Ui 2xj +2Uj 2xi & ' ( ) * + +C4K aimSmj + ajmSmi ! 2 3almSml"ij # $% & '(+ C5K a

(

imWmj + ajmWmi

)

! 2 3)"ij !" !t + Uj !" !xj = 1.44 " KPK # 1.92 "2 K + ! !xj $ +$T %_" & '( ) *+ !" !xj , -. /. 0 1 . 2. 第５章 計算手法 5.1 テイラー展開  ある距離!x はなれた関数値を評価する方法 f (x + !x) = f (x) + !x"f "x+ !x2 2 "2f "x2 + !x3 6 "3f "x3 +! = !xm m! m=0 #

$

" m f "xm 5.2 計算法の評価  流れの数値計算のプログラムは適合性、安定性、収束性を満足する必要がある。  ○ 適合性

  時間間隔と空間間隔を縮めていくと打ち切り誤差の項が消えて、元の微分方程式 に近づいていくこと。  例 波動方程式 !u !t " c !u !x = 0   テイラー展開 !u !t = u_in+1 " ui n #t " #t 2 !2u !t2! !u !x = u_i+1n " ui n #x " #x 2 !2u !x2!   差分化された方程式 ui n+1 ! ui n "t ! c u_i+1n ! ui n "x = "t 2 #2u #t2 ! c"x 2 #2u #x2! = O("t,"x)       # !x,!t"0# # # # 0"  ○ 安定性   打ち切り誤差や丸め誤差（コンピューターで計算する際に避けられない誤差）な  どが成長しないスキームであるかどうか。  ○ 収束性   時間間隔と空間間隔を十分小さく取った場合、差分方程式の数値解が元の微分方 程式の解析解に近づいていくかどうか。 5.3 時間発展  時間発展１階微分方程式 !" t

( )

!t = f

(

" t

( )

)

 離散化式 ! t

( )

n+1 "! t

( )

n #t = f ! t l

( )

( )

! t

( )

n+1 =! t

( )

n + "tf

( )

! t

( )

l  ○ 陽解法    l = n というように過去のデータのみを使って時間発展させていく方法    過去量は既知量であるため、代入操作だけで計算できる簡便な方法    計算安定性は比較的低いため、時間間隔を大きくはとれない  ○ 陰解法    l = n +1というように現在の時刻を含めて時間発展させていく方法    行列計算等を必要とするため、時間を１スッテプ進めるためにかかる計算コス

   トは陽解法に比べると大きい。    計算安定性が高く、時間間隔を大きくはとれる  ○ 半陰解法    陽解法と陰解法の両方を使って時間発展させていく方法 5.3.1 陽解法  ○ Euler 陽解法 ! t

( )

n+1 =! t

( )

n + "tf

(

! t

( )

n

)

  １次精度計算 精度評価 ! t

( )

n+1 "! t

( )

n #t " f ! t n

( )

(

)

= ! t

( )

n + #t$! t n

( )

$t + 1 2#t 2$ 2 ! t

( )

n $t2 +!" ! t n

( )

#t " $! t

( )

n $t         = 1 2!t "2# t

( )

n "t2 +! = O !t

( )

 ○ Adams-Bashforth 法   2 次精度 Adams-Bashforth 法 ! t

( )

n+1 =! t

( )

n +3 2"tf ! t n

( )

(

)

#1 2"tf ! t n#1

( )

(

)

精度評価 ! t

( )

n+1 "! t

( )

n #t " 3 2 f ! t n

( )

(

)

+1 2 f ! t n"1

( )

(

)

= ! t

( )

n + "t#! t n

( )

#t + 1 2"t 2# 2 ! t

( )

n #t2 + 1 2"t 3# 3 ! t

( )

n #t3 +!$ ! t n

( )

"t $ 3 2 #! t

( )

n #t + 1 2 #! t

( )

n$1 #t = 1 2!t 2" 2 # t

( )

n "t2 + 1 2!t 3" 3 # t

( )

n "t3 ! !t $ 1 2 "# t

( )

n "t + 1 2 "# t

( )

n "t $ !t "2# t

( )

n "t2 + 1 2!t 2" 3 # t

( )

n "t3 +! % & ' (' ) * ' +' = 3 4!t 2" 3 # t

( )

n "t3 +! = O !t 2

( )

  4 次精度 Adams-Bashforth 法 ! t

( )

n+1 =! t

( )

n +55 24"tf ! t n

( )

(

)

#59 24"tf ! t n#1

( )

(

)

+37 24"tf ! t n#2

( )

(

)

# 9 24"tf ! t n#3

( )

(

)

 ○ Runge-Kutta 法   2 次精度 Runge-Kutta 法   ●修正 Euler 法

!*=! t

( )

n +1 2"tf ! t n

( )

(

)

! t

( )

n+1 =! t

( )

n + "tf

( )

!*   ●改良 Euler 法 !*=! t

( )

n + "tf

(

! t

( )

n

)

! t

( )

n+1 =! t

( )

n +1 2"t f ! t n

( )

(

)

+ f

( )

!*

{

}

  4 次精度 Runge-Kutta 法 !1=! t

( )

n + "tf

(

! t

( )

n

)

!2 =! t

( )

n + "tf ! t n

( )

+!1 2 # $ % & ' ( !3=! t

( )

n + "tf ! t n

( )

+!2 2 # $ % & ' ( !4 =! t

( )

n + "tf

( )

!3 ! t

( )

n+1 =1 6 ! 1 + 2!2+ 2!3+!4

(

)

5.3.2 陰解法  ○ Euler 陰解法（完全陰解法） ! t

( )

n+1 =! t

( )

n + "tf

(

! t

( )

n+1

)

! t

( )

n+1 " #tf

(

! t

( )

n +1

)

=! t

( )

n   １次精度計算 精度評価   ! t

( )

n+1 "! t

( )

n #t " f ! t n+1

( )

(

)

= ! t

( )

n+1 "! t

( )

n+1 " #t$! t n+1

( )

$t " 1 2#t 2$ 2 ! t

( )

n+1 $t2 +! #t " $! t

( )

n+1 $t         = ! 1 2"t #2$ t

( )

n+1 #t2 +! = O "t

( )

 ○ Crank-Nicolson 法 ! t

( )

n+1 =! t

( )

n +1 2"t f ! t n

( )

(

)

+ f

(

! t

( )

n+1

)

{

}

! t

( )

n+1 "1 2#tf ! t n+1

( )

(

)

=! t

( )

n +1 2#tf ! t n

( )

(

)

  2 次精度計算 精度評価 ! t

( )

n+1 "! t

( )

n #t " 1 2 f ! t n

( )

(

)

"1 2 f ! t n+1

( )

(

)

= ! t

( )

n + "t#! t n

( )

#t + 1 2"t 2# 2 ! t

( )

n #t2 +!$ ! t n

( )

"t $ 1 2 #! t

( )

n #t $ 1 2 #! t

( )

n+1 #t = 1 2!t 2" 2 # t

( )

n "t2 + 1 2!t 3" 3 # t

( )

n "t3 +! !t + 1 2 "# t

( )

n "t $ 1 2 "# t

( )

n+1 "t = 1 2!t 2" 2 # t

( )

n "t2 + 1 2!t 3" 3 # t

( )

n "t3 +! !t + 1 2 "# t

( )

n "t $ 1 2 "# t

( )

n "t + !t "2# t

( )

n "t2 + 1 2!t 2" 3 # t

( )

n "t3 +! % & ' (' ) * ' +' = 1 4!t 2" 3 # t

( )

n "t3 +! = O !t 2

( )

5.4 空間解析法  ○ スペクトル法・擬スペクトル法   従属変数を適当な関数で級数展開を行い、その展開係数の従う方程式を解く方法   がスペクトル法である。概念的には格子点を考慮する必要はなく、離散化も無用   である。   特に三角関数によるフーリエ級数展開、チェビシェフ多項式によるチェビシェフ   級数展開による解析が主流である。   Taylor 展開における打ち切り誤差に無縁なため非常に高い精度で計算が可能で   あるが、対象の流れ場に適合する関数展開の選択が通常では困難が伴う。   DNS および LES 向きである。  ○ 有限差分法   格子点上に変数を配置し、数学的表現である微分を離散化された有限個の点にお   ける従属変数を用いた代数式（Taylor 展開を基盤とする）で計算していく。   Taylor 展開の打ち切り誤差を自分の必要とする精度まで簡単に高精度化するこ   とができるが、複雑化は免れない。流れ場の適用性も汎用的である。保存性を確   保する場合には注意が必要である。   DNS および LES 向きである。  ○ 有限体積法   解析方程式を有限な体積（コントロールボリューム）にわたって積分を実行して、

  コントロールボリューム界面からのフラックスのバランスを考慮して計算して   いく。   コントロールボリューム界面からの出入りを考慮しているため、保存性を常に確   保することが可能である。しかし、積分表現の精度をあげられず、2 次精度計算   が限界である。   RANS 向きである。  ○ 有限要素法   変数を格子頂点に配置し、その頂点を囲む要素（面や体）における変数分布を頂   点の変数から内挿する。変数分布を支配方程式に代入し、重み関数かけて積分す   ると、頂点上の変数の関係式が導出できる。この方程式を解析していく。   計算精度をあげるのは困難を伴う。しかし、複雑形状の流れ場の解析への適用性   は非常に高い。 5.5 スペクトル法・擬スペクトル法   スペクトル法と擬スペクトル法の違いは非線形項の取り扱い方による。  Navier-Stokes 方程式   

!u

_i

!t

= "u

j

!u

_i

!x

_j

"

!p

!x

_i

+

#

!

2

u

_i

!x

_j

!x

_j   

!u

_j

!x

_j

= 0

 周期境界条件   

u

i

(

x

m

+ 2!e

m

,t

)

= u

i

(

x

m

,t

)

 フーリエ変換   

u

_i

_{( )}

x,t

=

ˆu

i

( )

k,t

k_z=!" "

#

k_y=!" "

#

k_x=!" "

#

e

ikxx+ikyy+ikzz   

p x,t

_{( )}

=

ˆp k,t

_{( )}

k_z=!" "

#

k_y=!" "

#

k_x=!" "

#

e

ikxx+ikyy+ikzz  波数空間での Navier-Stokes 方程式   

d

dt

+

!

k

2

"

#

$

%

u

ˆ

i

= &ik

i

p & u

ˆ

j

'u

i

'x

_j   

i k

j

u

ˆ

j

= 0

 非線形項

  

ˆ

f

_i

! "u

j

#u

i

#x

_j  変形   

d

dt

+

!k

2

"

#

$

%

u

ˆ

i

= &ik

i

p + ˆ

ˆ

f

i  

i k

i_をかける   

i k

i

d

dt

+

!k

2

"

#

$

%

u

ˆ

i

= &ik

i

i k

i

p + i k

ˆ

i

ˆ

f

_i  圧力を解く   

p = !

ˆ

i k

_i

k

2

ˆ

f

_i  上式に代入   

d

dt

+

!

k

2

"

#

$

%

u

ˆ

i

= ˆ

f

i

&

k

_i

k

_j

k

2

ˆ

f

_j 非線形項の取り扱い方 ○スペクトル法  非線形項を単純にフーリエ変換して和則にしたがって計算を実行する。   

ˆ

f

_i

! "ik

_j

u

ˆ

_j

₍

m,t

₎

u

ˆ

_i

_{( )}

n,t

m +n =k

#

 利点 理論では有効。  欠点 計算時間がかかる。 ○擬スペクトル法  フーリエ変換を利用して、実空間で非線形項を計算し、さらにフーリエ変換で波数 空間に戻り、計算を進める。      逆変換       積を取る    変換  

u

ˆ

j_→

u

j        

u

_j

!u

i

!x

_j     →   

u

_j

!u

i

!x

_j  

i k

j

u

ˆ

i_→

!u

i

!x

_j  利点 計算時間が短縮できる。  欠点 特になし。

 （注意）擬スペクトル法とは非線形項で上述の和則に従わない部分（エリアイジング     エラー）を取り込んだものを示唆する場合もある。しかし、通常は 2/3 ルー     ル等の処理でこのエラーは擬スペクトル法でも除去している。 5.6 有限体積法   次の図にあるようなコントロールボリュームにわたる積分オペレーター ! =_"wedx_"sndy_"btdz  を Navier-Stokes 方程式にオペレートする。   例えばu の方程式 !u !t + !uu !x + !uv !y + !uw !z = " !p !x+# !2u !x2 +# !2u !y2 +# !2u !z2  では !u_P

!t dxdydz + u

(

eue" uwuw

)

dydz + u

(

nvn" usvs

)

dxdz + u

(

twt" ubwb

)

dxdy = " p

(

_e" p_w

)

dydz + # !u !x _e " !u !x_w $ % & ' ( ) dydz + # !u !y _n" !u !y_s $ % & ' ( ) dxdz + # !u !z_t" !u !z_b $ % & ' ( ) dxdy  となる。このレベルで表現は Taylor 展開の 2 次精度の表現となっている。  粘性拡散項に表れている微分表現は 2 次精度中心差分法により !u !x_e = u_E" u_P dx  を適用する。これにより !u_P

!t dxdydz + u

(

eue" uwuw

)

dydz + u

(

nvn" usvs

)

dxdz + u

(

twt" ubwb

)

dxdy = " p

(

_e" p_w

)

dydz + # u

(

_E " 2u_P + u_W

)

dydz dx +# u

(

N " 2uP + uS

)

dxdz dy +# u

(

T " 2uP + uB

)

dxdy dz  となる。小文字の添え字は補間法による評価が必要である。

P

b

t

e

w

s

n

dx

dy

dz

コントロールボリューム（大文字：定義点、小文字：補間点）  対流項の補間法   方程式を線形化するため u_eu_e ! u_eoldu_e  として速度の片方を繰り返し計算の一回前の値を利用する。  ○ 1 次精度風上補間 u_eoldu_e = uP, ue old > 0 u_E, u_eold < 0 ! " #   流れ場のレイノルズ数自体が不明になる人工粘性が生じるので、高い安定性を有   してはいても決して利用してはいけない方法である。名前をかえていまだに秘   かに生き残っているので精力的に撲滅する必要がある。  ○ 2 次精度中心補間 u_eoldu_e = u_eolduP + uE 2   精度は高いが、頂点衝突がある流れ場では数値振動（wiggle）が生じて計算でき   なくなる。  ○ 3 次精度風上補間（QUICK） u_eoldu_e = 3uE / 8 + 3uP / 4 ! uW / 8, ue old > 0 !u_EE / 8 + 3u_E / 4 + 3u_P/ 8, u_eold < 0 " # $   精度が高く、数値振動を防ぐことができ、安定性が 2 次精度中心補間より高い。 5.7 有限差分法  格子点上に変数を配置し、数学的表現である微分を離散化された有限個の点にお ける従属変数を用いた代数式（Taylor 展開を基盤とする）で計算していく。

x

x-2dx

x-dx

x+dx

x+2dx

Navier-Stokes 方程式では 1 階微分と 2 階微分の差分表現が必要になる。上図のよ うにターゲットの格子点の外にそれぞれ 4 点考慮して得られる 1 階微分と 2 階微分 の差分表現を導出する。それぞれの点の値はターゲットの格子点の量による Taylor 展開により u(x+ 2dx) = u(x) + 2dx!u !x+ 2dx 2! 2 u !x2 + 4dx3 3 !3 u !x3 + 2dx4 3 !4 u !x4 + 4dx5 15 !5 u !x5 + 4dx6 45 !6 u !x6 +!

u(x+ dx) = u(x) + dx!u !x + dx2 2 !2 u !x2 + dx3 6 !3 u !x3 + dx4 24 !4 u !x4 + dx5 120 !5 u !x5 + dx6 720 !6 u !x6 +! u(x)= u(x)

u(x! dx) = u(x) ! dx"u

"x + dx2 2 "2 u "x2 ! dx3 6 "3 u "x3 + dx4 24 "4 u "x4 ! dx5 120 "5 u "x5 + dx6 720 "6 u "x6 +!

u(x! 2dx) = u(x) ! 2dx"u

"x+ 2dx 2" 2 u "x2 ! 4dx3 3 "3 u "x3 + 2dx4 3 "4 u "x4 ! 4dx5 15 "5 u "x5 + 4dx6 45 "6 u "x6 +! で表される。 コンパクトさを考慮し、1 階微分の差分表現は以下のように書くことができる。 1 次精度 u(x+ dx) ! u(x) dx = "u "x+ dx 2 "2 u "x2 +! u(x)! u(x ! dx) dx = "u "x! dx 2 "2 u "x2 +! 2 次精度 u(x+ dx) ! u(x ! dx) 2dx = "u "x+ dx2 12 "3 u "x3 +! 3 次精度

!u(x + 2dx) + 6u(x + dx) ! 3u(x) ! 2u(x ! dx)

6dx = "u "x ! dx3 12 "4 u "x4 +!

2u(x+ dx) + 3u(x) ! 6u(x ! dx) + u(x ! 2dx)

6dx = "u "x + dx3 12 "4 u "x4 +! 4 次精度

!u(x + 2dx) + 8u(x + dx) ! 8u(x ! dx) + u(x ! 2dx)

12dx = "u "x! dx4 30 "5 u "x5 +! 2 階微分は以下のようになる。 2 次精度

u(x+ dx) ! 2u(x) + u(x ! dx)

dx2 = "2 u "x2 + dx2 12 "4 u "x4 +! 4 次精度

!u(x + 2dx) + 16u(x + dx) ! 30u(x) + 16u(x ! dx) ! u(x ! 2dx)

12dx2 = "2 u "x2 ! dx4 90 "6 u "x6 +! 2 階微分表現では非対称な 3 次精度の表現は一般には使われないので省略している。 それぞれ右辺第２項以降が差分誤差を意味している。精度は差分誤差における格子 間隔のべき乗に対応している。

 次に実際の計算すべき方程式に差分法を導入した場合の解説を行う。サンプルと して Navier-Stokes 方程式の特質を兼ね備えた１次元バーガース方程式を例に取る。 1 次元バーガース方程式は !u !t = " !uu !x +# !2 u !x2 である。拡散項（右辺第 2 項）は非線形移流項（右辺第 1 項）に比べて差分法依存 性は通常小さいので、非線形移流項について注目する。差分法には  ○ 保存型 !uu !x  ○ 非保存型 2u!u !x の取り扱いがある。数学的にはこの両者は同一のものであるが、数値計算では離散 化により数学的連続性・微分可能等の厳密性が成立していないため両者は異なる結 果を与える。保存型は外部からの流入が無ければ、積分量 dx

!

"uu_"x = uu

[ ]

_B1B2 = 0 は代数的に常に 0 として内部量が保存することを意味している。しかし、非保存型 では必ずしも 0 にならず、内部量が代数的に保存しないことを意味している。 5.7.1 １次精度風上差分法（非保存型）  風上法とは流れの上流の情報を利用して差分化する手法である。 u!u !x =

u(x)u(x)" u(x " dx)

dx + dxu(x) 2 !2 u !x2 +!

u(x)u(x+ dx) " u(x)

dx " dxu(x) 2 !2 u !x2 +! # $ %% & % % u(x)> 0 u(x)< 0    = u(x) ! !u(x)

(

)

u(x)! u(x ! dx) 2dx + u(x) ! u(x)

(

)

u(x+ dx) ! u(x) 2dx + dx u 2 "2 u "x2 +! この差分法を導入すると 1 次元バーガース方程式は !u !t = "2 u(x) " "u(x)

(

)

u(x)" u(x " dx) 2dx + u(x) " u(x)

(

)

u(x+ dx) " u(x) 2dx # $ % & ' (+) !2 u !x2 となるが、この方程式は結果として !u !t = "2u !u !x+

(

# + dx u

)

!2 u !x2 +! を解析していることになる。1 次精度風上差分法は右辺第 2 項が示しているように、 分子粘性に差分法による人工粘性が加わることを意味しており、常に安定な計算を

実行することはできるが、有限なコンピューターでは決して正しい計算を行うこと はできない。それゆえ、決して通常使ってはいけない計算方法である。また、この 手法は 2 次精度中心差分法とカップリングさせて名前を変え、ハイブリッド法、指 数法、Skew Upwind、EXQUISTE、流束分離法、TVD スキームなどに生き残って いる。注意すべし。ただし、超音速流のように、不連続な分布である衝撃波が生じ るケースでは今日の計算手法ではこのような方法論を使わざるおえないという限界 も存在している。 5.7.2 ２次精度中心差分法（非保存型）  非線形移流項に 2 次精度中心差分法を導入すると 2u!u !x = u(x) u(x+ dx) " u(x " dx) dx " dx2u(x) 6 !3 u !x3 +! !u !t = "u(x) u(x+ dx) " u(x " dx) dx +# !2 u !x2 となる。誤差を考慮すると !u !t = "2u !u !x+# !2 u !x2 " dx2 u 6 !3 u !x3 +! のように 3 階微分の影響が方程式に入ってくる。3 階微分項には分散効果があり、 絶対的な大きさを損なうことはないが、振動解を引き起こす原因となる。これが 2 次精度中心差分法における頂点衝突を含む流れ場での非物理的振動(wiggle)と対応し ている。 5.7.3 ２次精度中心差分法（保存型）  非線形移流項に以下のように 2 次精度中心差分法を導入する。 !uu !x = u2(x+ dx) " u2(x" dx) 2dx " dx2 12 !3 uu !x3 +! !u !t = " u2 (x+ dx) " u2 (x" dx) 2dx +# !2 u !x2 !u !t = " !uu !x +# !2 u !x2 + dx2 12 !3 uu !x3 +! 計算精度と保存性という点で比較的優れた計算スキームである。 5.7.4 ３次精度風上差分法（非保存型）  非線形移流項を以下のように風上差分法を導入すると 3 次精度風上差分法の一般 型がえられる。

2u!u

!x = u(x)

"u(x + 2dx) + 8u(x + dx) " 8u(x " dx) + u(x " 2dx) 6dx

+# u(x)u(x+ 2dx) " 4u(x + dx) + 6u(x) " 4u(x " dx) + u(x " 2dx) 2dx +dx3# u(x) 2 !4 u !x4 " dx4 u(x) 15 !5 u !x5 +! ここで η=1 は UTOPIA スキーム、η=3 は KK スキームと呼ばれる。この方法は 1 次精度風上差分法に比べると精度が良く、2 次精度中心差分法に比べると 4 階微分 項（超粘性拡散項）が加わるため安定性が高い。以前はこの安定スキームを利用し て直接数値計算を行ってきていた時期もあったが、超粘性も人工粘性であることに は変わりがないため、今日では有識者には DNS としては利用されていない。 5.7.5 ４次精度中心差分法（非保存型）  3 次精度風上差分法でとおくと 2u!u !x = u(x)

"u(x + 2dx) + 8u(x + dx) " 8u(x " dx) + u(x " 2dx)

6dx " dx4u(x) 15 !5 u !x5 +! 4 次精度の表現となる。  近年では更なる高次の中心差分法の利用などが進んでいるが、個人的には、微分 はローカルな概念であり高次にすると広範囲のデータが必要になるという点は矛盾 していると思われる。そこで、岡本研では IDO 法のようなコンパクト差分系の研究 を進めている。