針金内の温度分布の時間変化

こんにちは。おととしです。
針金内の温度分布について、与えられた境界条件から定常状態となった時の温度分布については、以前やってみました。

今度は、針金の両端の温度のみ固定し、徐々に定常状態へ遷移していく様子をシミュレーションしてみたいと思います。
考えるのは、両端をそれぞれ0℃と100℃に固定した針金です。時間変化を考えるには、まず針金をいくつかの細かい部分に分割します。

次にそれぞれの針金の部分について、隣り合ったところからどのように熱が伝わってくるのかを考えます。
点Cへ、点Aおよび点Bから伝わる熱は次の式のようになると考えます。

これは距離あたりの温度の差に対して、温度の伝わる度合いを示す比例定数kをかけたものとなっています。左からは熱が流出し、右から熱が流入すると考えています。
今度はこれを温度変化と対応させます。熱が伝わってきたときの温度の上がり方は、cΔxに比例すると考えられます。
cは比熱に相当する係数です。ここでΔxを掛けるのは、針金の質量は長さに比例するからです。以上のことを踏まえると、Δtの時間が経った後の点Cにおける温度について、次の式が成り立ちます。

変形すると、

と書くことができます。この式をすべての分割に対して行うことで、針金内の温度分布変化を計算することができるのです。
実際にプログラムを実行すると、次のような図が得られます。

# Const Definition
d=21			# split
a=0			# start point
b=1			# end point
x=np.linspace(a,b,d)

itr=500			# iterations
dt=0.01			# Time interval
dx=x[1]			# X interval
k=0.1			# Temperature Rise Rate
c=1			# Heat Capacity
r=k*dt/(c*dx**2)

u=np.zeros((d,itr,))

# Boundary Condition
u[0,0]=0
u[-1,0]=100
t=0

# Matrix Definition
C=np.zeros((d,d,))
C[0,0]=1
C[-1,-1]=1
for i_row in range(1,d-1):
    C[i_row,i_row-1]=r
    C[i_row,i_row]=1-2*r
    C[i_row,i_row+1]=r

# Calculation
for i in range(itr-1):
    u[:,i+1]=np.dot(C,u[:,i])

# Final Plot
fig=plt.figure(figsize=(4,4))
def update(i, x, u, dt):
    if i !=0:
        plt.cla()

    t=i*dt
    plt.plot(x, u[:,i], color = 'orange')
    plt.title('t= %.1f' %t)

# Animation
def animation_play(x, u, dt, itr):
	params={
		'fig':fig,
		'func':update,
		'fargs':(x, u, dt),
		'interval':30,
		'frames':itr,
		'repeat':False,
	}
	ani=anm.FuncAnimation(**params)
	plt.show()

animation_play(x, u, dt, itr)

kやcには本来、物理的な意味があるはずですが、ここではよくわからないまま使っています。
今回行ったシミュレーションは陽解法と呼ばれる種類のものです。t+Δtの様子を、tのときの値のみを使って表そうとしています。これに対し、t+Δtの様子をtのときの値だけでなく、t+Δtのときの値も用いて計算する方法もあるようです。今回の陽解法では、計算条件によってはうまくシミュレーションが行えません。具体的にいうと、途中で1-2rという数がでてきます。この値が負になるとうまくいかなくなります。ですので、rが0.5を超えないようにdt、dx、k、cを設定しなければならないという縛りがあるのです。

時間変化を計算するので、針金内の温度を何度も計算することになります。その回数をitrで指定しています。uという「箱」を作っておいて、計算するたびに結果を書き込んでいきます。Calculationの部分が終わった段階で、itr回の計算はすべて終了しています。

アニメーションについて、図はfigへと引き渡されます。この時にupdateと呼ばれる関数が呼ばれます。関数は30msごとに呼び出され、引数としてx、u、dtを渡します。呼び出しはitr回行われ、itr回の呼出し後の繰り返しは行いません。

プロットについて、animation_playによって関数updateが呼び出されると、引数としてi、x、u、dtが渡されます。iはアニメーションの関数によって指定されているもので、呼び出した回数（0, 1, 2, ..., itr-1)を表しています。updateが呼び出されるたびに一度プロット画面(fig)をリセットし、プロットするということを行います。

あめつちひらく

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

針金内の温度分布の時間変化