1次元熱伝導問題を有限要素法,pythonで解いてみる

1. モチベーション

この「計算力学（第2版）有限要素法の基礎、日本計算工学会 (編)、森北出版、2012」(以下教科書)のp93、図5.9に有限要素法の要素分割数を変えた場合の計算結果が示されています。
ただし具体的な計算は省略してあるので実際にpythonを用いて計算してみます。numpyのversionは'1.20.3'です。

問題としては1次元熱伝導問題で、教科書の図5.9では要素分割数が2~5の場合のみが図示されていますが、ここでは要素分割数20まで計算した結果を示してみます。
数式については実際の計算に使用した最小限のみ,教科書の数式を分かりやすく変更しています。

2.問題設定

図1のような状態を考えます。方程式と境界条件は以下で与えます。

式(1)-(3)の微分方程式を解くと以下のようになります。

pythonで関数で表現すると以下になります。

import numpy as np

def calc_analytical(x, u0):
    #解析解を求める
    return u0 * np.cosh(1 - x) / np.cosh(1)

3.有限要素法の計算

式(1)~(3)を有限要素法(区分的関数を用いたガラーキン法)で解きます。

例えば分割数3の場合は以下のようになります。ただし$u_0=10$です。

3分割(h=1/3)の場合、有限要素法の方程式は以下のようになります。

4分割(h=1/4)の場合、有限要素法の方程式は以下のようになります。

教科書では頑張ってこの連立方程式を解いていますが、numpyでは逆行列を求める関数が用意されているので、ここではこの関数を使って解きます。

例えば4分割の場合

と置くと,式(6)は

で計算することができます。

5分割以上も同様の考え方で計算できます。これをpythonの関数で表現すると,以下のようになりました。

def element_matrix(h):
    res = np.array([[2 * (h ** 2 + 3), h ** 2 - 6],
                    [h ** 2 - 6 , 2 * (h ** 2 + 3)]])
    return res / (6 * h)

def finite_element(s_mum, u0):
    h = 1 / s_mum
    
    #行列K(n_num*n_num)
    matrix = np.zeros([s_mum + 1, s_mum + 1])
    for i in range(s_mum):
        matrix[i:(i+2),i:(i+2)] += element_matrix(h)
    
    #行列X(s_num*1)を作成
    matrix_u0 = np.zeros([s_mum, 1])
    matrix_u0[0] = -u0 * matrix[1, 0]

    #逆行列を用いて変位行列Uを計算
    res_u = np.dot(np.linalg.inv(matrix[1:, 1:]), matrix_u0)
    return np.append(u0, res_u.flatten())

4. 計算結果

実際に計算してみます。

import matplotlib.pyplot as plt

S_Num = 3 #分割数
U0 = 10. #境界条件

#有限要素法
x_finite = np.linspace(0, 1, S_Num + 1)
res_finite = finite_element(S_Num, U0)

#解析解
x_analytical = np.linspace(0, 1, 100)
res_analytical = calc_analytical(x_analytical, U0)

#x=1(ノイマン境界における相対誤差(%))
print(100 * np.abs(res_finite[-1] - res_analytical[-1])/res_analytical[-1])

#matplotlibで可視化
plt.plot(x_analytical, res_analytical, c="blue", 
        linestyle="dashed", linewidth=2, label="Analytical")
plt.scatter(x_finite, res_finite, c="red")
plt.plot(x_finite, res_finite, c="red", label="Finite Element")
plt.xlabel("Displacement x")
plt.ylabel("Temperature u")
plt.legend()
plt.show()

matplotlibで可視化した結果を示します。

横軸が位置x,縦軸が温度u, 赤が有限要素法, 青が解析解です。
3分割の場合でもソコソコ近似できています。

2分割の場合

20分割の場合

x = 1における温度の相対誤差を示します。

今回計算した結果と教科書の図を比較するとほぼ同じ値を取っているように見えます。
(注p93、図5.9では分割数2-5までです。)