スポンサーサイト
-- / -- / -- ( -- )
上記の広告は1ヶ月以上更新のないブログに表示されています。
新しい記事を書く事で広告が消せます。
--:--:-- | スポンサー広告 | 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↑
前ページ | ホーム | 次ページ
上記広告は1ヶ月以上更新のないブログに表示されています。新しい記事を書くことで広告を消せます。