熱伝導方程式を数値的に解く

1

今回は領域の縁で境界条件が与えられている熱伝導方程式を数値的に解く。

1.1

差分法による解法

これまでの講義と同様にまずは差分化によって解いてみる。

1.1.1

熱伝導方程式の差分化

以下の 2 次元熱伝導方程式を与えられた境界条件で解くことを考えよう。

∂ 2T∂T∂ 2T+=κ
∂t
∂x2
∂y 2!

(1)

ここで κ は熱伝導係数である。前回の波動方程式と同様に右辺の 2 回微分を差分で評価す
ると、∆x = ∆y = ∆ とおけば、R.H.S. =κ(Ti+1 j + Ti−1 j + Ti j+1 + Ti j−1 − 4Ti j )∆2(2)

左辺は一階微分なので時間方向のテーラー展開の式T (tk + ∆t, x, y) = T (tk , x, y) +∂T∆t + · · ·∂t(3)を用い、時間に関するグリッドの添え字を右肩に括弧つきで書くことにすると

∂T
∂t(k)

=T (k+1) − T (k)∆t(4)

が得られる。よって

T (k+1) − T (k)∆t
右辺も k で評価するべきなので、結局もとの熱伝導方程式は、
L.H.S. =
1  (k+1)κ  (k)
(k)
(k)
(k)
(k)
(k)
Ti j − Ti j = 2 Ti+1 j + Ti−1 j + Ti j+1 + Ti j−1 − 4Ti j
∆t
∆
(k+1)

これを Ti j(5)(6)

について解くと、(k+1)

Ti j(k)= Ti j +κ∆t  (k)
(k)
(k)
(k)
(k)
T+T+T+T−4Ti+1 ji−1 ji j+1i j−1i j∆2(7)この差分式は右辺に k + 1 の時間ステップの値を含んでいないので、やはり前回同様に反復する必要は無い。図は差分の様子を表しているが、y 方向は紙面に垂直方向なので割愛してある。ここで T 0 は初期の温度分布を表す。

66t(k+1)

Ti j 6I(k)T (t, 0, y) ≡ 100

Ti−1 j(k)Ti j(k)Ti+1 jT (0, x, y) ≡ T 0(x, y)
図 1: 格子点

T (t, L, y) ≡ 10x -1.1.2境界条件・初期条件

今回の問題では、まず境界条件は正方形の枠が T = 10 の熱浴に接しており、恒常的に端で T = 10 が満たされているとする。また初期条件はたとえば初期の温度分布を T (0, x, y) =T 0(x, y) と与えればよい。今回は時間に関して 1 階の微分方程式なのでこれで初期条件は決まる（初期「速度」の情報が物理的にいらない）。
それでは具体的にコードを見ていこう。このコードは直線の中心に最初 T = 100 の高温スポット (1 メッシュ分) があり、そのほかのメッシュ点は境界と同じ温度で T = 10 であるとしている。また簡単のために κ = 1 として解いている。
コードを理解して動かせたら• 異なる初期条件で計算してみよ。
• 高温スポットの計算で、最初に高温においた点も T = 100 の熱源に接しており、恒常的に T = 100 が保たれるとした場合を調べよ。

1.2

モンテカルロ法による解法

次に前項と同じ熱伝導方程式を、乱数を用いて確率的に解くモンテカルロ法という手法を用いて解いてみよう。

1.2.1

マルコフ過程次に起こる事象の確率が、現在の状態に至るまでの経過とは関係なく、現在の状態によってのみ決定される確率過程のことをマルコフ過程と呼ぶ。この場合、未来の状態が現在の状態のみに依存し過去の経過には依存しない。ちょうど図 1 の場合は時刻 (k + 1) の状態は時刻 (k) の状態のみに依存しており、(k) → (k + 1) の変化を確率的に扱えばマルコフ過程 となるであろう。一次元空間に粒子が存在しその粒子は時間 ∆t の間に ∆x 進むという状況を考える。そ れぞれの粒子のすすむ距離 ∆x について、x から ∆x だけ進み x + ∆x へ移動する確率が P (x, x + ∆x) と与えられるとする。すると、ある時刻の粒子の分布 f (t, x)
は、時間 ∆t 前の分布関数 f (t − ∆t, x) を用いて以下のように書ける:Z +∞f (t, x) =−∞f (t − ∆t, x − ∆x)P (x − ∆x, x)d(∆x)(8)∆t, ∆x がそれぞれ小さいとして、右辺の被積分関数の f (t − ∆t, x − ∆x) をテイラー展開し、微小量の 2 乗までとると、Z +∞ (∂f∂f1 ∂ 2f
1 ∂ 2f2f (t, x) −∆t −∆x +f (t, x) =(∆t)+(∆x)2∂t∂x
2 ∂t22 ∂x2−∞· P (x − ∆x, x)d(∆x))(9)となる。P (x − ∆x, x) は確率なので、
Z +∞

1=−∞

P (x − ∆x, x)d(∆x)

(10)

Z +∞

∆xP (x − ∆x, x)d(∆x)(11)(∆x)2 P (x − ∆x, x)d(∆x)

(12)

∂f
∂f
1 ∂ 2f1 ∂ 2f2∆t −⟨∆x⟩ +(∆t)+⟨(∆x)2 ⟩22∂t∂x
2 ∂t2 ∂x(13)⟨∆x⟩ =⟨(∆x) ⟩ =2−∞

Z +∞−∞に注目し、式 (9) を解くと、

f (t, x) = f (t, x) −となる。ここで、x から ±∆x だけ進む確率が等しい、つまり、 P (x−∆x, x) = P (x+∆x, x)とすると、∆x の平均は ⟨∆x⟩ = 0 となる。これを用いて、上の式をさらに書きなおし、∆t, ∆x についてリーディングオーダーの項のみ残すと、∂f⟨(∆x)2 ⟩ ∂ 2 f=∂t2∆t ∂x2となり、前項で解いた熱伝導方程式と同じ形になる。これは、前項の熱伝導係数 κ =としたものに対応している。

1.2.2(14)⟨(∆x)2 ⟩2∆t

モンテカルロ法による考え方

モンテカルロ法では、なんらかの粒子の集まりが確率的に運動していると考える。ここでは温度の代わりにある単位エネルギーをもった粒子が存在していると考えよう。実際に、N 個の一次元単原子分子からなる温度 T の系の平均エネルギーは、 E = N kB T と与えられるので、エネルギー粒子が移動していると考えてもよいであろう。たとえば、全系を nN kB T個のエネルギー粒子で表すとすると、一つのエネルギー粒子は e =のエネルギーをn持っている。エネルギー粒子は ∆t の間にそれぞれランダムに ∆x だけ移動すると考よう。
このq∆t と ∆x は、熱伝導係数 κ と上の式で関係づけられている。ここで、この粒子が分散σ = ⟨(∆x)2 ⟩ の正規分布で分布する ∆x に従って移動すると考えよう。境界条件、初期条件、熱伝導係数は前項と同じとしている。ここでコードを見て みよう。ここでは、温度 1に対して 10 個の粒子を与えている。
コードを理解して動かせたら• 温度 1 に対して 100 個、1000 個の粒子を与えた場合にどうなるか見てみよ。
• 差分化を用いて計算した結果と比べてみよ。二つの異なる計算を比較する際には同じ時間で比較しなければならないことに注意すること。
• 異なる初期条件で計算してみよ。

1.2.3

モンテカルロ法による二次元熱伝導方程式の解法

次に、モンテカルロ法を用いて

∂T
∂ 2T
∂ 2T
=κ+∂t
∂x2
∂y 2!

(15)

を解くことを考えよう。計算コードは web に置いてあり、中心の 1 メッシュに T = 100、その他のメッシュと境界条件が T = 10 の計算が実行できる。ここでは、温度 1 に対して30 個の粒子を与えている。二次元なので、エネルギー粒子の位置は x, y 成分の座標で表されることに注意しよう。

1.3

提出課題

シュレーディンガー方程式の回の課題を参考に、差分法を用いて二次元熱伝導方程式

∂ 2T
∂ 2T
∂T=κ
+∂t
∂x2
∂y 2!

(16)

を解け。中心の 1 メッシュに T = 100、その他のメッシュと境界条件が T = 10 であるような場合を計算し、モンテカルロ法を用いた計算結果と差分法を用いた結果と比較せよ。また、温度 1 に対して 30 個の粒子を与え、モンテカルロ法を用いて一次元熱伝導方程式を解け。二次元計算は一次元計算と比べてどのくらい時間がかかるか、理由を考え議論せよ。

提出期限: 次回授業開始時