あめつちひらく

物理が好きです。数値シミュレーション始めました。資格の勉強とか。

【よくあるやつ】温度分布シミュレーション


こんにちは。おととしです。数値シミュレーションといえばこれかな、というのがあるのでやってみました。それは、熱が伝わっていくのをシミュレーションするものです。今回、数値シミュレーションというものがどういうものかよくわかっていないので、本を1冊購入しました。この本をもとにやっていきたいと思います。コードもそれっぽいものを考えて書いてきます。ただ、この本には考え方が書いてあるのみで、実際のコードがどのようなものかは書いてありません。ですので、この記事で紹介するコードは正解であるとは限りません。

最初にやるのは、針金(1次元)の温度分布です。一端を0℃、もう一方の端を100℃に保って放置すると、最終的にどのような温度分布になるのかというものです。

針金の両端の温度を固定する

数値シミュレーションでは、針金をいくつかに区切り、それぞれの点においてどのような温度となるかを計算していきます。上の図でいうと、点Cの温度は隣り合う点Aおよび点Bの温度の影響を受けて決定されます。これをプログラムしていきます。

用いる式は次の式です。点Cの温度は、点Aと点Bの温度の平均であるという式です。

これは注目している1つの点に対して成り立ちます。これを針金内のすべての点について当てはめると、連立方程式ができます。また、針金の両端は0℃と100℃に固定されています。これらのことをまとめると、次のような行列を用いた式ができます。

あとは、逆行列を求めて、両辺に左から掛けると終了です。

import numpy as np
import matplotlib.pyplot as plt
# Const Difinition
num=11      # split
a=0         # start point
b=10        # end point
x=np.linspace(a,b,num)

u0=np.zeros((num,1,))

# Initial state
u0[-1]=100

# Matrix Definition
A=np.zeros((num,num,))
for i_row in range(1,num-1):
    A[i_row,i_row-1]=1
    A[i_row,i_row]=-2
    A[i_row,i_row+1]=1
A[0,0]=1
A[num-1,num-1]=1

B=np.linalg.inv(A)

# Solve
u=np.dot(B,u0)

fig=plt.figure(figsize=(3,3))
plt.plot(x, u)
plt.show()

この結果は次のように出てきます。

定常状態における針金の温度分布

温度は針金内に線形に分布していることがわかります。