スポンサーサイト
-- / -- / -- ( -- )
上記の広告は1ヶ月以上更新のないブログに表示されています。
新しい記事を書く事で広告が消せます。
--:--:-- | スポンサー広告 | page top↑
具体的な三次元ポテンシャル流れの数値計算
2008 / 02 / 19 ( Tue )
以前、とある教科書を分担執筆した際に、立方体周りのポテンシャル流れを数値計算したことがありますので、それをここでは簡単に紹介します。

格子点数は21x21x21とし、立方体はその表面がちょうど、i = j = k = 8  な らびに
i = j = k = 14 に位置するように設定します。境界条件は、入口境界i = 1 面 でφ = 0、 出口境界 i = IF 面 でφ = 1を与え、それ以外の外部境界および立方体表面には、法線方向 にφ/∂n=0 となるようノイマン境界条件を与えます。
なお、格子間隔は
∆x=∆y =∆z=0.1 で、過緩和係数は1.5とします。

「三次元ポテンシャル流れを反復法で解く」で導出した式をSOR法で反復計算しますと、ポテンシャル
φ は次の図のように求められます。

このポテンシャルを各座標方向に1階偏微分したものが、それぞれの方向の流速になりますから、
φ の値を単純に差分近似すれば、流速ベクトルを次の図のように求めることができます。

ここで作成した
Fortranによる計算プログラムは、研究室ホームページの下記場所に掲載してありますので、興味のある方はご自由にダウンロードしてください。

三次元ポテンシャル流れのFortran プログラム

スポンサーサイト

テーマ:自然科学 - ジャンル:学問・文化・芸術

21:40:08 | 数値流体力学 | トラックバック(0) | コメント(0) | page top↑
3次元ポテンシャル流れをSOR法で解く
2008 / 02 / 18 ( Mon )

ポテンシャル流れが、ラプラス方程式を解いて求められることを紹介しました。
ここでは具体的に、三次元のポテンシャル流れにSOR法を適用した際の式を導出します。

まず、ポテンシャル
φ の三次元ラプラス方程式は、

2φ /∂x2 + ∂2φ /∂y2 + ∂2φ /∂z2 = 0

になります。ただし、流速u, v, w はポテンシャルφ  

u = ∂φ /∂x , v = ∂φ /∂y , w = ∂φ /∂z

 のように定義されます。もともとこれを連続の式に代入すればラプラス方程式が得られます。

三次元ラプラス方程式を差分近似すると

(φi+1,j,k - 2φi,j,kφ i-1,j,k) / (∆x)2

+ (φi,j+1,k - 2φi,j,kφ i,j-1,k ) / (∆y)2

+ (φi,j,k+1 - 2φi,j,kφi,j,k-1 ) / (∆z)2= 0

のように多少長い式になります。z方向に新たに計算格子点k を用います。

これまでは、
∆x
, ∆yを簡単にして式を導出してきましたが、ここではこれらの値に忠実に式を導出します。

上記の式を、
φi,j,k の式に変形すれば、

φi,j,k = {(φi+1,j,kφ i-1,j,k) / (∆x)2 + (φi,j+1,kφ i,j-1,k ) / (∆y)2

+ (φi,j,k+1φi,j,k-1 ) / (∆z)2}/{2/(∆x)2+2/(∆y)2+2/(∆z)2}

のようになります。さらにSOR法を適用すれば、

φn+1i,j,k (1 - ω )φni,j,k+ω{(φni+1,j,kφn+1 i-1,j,k) / (∆x)2

+ (φni,j+1,kφn+1 i,j-1,k ) / (∆y)2 + (φni,j,k+1φn+1i,j,k-1 ) / (∆z)2}

/{2/(∆x)2+2/(∆y)2+2/(∆z)2}

となります。

この式を反復計算することにより、三次元ポテンシャル流れを数値計算することができます。

テーマ:自然科学 - ジャンル:学問・文化・芸術

21:15:20 | 数値流体力学 | トラックバック(0) | コメント(0) | page top↑
流れ関数
2008 / 02 / 09 ( Sat )

「渦度方程式」では、渦度を直接計算できる渦度方程式を紹介しました。

しかしながら、あいかわらず流速u,v の値が必要です。

これを解決するために、流れ関数(Steam Function)と呼ばれる新たなパラメータを導入します。流れ関数には ψ を用います。流速u,v は、ψの関数として次のように定義されます。

u = ψ/∂y,  v = - ψ/∂x

これを、渦度方程式に代入すれば、

ω/∂t + ∂ψ/∂y ・∂ω/∂x - ψ/∂x ・∂ω/∂y =1/Re(∂2ω/∂x2 + ∂2ω/∂y2)

一方、渦度は、

ω = v/∂x- u/∂y

でしたが、これに流れ関数で定義された流速を代入すると、

ω = - ∂(∂ψ/∂x)/∂x- ∂(∂ψ/∂y)/∂y

もう少し変形して、

2ψ/∂x2 + ∂2ψ/∂y2 = - ω

となり、結局、渦度方程式とこの式の2式を連立して解くことにより、二次元非圧縮性粘性流れを数値計算することができます。渦度方程式は時間進行法で解き、流れ関数のポアソン方程式は、「ラプラス方程式を差分法で解く」で紹介した方法で解きます。

通常、二次元非圧縮性粘性流れの基礎方程式は、連続の式と2つの運動方程式の計3つの式からなりますが、渦度と流れ関数を用いることで、計算する式の数を一つ減らすことができるため、計算量を削減することができます。

ただし、たとえば圧力の情報が要らないことが、逆に境界条件の設定に不都合だったり、三次元の渦度の定義が複雑なことから、三次元流れに拡張するのが難しいなど、問題点もいくつかあります。


テーマ:自然科学 - ジャンル:学問・文化・芸術

16:25:38 | 数値流体力学 | トラックバック(0) | コメント(0) | page top↑
渦度方程式
2008 / 02 / 08 ( Fri )

二次元非圧縮性粘性流れの運動方程式は、

u/∂t + uu/∂x + vu/∂y = -∂p/∂x +1/Re(∂2u/∂x2 + ∂2u/∂y2)

v/∂t + uv/∂x + v v/∂y = -∂p/∂y +1/Re(∂2v/∂x2 + ∂2v/∂y2)

のようになります。

さて、「ポテンシャル流れ」では、渦なし流れはラプラス方程式で表されることを示しました。

いま、流速 v の運動方程式をx でさらに偏微分した式から、流速 u の運動方程式をy でさらに偏微分した式を引いて整理すると、

ω/∂t + uω/∂x + v ω/∂y =1/Re(∂2ω/∂x2 + ∂2ω/∂y2)

のような式ができます。

ここで、ω はもともと、ω = ∂v/∂x - u/∂y で、二次元の渦度の式です。

「ポテンシャル流れ」では、ω =0 を仮定していましたが、二次元の二つの運動方程式から、渦度の方程式、すなわち、渦度方程式(Vorticity Equation)を導出することができます。これにより流れの渦度を直接計算することができます。

よく見ると、運動方程式にあった圧力項が渦度方程式にはありません。式の導出の過程で消えました。したがって、圧力場の情報がなくても渦度方程式は計算できます。しかしながら、相変わらず流速u, v の情報は必要です。

テーマ:自然科学 - ジャンル:学問・文化・芸術

21:53:31 | 数値流体力学 | トラックバック(0) | コメント(0) | page top↑
ポテンシャル流れ
2008 / 02 / 07 ( Thu )

ラプラス方程式が差分法で簡単に計算できることを紹介しましたが、流れ問題をラプラス方程式で近似することができます。

いま、流速u,v をポテンシャルと呼ぶ記号 φ で次のように定義します。

u = ∂φ/∂x  ,  v = ∂φ/∂y

ところで、流れの渦度 ω は二次元で、

ω=∂v/∂x - u/∂y

となりますが、これにポテンシャルで定義された流速を代入すると、

ω=∂(∂φ/∂y)/∂x - ∂(∂φ/∂x)/∂y = 0

です。すなわち、ポテンシャルで表わされた流れ場は渦なし流れになります。

さて、連続の式は、

u/x + ∂v/∂y  = 0

ですので、これに渦なし流れの流速を代入すると、

∂(∂φ/∂x/x + ∂(∂φ/∂y)/∂y  = ∂2φ/∂x2 + ∂2φ/∂y2 = 0

となります。

これより、渦なし流れはポテンシャルφ のラプラス方程式を解いて求めることができます。渦なし流れはポテンシャル流れ(Potential Flow)とも呼ばれます。

テーマ:自然科学 - ジャンル:学問・文化・芸術

22:18:48 | 数値流体力学 | トラックバック(0) | コメント(0) | page top↑
非圧縮性粘性流れの基礎方程式
2008 / 02 / 03 ( Sun )

無次元化された非圧縮性粘性流れの基礎方程式は、これまで紹介した式をまとめれば、結局次のようになります。

u/x + ∂v/∂y + ∂w/∂z = 0  

u/∂t + u u/∂x + v u/∂y + w u/∂z = -∂p/∂x +1/Re(∂2u/∂x2 + ∂2u/∂y2+ ∂2u/∂z2

v/∂t + u v/∂x + v v/∂y + w v/∂z = -∂p/∂y +1/Re(∂2v/∂x2 + ∂2v/∂y2+ ∂2v/∂z2

w/∂t + u w/∂x + v w/∂y + w w/∂z = -∂p/∂z +1/Re(∂2w/∂x2 + ∂2w/∂y2+ ∂2w/∂z2

T/∂t + u T/∂x + v T/∂y + w T/∂z = κ/ρCp(∂2T/∂x2 + ∂2T/∂y2+ ∂2T/∂z2

これらの式を連立して解くことにより、非圧縮性粘性流れを計算します。

 

実は、変数の数が方程式の数より1つ多いため、これらの式を解いても解は求まりません。圧力p が不定です。

 

p を求める方法には大きくわけて二つあります。一つは理想気体を仮定して、圧力の状態方程式 p = ρRT を合わせて計算する方法です。空気の遅い流れが、非圧縮性粘性流れと仮定して解くことができます。

 

それでは水の場合にはどうでしょう。理想気体の状態方程式は使えません。

 

連続の式が直接解けないの加えて、圧力場を如何に求めるかが、非圧縮性粘性流れの数値流体力学にとっては大きな問題です。           

テーマ:自然科学 - ジャンル:学問・文化・芸術

10:10:16 | 数値流体力学 | トラックバック(0) | コメント(0) | page top↑
エネルギー式
2008 / 02 / 02 ( Sat )

エネルギー式として

∂(ρh)/∂t + u ∂(ρh)/∂x + v ∂(ρh)/∂y + w ∂(ρh)/∂z = κ(∂2T/∂x2 + ∂2T/∂y2+ ∂2T/∂z2

を紹介しましたが、いま非圧縮性流れで、さらに hCpT という熱力学の関係式を導入しますと、

T/∂t + u T/∂x + v T/∂y + w T/∂z = κ/ρCp(∂2T/∂x2 + ∂2T/∂y2+ ∂2T/∂z2

となります。ただし、Cpは定圧比熱と呼ばれる物質固有の物性値です。温度を変数とした式に変形することができます。この式を見ると、流れがない場合には対流項はすべてなくなり、結局、熱伝導方程式と同じ形の式になることがわかります。

テーマ:自然科学 - ジャンル:学問・文化・芸術

10:29:03 | 数値流体力学 | トラックバック(0) | コメント(0) | page top↑
運動方程式
2008 / 02 / 01 ( Fri )

「ナビエ・ストークス方程式」で紹介した運動方程式のうち、x方向の流速 u の運動方程式は、

∂(ρu)/∂t + u ∂(ρu)/∂x + v ∂(ρu)/∂y + w ∂(ρu)/∂z

                              = -∂p/∂x +μ(∂2u/∂x2 + ∂2u/∂y2+ ∂2u/∂z2

 

でした。ただし、非圧縮性流れを仮定すれば、

 

u/∂t + u u/∂x + v u/∂y + w u/∂z

                              = -1/ρp/∂x +μ/ρ(∂2u/∂x2 + ∂2u/∂y2+ ∂2u/∂z2

 

となります。

 

運動方程式の左辺は、時間微分項(TIme-derivative Term)と対流項(Convection Term)、右辺は、圧力項(Pressure Term)と粘性項(Viscosity Term)からなっています。

非圧縮性流体は、圧力の高いところから低いところへ移流します。また粘性があり、物体壁などに引きずられ移流に影響を与えます。上の式は、これらの項のバランスから成り立っています。y方向、z方向についても同様です。

 

無次元化(Nondimensionalization)という操作をしますと、

 

u/∂t + u u/∂x + v u/∂y + w u/∂z

                              = -∂p/∂x +1/Re(∂2u/∂x2 + ∂2u/∂y2+ ∂2u/∂z2

 

となります。1/ρ がなくなり、μ/ρ 1/Re に変わりました。このときの、Re がレイノルズ数(Reynolds number)と呼ばれるもので、流体力学では一番有名な無次元パラメータです。レイノルズ数は慣性力と粘性力の比で定義されます。レイノルズ数が極めて大きい場合には、粘性力は慣性力に比べて取るに足りないと仮定されて、粘性項を省略した式を計算する場合もあります。ところで無次元化の操作ですが、説明が少し面倒ですので、ここでは省略します。

テーマ:自然科学 - ジャンル:学問・文化・芸術

21:57:33 | 数値流体力学 | トラックバック(0) | コメント(0) | page top↑
連続の式
2008 / 01 / 31 ( Thu )

流れの連続の式として、

ρ/∂t + ∂(ρu)/∂x + ∂(ρv)/∂y + ∂(ρw)/∂z = 0

を紹介しました。仮に流体が非圧縮性ですと、密度 ρ は変化しませんので、連続の式は次のように簡略化されます。

u/x + ∂v/∂y + ∂w/∂z = 0

1階の偏導関数のみからなる見た目にも簡単な偏微分方程式になります。しかしながら、この方程式を直接解くことができる数値計算手法は見当たりません。事実上、この方程式は直接解けないと考えた方がよさそうです。

実は、非圧縮性流れの連続の式がこのような形をしていることが、非圧縮性流れの数値流体力学の展開にとって極めて重要な意味を持ちます。

これまで紹介した、反応方程式、移流方程式、そして拡散方程式にはいずれも時間微分項、u/∂t が含まれており、差分法では

un+1 = unt f n

のように反復計算しました。このような計算手法を時間進行法(Time-marching Method)と呼びます。時間微分項がある方程式は時間進行法を用いて比較的簡単に解くことができます。しかしながら、この項がないと別の手法で解かなければなりません。連続の式を解く適当な手法が見当たらないというのが現状です。

テーマ:自然科学 - ジャンル:学問・文化・芸術

21:34:41 | 数値流体力学 | トラックバック(0) | コメント(0) | page top↑
ナビエ・ストークス方程式とは
2008 / 01 / 30 ( Wed )

計算数理科学として、反応方程式、移流方程式、そして拡散方程式を紹介してきました。それでは次の方程式はどうでしょう。

 

ρ/∂t+ ∂(ρu)/∂x + ∂(ρv)/∂y + ∂(ρw)/∂z = 0 

 

 これは、三次元流れの連続の式(Continuity Equation)です。t は時間、(x,y,z)は三次元直交座標、(u,v,w)は(x,y,z)方向の流速、ρは流体の密度です。すべて、1階の偏導関数から成っており移流方程式に似ています。

 

次はどうでしょう。

 

∂(ρu)/∂t + u ∂(ρu)/∂x + v ∂(ρu)/∂y + w ∂(ρu)/∂z

                                     = -∂p/∂x +μ(∂2u/∂x2 + ∂2u/∂y2+ ∂2u/∂z2

 

∂(ρv)/∂t + u ∂(ρv)/∂x + v ∂(ρv)/∂y + w ∂(ρv)/∂z

                                      = -∂p/∂y +μ(∂2v/∂x2 + ∂2v/∂y2+ ∂2v/∂z2

 

∂(ρw)/∂t + u ∂(ρw)/∂x + v ∂(ρw)/∂y + w ∂(ρw)/∂z

                                     = -∂p/∂z +μ(∂2w/∂x2 + ∂2w/∂y2+ ∂2w/∂z2

 

複雑すぎてよく見れないかもしれません。これらは比較的遅い流れの運動方程式(Momentum Equations)です。三次元の場合には3つの方程式からなります。pは圧力で、μは粘性係数です。とりあえず、これまでの知識から1階と2階の偏導関数からできているのがわかります。移流方程式と拡散方程式が混ざったような形をしています。ただし、厳密には2階微分項の部分はもう少し複雑になります。これについては改めて説明いたします。

 

最後に、

∂(ρh)/∂t + u ∂(ρh)/∂x + v ∂(ρh)/∂y + w ∂(ρh)/∂z

                                     = κ(∂2T/∂x2 + ∂2T/∂y2+ ∂2T/∂z2

 

は比較的遅い流れのエネルギー方程式(Energy Equation)です。hは比エンタルピー、Tは温度、κは熱伝導係数です。これも、移流方程式と拡散方程式が混ざったような形をしているのがわかります。特に右辺は、熱伝導方程式の拡散項と同じです。

 

ここで紹介した5つの方程式は、流体力学の完成形である、三次元ナビエ・ストークス方程式(Navier-Stokes Equations, 略してN-S 式)のひとつの表記形です。いろいろな仮定の導入で、N-S式の記述方法も違ってきます。ところで、一般的に非圧縮性のN-S式とは運動方程式そのものを指すことが多いようです。

 

数値流体力学では、流体力学の基本知識は重要です。ただ、とりあえずN-S式を構成する各項、すなわち、時間微分項(TIme-derivative Term)、対流項(Convection Term)、圧力項(Pressure Term)、そして粘性項(Viscosity Term)がどのような偏導関数から成り立っているかがわかれば、数値計算することはできます。

 

テーマ:自然科学 - ジャンル:学問・文化・芸術

21:29:00 | 数値流体力学 | トラックバック(0) | コメント(0) | page top↑
| ホーム | 次ページ
上記広告は1ヶ月以上更新のないブログに表示されています。新しい記事を書くことで広告を消せます。