本記事では、有限要素法構造解析の剛性方程式$[K]\{U\}=\{F\}$を解くときなどに使用する、基本的な連立方程式解法のプログラムを紹介します。
有限要素法構造解析における連立方程式
たとえば上図のような2要素のばねモデルについて、節点1、2、3それぞれのつりあい方程式を記述すると、節点変位$(u_1, u_2, u_3)$についての3元1次連立方程式の形になります。（ここでは、ばねの物性により求められる全体剛性マトリクス$[K]$、および節点外力ベクトル$\{F\}$は既知であるとします。）
$$\left\{
\begin{array}{l}
k_1(u_2-u_1)=f_1 \\
-k_1(u_2-u_1)+k_2(u_3-u_2)=f_2 \\
-k_2(u_3-u_2)=f_3
\end{array}
\right.$$
$$\Leftrightarrow
\left\{
\begin{array}{l}
-k_1u_1+k_1u_2=f_1 \\
k_1u_1-(k_1+k_2)u_2+k_2u_3=f_2 \\
k_2u_2-k_2u_3=f_3
\end{array}
\right.$$
これを、マトリクスとベクトルの積の形をした全体剛性方程式に直すと、以下のようになります。
$$\begin{bmatrix}
-k_1 & k_1 & 0 \\
k_1 & -k_1-k_2 & k_2 \\
0 & k_2 & -k_2
\end{bmatrix}
\begin{Bmatrix}
u_1 \\
u_2 \\
u_3
\end{Bmatrix}
\begin{Bmatrix}
f_1 \\
f_2 \\
f_3
\end{Bmatrix}$$
有限要素法構造解析のソルバーでは、この剛性方程式を$\{U\}=\{u_1, u_2, u_3\}^T$について解く、ということを行っています。
解法①：逆行列を掛ける
もっともシンプルな方法は、両辺の左側から、方程式の係数マトリクスの逆行列を掛けるというものです。たとえば前述の剛性方程式について、左側から剛性マトリクス$[K]$を掛けると、以下のようになります。
$$\begin{bmatrix}
-k_1 & k_1 & 0 \\
k_1 & -k_1-k_2 & k_2 \\
0 & k_2 & -k_2
\end{bmatrix}^{-1}
\begin{bmatrix}
-k_1 & k_1 & 0 \\
k_1 & -k_1-k_2 & k_2 \\
0 & k_2 & -k_2
\end{bmatrix}
\begin{Bmatrix}
u_1 \\
u_2 \\
u_3
\end{Bmatrix}$$
$$=\begin{bmatrix}
-k_1 & k_1 & 0 \\
k_1 & -k_1-k_2 & k_2 \\
0 & k_2 & -k_2
\end{bmatrix}^{-1}
\begin{Bmatrix}
f_1 \\
f_2 \\
f_3
\end{Bmatrix}$$
ここで、$[K]^{-1}[K]=1$であるため、
$$\begin{Bmatrix}
u_1 \\
u_2 \\
u_3
\end{Bmatrix}
=
\begin{bmatrix}
-k_1 & k_1 & 0 \\
k_1 & -k_1-k_2 & k_2 \\
0 & k_2 & -k_2
\end{bmatrix}^{-1}
\begin{Bmatrix}
f_1 \\
f_2 \\
f_3
\end{Bmatrix}$$
となり、節点変位ベクトルを算出することができます。Pythonでは、行列Aの逆行列は科学計算ライブラリNumpyのモジュールを使って、以下のコードで算出できます。
A_inv = np.linalg.inv(A)試しに、以下の連立方程式を解かせるPythonプログラムを作成してみます。
$$\begin{bmatrix}
1 & 1 & 2 \\
2 & 1 & -1 \\
1 & -3 & -1
\end{bmatrix}
\begin{Bmatrix}
x \\
y \\
z
\end{Bmatrix}
=
\begin{Bmatrix}
2 \\
-3 \\
8
\end{Bmatrix}$$
import numpy as np
A = ([[ 1, 1, 2],
[ 2, 1, -1],
[ 1, -3, -1]])
b = ([2, -3, 8])
# Getting inverce matrix of A
A_inv = np.linalg.inv(A)
# Multiplying A_inv with b
x = np.dot(A_inv, b)
print(x)このプログラムの実行結果は以下のとおりです。きちんと正解が出力されています。
[ 1. -3. 2.]解法②：Numpyのlinalg.solve関数を使う
ここまで書いておいてなんですが、わざわざ逆行列を求めなくても、Numpyのlinalg.solve関数を用いれば、連立方程式の係数行列をそのまま入力するだけで解を求めることができます。
x = numpy.linalg.solve(A, b)こちらも、解法①で扱った連立方程式を解かせるプログラムを作成します。
import numpy as np
A = ([[ 1, 1, 2],
[ 2, 1, -1],
[ 1, -3, -1]])
b = ([2, -3, 8])
# Solving a simultaneous equation
x = np.linalg.solve(A, b)
print(x)こちらも実行結果は以下の通り。
[ 1. -3. 2.]おまけ：ガウスの消去法
行列方程式の数値解法として最も一般的な手法が、ガウスの消去法といものです。ガウスの消去法の手順については書籍や他サイトで豊富に取り上げられているのでここでは省略しますが、せっかく作ったので以下に載せておきます。なお、このプログラムは書籍「有限要素法のつくり方！（日刊工業新聞社）」の、VBAで作成されたプログラムを書き換えたものです。
import numpy as np
def gaussian_elimination(A, b):
# For equation Ax = b
dof_total = len(b)
A = np.array(A).astype(np.float64)
b = np.array(b).astype(np.float64)
# Forward elimination
for i in range(dof_total):
pivot = A[i, i]
if abs(pivot) < 1e-6:
print('pivot = ', pivot)
raise ZeroDivisionError
A[i] = A[i] / pivot
b[i] = b[i] / pivot
# Computation for line i+1 thru dof_total
for j in range(i+1, dof_total):
p = A[j, i]
A[j] = A[j] - p*A[i]
b[j] = b[j] - p*b[i]
# Backword substitution
x = np.zeros(dof_total)
for i in range(dof_total-1, -1, -1):
x[i] = b[i]
for j in range(i+1, dof_total):
x[i] = x[i] - A[i, j] * x[j]
return x
A = ([[ 1, 1, 2],
[ 2, 1, -1],
[ 1, -3, -1]])
b = ([2, -3, 8])
# Solving a simultaneous equation
# with gaussian elimination
x = gaussian_elimination(A, b)
print(x)こちらも実行結果は以下の通り。前述の、逆行列を掛ける方法、およびsolve関数を使う方法と同一の解が得られました。
[ 1. -3. 2.]