二次元の温度分布 - あめつちひらく

以前の記事では、一次元の温度分布をシミュレーションしました。
【よくあるやつ】温度分布シミュレーション - あめつちひらく

次は二次元の温度分布をやりたくなるというものです。
というわけで、自分のできる範囲でやってみました。
参考にしたのは前回と同じく、次の書籍です。

作者:河村哲也
サイエンス社

シミュレーションするのは、4辺を一定の温度に保った平板内の温度分布です。

やることはほとんど変わらないと思っています。二次元のものをどのように表し、計算できるようにするかを考えなければなりません。
まずは、平板を次のように分割することにします。dは分割の数、mはdの2乗です(m=d^2)。pythonで計算することを想定し、スタートを0にして描いています。

Pにおける温度は、周囲のA、B、C、Dの温度の影響を受けて、それら4点の平均の温度となると考えます。式に表すと次のようになります。

これは注目しているそれぞれの点について成り立ちます。ですので、平板内のすべての点において当てはめ、連立方程式をたてることができます。平板の4辺はそれぞれ一定の温度に保たれているので、次のような式ができます。

上の式では上手く表せず、申し訳ないです。今回は平板の温度を考えるのですが、計算の都合上、1列に並びかえています。uの部分が平板を分割したそれぞれの部分の温度を表しています。4辺の温度は一定なので、周囲には影響されず、80℃や40℃、0℃のままとなります。平板の内部は周囲の温度の影響を受けるので、行列に-4が表れています。つまり、上の式の1行1行は分割された小平板1つ1つに対応していて、その小平板が4辺にある時は一定の値となり、4辺でないときは周囲の温度の影響を受けるということを表しているのです。この計算を行うと、次のような図が得られます。

図の左下（0のところ）で0℃となっており、10℃刻みで線が引かれています。一番上が80℃になっていることが確認できます。この図を出力するために書いたのが下のコードとなります。

import numpy as np
import matplotlib.pyplot as plt

# Const. Definition
d=25			# split
m=d**2
u0=np.zeros((m,1,))

a=0			# start point
b=1			# end point
x=np.linspace(a,b,d)
y=np.linspace(a,b,d)
X,Y=np.meshgrid(x,y)

# Boundary Condition
for i in range(m):
    if i<d-1:
        u0[i]=80
    elif i>m-d-1:
        u0[i]=0
    elif i%d==0:
        u0[i]=0
    elif i%d==d-1:
        u0[i]=40

# Matrix Definition
A=np.zeros((m,m,))
for i_row in range(m):
    if i_row<d-1:
        A[i_row,i_row]=1
    elif i_row>m-d-1:
        A[i_row,i_row]=1
    elif i_row%d==0:
        A[i_row,i_row]=1
    elif i_row%d==d-1:
        A[i_row,i_row]=1
    else:
        A[i_row,i_row]=-4
        A[i_row,i_row-1]=1
        A[i_row,i_row+1]=1
        A[i_row,i_row-d]=1
        A[i_row,i_row+d]=1

B=np.linalg.inv(A)

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

# Final Plot
fig=plt.figure(figsize=(3,3))
u_Reshape=np.reshape(u,[d,d])
u_UpsideDown=np.flipud(u_Reshape)
plt.contour(X, Y, u_UpsideDown)
plt.show()

まず、numpyとmatplotlibをインポートします。

次に定数の定義を行います。dが分割数で、1辺をいくつに分割するかを表します。私のPCでは、100を超えると目に見えて計算が遅くなりました。50では余裕でした。今回は平板なので、その2乗の数の計算が行われます。分割数は小さくて良いならなるべく小さい数にしたいところです。平板のx方向とy方向をそれぞれ分割し、後にプロットするためのメッシュをここで定義しておきます。

次に、境界条件を付与します。u0には温度の初期状態を入力します。境界すなわち平板の4辺をそれぞれの条件で表しています。i＜d-1は一番上の辺(80℃)、i＞m-d-1は一番下の辺(0℃)を表します。i%d==0では、dで割り切れる場所、ここでは左の辺(0℃)を表しています。i%d==d-1では、dの倍数より1小さい場所、ここでは右の辺(40℃)を表しています。このとき、頂点の部分を何℃の境界にするかによって微妙にシミュレーション結果が異なるので注意が必要です。

次に行列の定義を行います。初期条件のときと同様、対応する小平板が4辺にあるかどうかで形が異なります。4辺にあるかどうかの条件はu0のときと同様で、そのときA(i_row,i_row)=1となります。4辺にないとき、A(i_row,i_row)=-4で、その上下左右が1となります。上下は列の値に±1をした要素、左右は列の値に±dをした要素が対応します。逆行列をとって左から掛けることで温度分布uが求まります。

最後に、プロットを行います。このままでは直線になっているので、reshape関数で正方形に戻します。また、上下が逆転してしまったので、上下を戻す関数を使っています。この辺はpythonに詳しくないのでもっと良い方法があるかと思います。メッシュの場合は、contour関数を使ってプロットを行います。

2次元になると急に計算量も手間も増えたので大変でした。次は非定常状態から平衡状態となる過程を描写していきたいと思います。