スポンサーサイト
-- / -- / -- ( -- )
上記の広告は1ヶ月以上更新のないブログに表示されています。
新しい記事を書く事で広告が消せます。
--:--:-- | スポンサー広告 | page top↑
Excelで2次元熱伝導方程式を解く(その2)
2008 / 02 / 29 ( Fri )

Excelのシートを最大限活用します。ここで紹介する手法は、すでにインターネット上では、いくつか紹介されているものです。


まず、今回用意したExcelシートを、次に示します。

Excelsheet1.png

シートは大きく分けて、初期値を与える「領域1」、n 時間ステップの計算結果を保管する「領域2」、n+1時間ステップの計算ならびに計算結果を保管する「領域3」、そして計算結果を可視化する「領域4」の4つの領域からなります。

領域1では4つの変数、∆t = 0.25、x = 0.1、κ = 0.01そして初期時間、 t = 0を規定しています。

次に、領域2には、計算格子の各点におけるn 時間ステップの値が表示されます。ただし、ここでは計算格子点数を、IF = 11、JF = 11としました。計算初期値として、セル番号C7からC17まで1を代入し、その他のセルにはすべて0が代入されています。

領 域3がn+1時間ステップの計算をしている中枢部分であり、境界以外の内点のセルで、差分近似式が解かれます。

たとえば、セル番号D22をクリックすると、

fx=D8+$B$4*$D$4*(E8+C8+D9+D7-D8*4)/$C$4/$C$4

という演算式が現れます。Excelでセルの行番号と列番号にそれぞれ$マークが付いた場合には、固定されたセル番号を参照しますので、$B$4、$C$4、$D $4はそれぞれ領域1のDelta t、Delta x、Kappaの値である0.25、0.1、0.01を参照します。
また、セル番号D8、E8、C8、D9、D7は、それぞれ、
uni,j , uni+1.j , uni-1,j , uni,j+1  ,  uni,j-1 の値が保 管されている領域2のセル番号に相当します。
したがって、たとえばセル番号D22においては、格子点 i = 2、j = 2における値が計算されます。

一方、境界では、ノイマン境界条件が与えられている上辺と底辺に相当する格子点において計算が必要になります。

ここでは、計 算する差分近似された式は陽解法ですので、n+1 時間ステップで計算された1つ内点の値を用いて外挿します。たとえば、上辺 境界に位置するセル番号D21の計算式は、fx= D22と置きます。

境界と内点の計算式を用いた計算がすべて終了した段階で、領域3にはn+1時間ステップの計算結果 が保管されます。

スポンサーサイト

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

21:54:32 | 計算数理科学 | トラックバック(0) | コメント(0) | page top↑
Excelで2次元熱伝導方程式を解く(その1)
2008 / 02 / 28 ( Thu )
今回から何回かに分けて、2次元熱伝導方程式を Excel を用いて解く方法について紹介します。

1次元熱伝導方程式については、すでにExcelを用いて計算した例を紹介しました。1次元の場合には、画面上でExcelを手動で操作しながら簡単に計算することができます。しかしながら、2次元の場合には無理かもしれません。

そこで、簡単なマクロプログラムを追加するだけで計算できる方法について紹介します。ここでの内容は、すでに、ある学会から依頼されて執筆し掲載された原稿に基づいています。

まず、2次元熱伝導方程式を次のように定義します。

u/∂t = κ(2u/∂x2 + ∂2u/∂y2)

これを陽解法で差分近似すると、

( un+1i,juni,j)/t = κ {(ui+1.j - 2ui,jui-1,j) / (∆x)2 + (ui,j+1 - 2ui,j ui,j-1 ) / (∆y)2 }

簡単にするため、∆x = ∆y =κt/ (x )2 とおけば、

un+1i,j = uni,j + (uni+1.juni-1,j + uni,j+1uni,j-4uni,j)

のように比較的簡単な式を導出することができます。

具体的な熱伝導問題として、次の図に示す境界条件を設定します。

境界条件

クランク・ニコルソン法を、ガウス・ザイデル法やSOR法により数値計算する手法を紹介しましたが、マクロプログラムにすると結構複雑になり、Excel 本来の利点が生かせないので、ここでは最も単純な陽解法による方法を紹介していきます。

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

21:31:08 | 計算数理科学 | トラックバック(0) | コメント(0) | page top↑
熱伝導方程式をSOR法で解く
2008 / 02 / 27 ( Wed )
1次元熱伝導方程式をクランク・ニコルソン法で差分近似して、ガウス・ザイデル法を適用した式は、
(1+
ℓ)(un+1i ) m+1= {(un+1i+1 )m(un+1i-1 )m+1} / 2 + di

のようになりました。

ラプラス方程式の反復法の場合と同じく、
ガウス・ザイデル法を、SOR法に改良することができます。SOR法は、ガウス・ザイデル法の値を過大評価するというものでしたので、その式は、

(1+ ℓ)(un+1i ) m+1= (1 - ω ) (1+ ℓ)(un+1i ) m + ω [ℓ{(un+1i+1 )m(un+1i-1 )m+1]} / 2 + di ]

となります。

一般的に、1次元のみならず、2次元、3次元の熱伝導方程式は、クランク・ニコルソン法により差分近似して、SOR法による反復計算により解を求めます。

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

21:13:30 | 計算数理科学 | トラックバック(0) | コメント(0) | page top↑
熱伝導方程式をガウス・ザイデル法で解く
2008 / 02 / 26 ( Tue )
1次元熱伝導方程式をクランク・ニコルソン法で差分近似して、ヤコビ法を適用した式は、
(1+
ℓ)(un+1i ) m+1= {(un+1i+1 )m(un+1i-1 )m} / 2 + di

のようになりました。

ラプラス方程式の反復法の場合と同じく、ヤコビ法を、ガウス・ザイデル法に改良することができます。
ガウス・ザイデル法では、すでに計算済みの値を有効活用するというものでしたから、

すなわち、

(1+ ℓ)(un+1i ) m+1= {(un+1i+1 )m(un+1i-1 )m+1} / 2 + di

と簡単に式を導出することができます。

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

20:20:18 | 計算数理科学 | トラックバック(0) | コメント(0) | page top↑
熱伝導方程式をヤコビ法で解く
2008 / 02 / 25 ( Mon )

1次元熱伝導方程式にクランク・ニコルソン法を適用した式は、

un+1iuni)/t = κ {(uni+1 -2uni uni-1 ) / (x )2 + (un+1i+1 -2un+1i un+1i-1 ) / (x )2}/2

でした。=κt/ (x )2 として、さらに変形すると、

un+1i = (un+1i+1 -2un+1i un+1i-1 ) / 2 + di

ただし、

di = uni + (uni+1 -2uni uni-1 ) / 2

になります。

ラプラス方程式の反復計算に用いた
時間ステップ n における解は、収束解が得られるまでは意味のないものでしたが、熱伝導方程式における時間ステップ n により求まった解は実時間の解を与えています。したがって、熱伝導方程式に反復法を適用するためには、時間ステップ n とは別に、反復法のための計算ステップが必要になります。

反復法の時間ステップをm として、上記式にヤコビ法を適用すると、

(1+
ℓ)(un+1i ) m+1= {(un+1i+1 )m(un+1i-1 )m} / 2 + di

のような式を導出できます

di 
は既知量ですから、反復法の計算からは、はずされます。

un+1i の値を求めるために、m を1, 2, 3,・・・と増やして計算し、mm+1の値が同じになったとき、その値が求まるという式です。




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

22:36:47 | 計算数理科学 | トラックバック(0) | コメント(0) | page top↑
熱伝導方程式の直接法と反復法
2008 / 02 / 24 ( Sun )
クランク・ニコルソン法が連立1次方程式の逆行列計算に帰着されることを紹介しました。
このように各時間ステップ n ごとに逆行列を1回だけ計算する方法のことを、直接法(Direct Method)と呼びます。

連立1次方程式の計算は、中学・高校の数学の知識があればできます。

直接法には、いくつかの方法がありますが、
連立1次方程式の逆行列計算を体系化した方法を、ガウス消去法(Gauss Elimination Method)と呼びます。また、逆行列を上三角形、下三角形領域に分割して計算する方法は、LU分解法(Lower-upper Decomposition Method)と呼びます。

たしかに、1次元の熱伝導方程式の場合には、直接法で解くことができますが、2次元以上の熱伝導方程式になると、行列計算式自体が作れませんので、そのままでは直接法は用いることができません。さらに、計算格子点の数が多くなると、直接法は、多くの記憶領域を必要とするため、効率的ではありません。

したがって一般的には、熱伝導方程式はクランク・ニコルソン法を適用して、反復法により計算します。

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

10:20:26 | 計算数理科学 | トラックバック(0) | コメント(0) | page top↑
クランク・ニコルソン法は逆行列計算
2008 / 02 / 23 ( Sat )

熱伝導方程式

u/∂t = κ 2u/∂x2

をクランク・ニコルソン法により差分近似すると

un+1unj )/t = κ {(un+1 -2un un-1 ) / (x )2 + (un+1+1 -2un+1 un+1-1 ) / (x )2}/2

でした。いま、=κt/ (x )2 とおいて式を変形すると、

-un+1-1 /2+ (1+ℓ)un+1  + un+1+1 /2 = un-1 /2+ (1-ℓ)un  + un+1 /2

となります。

この式の右辺は、すでに既知の値です。既知量を簡単な記号に置き換えて、この式をよりわかりやすくすると、

a un+1-1 + b un+1  + c un+1+1  = d

ただし、右辺は格子点 j ごとに値が違いますので、右辺のみ d  とします。

この式は、

         ・

         ・

         ・
a u
n+1-2
+ b un+1 -1 + c un+1j   = d
-1

a un+1-1 + b un+1  + c un+1+1  = d

a un+1j  + b un+1 +1 + c un+1+2  = d +1

         ・

         ・

         ・

のように、格子点の数だけある連立1次方程式であることがわかります。 したがって、1次元熱伝導方程式にクランク・ニコルソン法を適用した差分近似式は、この連立1次方程式を解けば計算することができます。

連立1次方程式を解くとは、係数 a, b, cからなる行列Aと、既知量d  からなるベクトルから、un+1 からなるベクトルを求めることに等しいですから、結局、行列Aの逆行列を計算することに帰着されます。


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

16:11:14 | 計算数理科学 | トラックバック(0) | コメント(0) | page top↑
クランク・ニコルソン法の線形安定性
2008 / 02 / 22 ( Fri )
熱伝導方程式

u/∂t = κ 2u/∂x2

にクランク・ニコルソン法を適用した場合の、線形安定限界を同様に求めてみます。

クランク・ニコルソン法により差分近似された式は、

un+1unj )/t = κ {(un+1 -2un un-1 ) / (x )2 + (un+1+1 -2un+1 un+1-1 ) / (x )2}/2

のようになります。陽解法のときと同様に、

un  = Gn exp( j i θ)

を代入して整理すれば、

G = {1-κ t(1-cosθ)/ ( x )2} / {1+κ t(1-cosθ )/ ( x )2

になります。

この式は常に1以下の値になることがわかります。すなわち、クランク・ニコルソン法は線形安定性の制約を受けない、すなわち無条件安定な方法であることがわかります。

κ xは定数値ですから、t の値が無限大にまで大きくできるということになりますが、実際には、初期値や境界条件の制約により無限大にとれるわけではありません。

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

21:47:35 | 計算数理科学 | トラックバック(0) | コメント(0) | page top↑
フォン・ノイマンの線形安定性理論
2008 / 02 / 21 ( Thu )

フォン・ノイマン(von Naumann)の線形安定性理論については、以前その名前だけ紹介しました。

「熱伝導方程式を解く」や「クランク・ニコルソン法」では、1次元熱伝導方程式を陽解法で解くとき、t をある値より大きくすると数値振動が発生することを示しました。

ここでは、その理由をフォン・ノイマンの線形安定性理論から簡単に説明します。

いま、熱伝導方程式を

u/∂t = κ 2u/∂x2

とします。陽解法で差分近似した式は、

un+1j = un + κ t (un+1 -2unun-1 ) / ( x )2

のようになります。なお、先の説明のため、x 方向の格子点番号には、j を使用します。

フォン・ノイマンの線形安定性理論では、時間ステップnにおける格子点 j の解unが、有限の振幅を持った任意の位相の三角関数成分で与えられると仮定します。

すなわち、

un  = Gn exp( j i θ)

と定義されます。ただし、G=G(θ )は振幅で、増幅係数(amplitude factor)と呼ばれます。また、Gについている n は、べき乗です。θ =π / s ( s=±1, ±2, ±3, ±4,・・・)。i は虚数単位(-1の平方根 )。これを、差分近似式に代入すると、

G =1-2κ t (1-cos θ )/ ( x )

が得られます。フォン・ノイマンの線形安定性理論では、G の絶対値が、1より小さければ線形安定となります。したがって、

κ t /( x )1/2

であれば、上記熱伝導方程式の陽解法は、線形安定となります。この値が、1/2を超えてしまうと、「熱伝導方程式を解く」で紹介したような数値振動が発生していまします。

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

21:43:27 | 計算数理科学 | トラックバック(0) | コメント(0) | page top↑
クランク・ニコルソン法
2008 / 02 / 20 ( Wed )
「熱伝導方程式を解く」では、差分法による計算式を紹介しました。

熱伝導方程式をもう一度示せば、

 ∂u/∂t = κ 2u/∂x2

のようになります。

これを差分近似すれば、


un+1iuni)/t = κ (uni+1 -2uni uni-1 ) / (x )2

となり、さらに、

un+1i = uni + κ t (uni+1 -2uni uni-1 ) / ( x )2

のように変形することができます。

「熱伝導方程式を解く」では、この式を具体的に表計算ソフトで解きました。
また、計算結果をグラフで可視化しました。そのとき、
t をある値より大きくすると、数値振動が発生することを示しました。

n+1時間ステップの値を、
n 時間での値のみから計算する方法のことを、陽解法(Explicit Method)と呼びます。線形安定性理論から、陽解法には安定限界があることが証明されており、上記式の場合には、κ t /( x )2 が0.5 を超えると安定限界を超えて数値振動が発生するこが知られています。「熱伝導方程式を解く」で発生した数値振動は、まさにこれに該当します。

熱伝導方程式の安定限界を緩和する方法として、クランク・ニコルソン法(Crank-Nicolson Metho)が知られています。この方法では、熱伝導方程式を次の式のように差分近似します。

un+1iuni)/t = κ {(uni+1 -2uni uni-1 ) / (x )2 + (un+1i+1 -2un+1i un+1i-1 ) / (x )2}/2

これは、2階偏導関数
2u/∂x2を、n時間ステップとn+1時間ステップで差分近似して、それらを平均するという式です。ところで、このようにn+1時間ステップの値を計算するために、n+1時間ステップの値が必要にある方法のことを、陰解法(Implicit Method)と呼びます。たとえば、ガウス・ザイデル法の陰解法になります。

陰解法である
クランク・ニコルソン法を用いることにより、熱伝導方程式は理論的に無条件安定に計算することができます。t をいくらでも大きく取れるということを意味します。ただし、実際には境界条件や多次元化などによる制約から、t を無限大にできるというわけではありません。


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

21:59:19 | 計算数理科学 | トラックバック(0) | コメント(0) | 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↑
SOR法
2008 / 02 / 17 ( Sun )
ガウス・ザイデル法の計算式は、

un+1i,j = (uni+1.jun+1i-1,j + uni,j+1un+1i,j-1) / 4 

でした。この式の反復回数をさらに短縮することができます。

いま、
ガウス・ザイデル法の un+1i,j を、(un+1i,j)GSとおいて、その計算式を導出すれば、

un+1i,j = (1 - ω ) uni,j + ω (un+1i,j)GS

のように定義されます。
この式は、n 時間ステップにおける
uni,jガウス・ザイデル法で求まった値のω による線形結合により、un+1i,j を求めるという形の式になっています。

通常、線形結合においては、緩和係数
ω は、0 < ω < 1なのですが、ここでのωは、1 < ω になります。そのため、このωのことを、過緩和係数(Over-relaxation paremeter)と呼びます。

1 < ω にするということは、ガウス・ザイデル法の値を過大評価することを意味します。

理論的に、
ω < 2 ですが、経験的には、ω は1.5近傍の値が上限値になります。仮にω = 1.5 とすれば、この方法は、ヤコビ法に比べて反復回数が3分の1程度で済むことになります。この方法のことを、SOR法(Successive Over-relaxation Method)と呼びます。

ラプラス方程式の差分法には、SOR法が最も広く用いられています。

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

16:10:58 | 計算数理科学 | トラックバック(0) | コメント(0) | page top↑
ガウス・ザイデル法
2008 / 02 / 16 ( Sat )
ヤコビ法として、

un+1i,j = (uni+1.juni-1,j + uni,j+1uni,j-1) / 4 

を紹介しました。
この方法では、un+1i,j を求めるために、周りの4点から計算します。
n+1 時間ステップを計算するとき、n 時間ステップの値は既知ですので、任意のi,j 点における計算はすべて独立に、いわゆる並列計算が可能です。
すべての計算格子点にCPUを割り当てて計算すれば、1時間ステップの計算は1回の計算で済んでしまいます。

ただし、手元のパソコンでは並列計算は無理かもしれません。

ヤコビ法は、上記の式をひたすら繰り返して、最終的に、
n+1 時間ステップと n 時間ステップの値が同じになれば、解が求まったと判断されます。
しかしながら、かなりの反復回数が必要であることが知られています。

この反復回数を減らすことができる反復法として、ガウス・ザイデル法(Gaus- Seidel Method)があります。

ガウス・ザイデル法では、

ガウスザイデル

に示すような、ハイパーライン(Hyper Line)と呼ばれる、いわゆる掃引により計算します。すなわち、図中のハイパーラインに付いた大きな矢印の方向に向かって計算します。すると、ハイパーラインが通過した点での値は、updateされますから、赤で示したように
n+1 時間ステップの値になります。

これを式で記述すれば、

un+1i,j = (uni+1.jun+1i-1,j + uni,j+1un+1i,j-1) / 4 

のようになります。

すでに
updateされた点の値を有効活用することにより、反復回数を約半分に減らすことができます。


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

17:10:57 | 計算数理科学 | トラックバック(0) | コメント(0) | page top↑
ヤコビ法
2008 / 02 / 15 ( Fri )
ラプラス方程式の差分近似式については、すでに「ラブラス方程式を差分法で解く」のところで紹介しました。 いま、任意の計算格子点 i, j に対して、改めて差分近似式を定義すれば、

(ui+1.j - 2ui,jui-1,j) / (∆x)2 + (ui,j+1 - 2ui,j ui,j-1 ) / (∆y)2 = 0

と記述することができます。ここで、∆x = ∆y とすれば、結局、

ui,j = (ui+1.jui-1,j + ui,j+1ui,j-1) / 4

になります。

計算格子点 i, j における ui,j は、近接する4つの計算格子点における 
ui+1.j , ui-1,j , ui,j+1 , ならびにui,j-1 を足して4で割ることにより求めるという単純な四則演算式です。

この式にさらに、時間ステップ n を当てはめた式を、

un+1i,j = (uni+1.juni-1,j + uni,j+1uni,j-1) / 4 

のように記述します。

この式は、n の値を1, 2, 3, ・・・と増やして反復計算することにより、
un+1i,j の収束解を求める式で、このような反復計算する方法のことを、反復法(Relaxation Method)と呼びます。また、この反復式は、それを考えた研究者の名前で、ヤコビ法(Jabobi Method)と呼ばれ、極めて単純な式ですが、立派な反復法のひとつです。

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

21:00:56 | 計算数理科学 | トラックバック(0) | コメント(0) | page top↑
数値計算における物理的センス
2008 / 02 / 14 ( Thu )

「ラプラス方程式を手計算で解く」のところで紹介した問題の答えは、0.5でした。

実は、物理的なセンスがある人は、解く前に答えがわかります。

いま、無限に広がった一様な厚さのある金属平板を想像してみてください。

上面を温度1度、下面を温度0度としたらば、金属平板の厚さ方向中心部の温度は何度になるでしょうか?

よほど不均一な物質でない限り中間の温度になるはずです。すなわち、0.5度です。

ノイマン境界条件として与えた

u/∂x = 0

は、x 方向に u の勾配がないということですので、x 方向に u は同じ値になります。

計算数理科学が目指す数理モデルの構築は、もともと物理現象を数学的に記述することでしたから、もとの物理現象に立ち返って数理モデルを考察することが重要です。

 

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

21:16:54 | 計算数理科学 | トラックバック(0) | コメント(0) | page top↑
ラプラス方程式を手計算で解く(その4)
2008 / 02 / 13 ( Wed )

格子点(1,2), (2,2), (3,2) における差分近似式はそれぞれ、

-4u1,2+ u1,3u1,1 + 2u2,2 = 0

u3,2u1,2 + u2,3u2,1  - 4 u2,2 = 0

- 4u3,2+ u3,3u3,1 + 2u2,2 = 0  

のように求められました。

ところで、格子点(1,3),(2,3),(3,3)では、u = 1、格子点(1,1),(2,1),(3,1)では、u = 0 でしたから、

-4u1,21 + 0 + 2u2,2 = 0

u3,2u1,2 1 +0  - 4 u2,2 = 0

- 4u3,21 +0 + 2u2,2 = 0

となります。

結局、これら3つの式からなる連立1次方程式を解けば、

u1,2 = u2,2 = u3,2 = 0.5

となり、ラプラス方程式を手計算で解くことができました。

このことからわかりますが、ラプラス方程式の差分法による数値計算は、それをひも解いていけば、未知のu がある点の数だけの連立1次方程式の計算に最終的には帰着されます。言い方を変えれば、連立1次方程式から構成される行列の逆行列を求めることと同意義です。

ここでの問題はたった9点の差分計算でしたが、考え方は点の数が1億点になってもまったく同じです。ただし、1億点では手計算できませんので、代わりにコンピュータに計算を任せることになります。

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

18:18:51 | 計算数理科学 | トラックバック(0) | コメント(0) | page top↑
ラプラス方程式を手計算で解く(その3)
2008 / 02 / 12 ( Tue )

境界から点がはみ出さないようにするためには、片側差分近似式というものを新たに導出する必要があります。

「ラプラス方程式を手計算で解く(その1)」にある図には、黒塗りの点を2つ示しており、それぞれの格子点番号は、(3-1/2,2), (1+1/2,2)と定義されています。実際にはこれらの格子点はありません。あくまでバーチャルな格子点です

たとえば、格子点(1,2)では、x方向の差分近似で片側差分する必要があります。

そのために、(1,2)における 2u/∂x2 の片側差分近似式は、まず1階の偏導関数を用いて、

(∂2u/∂x2 )1,2 ={(∂u/∂x)1+1/2,2 - (∂u/∂x)1,2}/(x/2)

のように定義できます。ここでx= 1でしたので、

(∂2u/∂x2 )1,2 =2{(∂u/∂x)1+1/2,2 - (∂u/∂x)1,2}

になります。さらに、もともと境界条件として、(∂u/∂x)1,2= 0 でしたので、

(∂2u/∂x2 )1,2 =2(∂u/∂x)1+1/2,2

になります。(∂u/∂x)1+1/2,2u の二次精度中心差分近似で

(∂u/∂x)1+1/2,2 = (u2,2 - u1,2)/x = u2,2 - u1,2

のように計算されますので、結局

(∂2u/∂x2 )1,2 =2(u2,2 - u1,2

といった、ちょっと特殊な形の片側差分近似式が導出されます。

同様に、格子点(3,2)では、

(∂2u/∂x2 )3,2 =2{(∂u/∂x)3,2 - (∂u/∂x)3-1/2,2}

となります。境界条件として (∂u/∂x)3,2= 0 が与えられていますので、

(∂2u/∂x2 )3,2 = - 2(∂u/∂x)3-1/2,2

となり、結局

(∂2u/∂x2 )3,2 = - 2(u3,2 - u2,2

が得られます。

これより、(1,2)における差分近似式は

 2(u2,2 - u1,2)+ u1,3 - 2u1,2 u1,1 = 0

となり、同じ項をまとめれば、

-4u1,2+ u1,3u1,1 + 2u2,2 = 0

が得られます。

同様に、(3,2)における差分近似式は

- 2(u3,2 - u2,2)+ u3,3 - 2u3,2 u3,1 = 0

で、これを整理すれば、

- 4u3,2+ u3,3u3,1 + 2u2,2 = 0

となります。

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

21:37:59 | 計算数理科学 | トラックバック(0) | コメント(0) | page top↑
ラプラス方程式を手計算で解く(その2)
2008 / 02 / 11 ( Mon )

ラプラス方程式を手計算で解く(その1)では、主に境界条件を設定しました。

ここでは差分近似式を導出します。

この問題で、u が未知な計算格子点は、(1,2), (2,2), (3,2) の3点です。

「ラプラス方程式を差分法で解く」では、ラプラス方程式の差分近似式を

(u|x+dx - 2u|x u|x-dx ) / ( x )2 + (u|y+dy - 2u|y u|y-dy ) / ( y)2 = 0

のように与えました。

簡単にするため、x = y = 1と仮定します。

まず、(2,2)における差分近似式を求めてみます。そのために、上記差分近似式に具体的な計算格子点番号を当てはめますと、

u3,2 - 2u2,2 u1,2+ u2,3 - 2u2,2 u2,1  = 0

が求まります。ただし、式を簡単にするため | をはずしました。

同じ項で整理すれば、

u3,2u1,2 + u2,3u2,1  - 4 u2,2 = 0

次に、(1,2)はどうでしょうか。同様に格子点番号を当てはめますと、

u2,2 - 2u1,2 u0,2+ u1,3 - 2u1,2 u1,1  = 0

と求まりました。としたいところですが、よく見ると、u0,2 は求めることができません。もとの式のu|x-dx が左辺の境界からはみ出してしまいました。

二次精度差分近似式では、いま計算しようとする点と両隣の点の値が必要になるため、いま計算しようとする点が境界上にあると、両隣のいずれかの点が境界からはみ出してしまい、差分近似ができなくなります。

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

17:17:22 | 計算数理科学 | トラックバック(0) | コメント(0) | page top↑
ラプラス方程式を手計算で解く(その1)
2008 / 02 / 10 ( Sun )

ラプラス方程式の解析解、差分法、ならびに流れの方程式への応用などについて紹介してきました。ラプラス方程式はたいへん汎用性のある方程式で、多くの物理現象を模擬することができます。

差分法による数値計算がより身近に感じるように、手計算で解く方法を何回かに分けて紹介します。

その1として、まず具体的な問題を設定します。

次の図にあるような計算格子点がたった9点の問題を考えます。

Laplace手計算

たった9点の問題でも立派なラプラス方程式の差分計算ができます。

ラプラス方程式のような楕円型方程式の差分計算には、境界条件が欠かせません。そのため、楕円型方程式の問題は、境界値問題(Boundary Value Problem)と呼ばれます。

計算格子点にはぞれぞれ番号 i,j をつけています。ただしここでは、i = 1,2,3,  j = 1,2,3です。

境界条件として、

上辺にある、(1,3), (2,3), (3,3)には、u = 1

下辺にある、(1,1), (2,1), (3,1)には、u = 0

を与えます。このように未知変数であるu そのものを規定する境界条件のことを第一種境界条件、もしくは考えた研究者の名前で、ディリクレ境界条件(Dirichlet's Boundary Condition)と呼びます。これらの点の値はすでに求まっていますから、改めて計算する必要はありません。

一方、左辺と右辺にある、(1,2) ならびに( 3,2)には、

u/∂x = 0

を与えます。このように未知変数であるu の1階偏微分を規定する境界条件のことを第二種境界条件、もしくは考えた研究者の名前で、ノイマン境界条件(Neumann's Boundary Condition)と呼びます。Neumannとは、先に紹介したvon Neumannです。これら2点では、u 自体の値は未知ですので計算により求める必要があります。

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

18:14:31 | 計算数理科学 | トラックバック(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 / 06 ( Wed )

「ラプラス方程式の解析解」で紹介したラプラス方程式を差分法で解いてみましょう。

ラプラス方程式は、

2u/∂x2 + ∂2u/∂y2 = 0

で与えられます。これを中心差分で近似しますと、

(u|x+dx 2u|x u|x-dx ) / ( x )2 + (u|y+dy 2u|y u|y-dy ) / ( y)2 = 0

のようになります。簡単のため、x = y して、この式をさらに変形すると、

u|x (u|x+dxu|x-dx u|y+dyu|y-dy ) / 4

が得られます。これは、自分自身の値を隣接する4点を足して4で割って求めるという単純な計算式です。さらに反復計算を施せば、

un+1|x (un|x+dxun|x-dx un|y+dyun|y-dy ) / 4

のようになります。

いま、「ラプラス方程式の解析解」のところで紹介した問題の底辺の境界条件をf ( x )=1として実際に反復計算をしますと次のような計算結果を得ることができます。

ラプラス計算結果

計算に使用している計算点の数が、10x10=100点と少ないために、解が多少波打っている部分もありますが、比較的簡単に答えを求めることができます。

計算にはExcelを使用しました。若干のマクロプログラムとともにシート上の計算だけで、Excelでもラプラス方程式を差分法で解くことができます。

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

21:34:51 | 計算数理科学 | トラックバック(0) | コメント(0) | page top↑
ラプラス方程式の解析解
2008 / 02 / 05 ( Tue )

楕円型方程式の標準形であるラプラス方程式の解析解を求めてみましょう。

 

いま、次の図で示した正方形領域について、各辺にそれぞれ境界条件を与えます。

 

Laplace.png

 

 この問題は、典型的な偏微分方程式の計算問題ですが、変数分離法と呼ばれる方法により、解析解を求めることができます。計算の詳細についてはここでは省略します。詳細が知りたいときは、既存の偏微分方程式に関する教科書にはどれにでも同じことが書かれていますので、参照してください。

 

結局のところ、導出された解析解は次のようになります。

 

u (x,y)=2∑sinh((y-1))/sinh(-)*sin(nπx)∫f(z)sin(nπz)dz

 

ただし、∑ は、n = 1 から までの和、∫は、z = 0 から1までの積分です。また、

 

sinh(y) = ( ey - e-y )/2

 

です。かなり複雑な式になっていますが、上記の図に示した領域内の任意の点の解が求まりますので、たいへんエレガントな式です。

 

このように、ラプラス方程式の解析解は紙と鉛筆があれば求めることができます。ただし、仮に上記の図で示した境界条件のu = 0 のひとつでも0以外の値にしますと、途端に解析解を求めるのが極端にむずかしくなります。もしくは、もはや求めること自体できまくなります。そのときは、数値計算により求めるしかなくなります。

 

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

20:50:21 | 応用数学 | トラックバック(0) | コメント(0) | page top↑
偏微分方程式の型の分類
2008 / 02 / 04 ( Mon )

これまでに、移流方程式、拡散方程式、そしてナビエ・ストークス方程式と偏微分方程式を紹介しました。これらの偏微分方程式はその特徴から、大きく3種類に分類できます。

いま、次の二次元一般形で表された2階の偏微分方程式を考えます。

A2u/∂x2 +2B2u/∂xy+ C2u/∂y2 = f ( x, y, u, ∂u/∂x, ∂u/∂y )

ここで、A, B,ならびに C は高々x,y の関数で、右辺のf は高々1階の偏導関数からなる関数とします。

D = B2 - AC としますと、D の符号に応じて、

D > 0 なら、双曲型(Hyperbolic)

D = 0 なら、放物型(Parabolic)

D < 0 なら、楕円型(Elliptic)

という型に分類されます。なお、型の分類の詳細は少し複雑な話ですので、ここでは省略します。

それぞれの型には典型的な式の形があり、

双曲型では、

2u/∂x2 - ∂2u/∂y2 = 0

放物型では、

2u/∂x2 = ∂u/∂y

楕円型では、

2u/∂x2 + ∂2u/∂y2 = 0

などです。これらの式は標準形と呼ばれます。

 双曲型の式は、計算数理科学では波動方程式(Wave Equation)とも呼ばれます。

放物型の式はよく見ますと、yt になれば、「拡散方程式」で紹介した熱伝導方程式と同じであることがわかります。したがって、熱伝導方程式と呼ばれます。

そして、楕円型の式は、ラプラス方程式(Laplace Equation)と呼ばれます。また、この式の右辺がゼロでない場合には、ポアソン方程式(Poisson Equation)と呼ばれます。

これらの式はいずれも、物理現象を模擬する上でたいへん有用な式です。

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

21:25:29 | 応用数学 | トラックバック(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↑
| ホーム |
上記広告は1ヶ月以上更新のないブログに表示されています。新しい記事を書くことで広告を消せます。