今回の記事について
今回の記事は、コロナ社発行のJavaによるはじめての有限要素法を参考とし、書籍ではJavaで書かれているプログラムをpythonで作ってみたものとなります。また、今回は要素剛性マトリクスを作製するプログラムまでです。
有限要素法の概略
有限要素法とは固体力学の応力やひずみをシミュレーションすることなどに用いられます。固体力学では、物体をたくさんの小さな領域(要素)に分割し、それぞれの要素毎の状態を計算することで、物質全体の変形の様子を求めることになります。物体の状態は3つの微分方程式(応力の平行方程式、ひずみ-変位関係式、応力-ひずみ関係式(構成式))で表現されますが、有限要素法ではそれらの微分方程式を等価な積分表現(仮想仕事の原理)に置き換えることで近似的に変位やひずみを求めます。積分表現により微分方程式は連立一次方程式である剛性方程式に置き換えられます。
また、有限要素法の計算の流れは、物体を各要素に分け、それぞれの要素で剛性方程式(要素剛性マトリクス)を立てたのちに、全ての要素の剛性方程式を統合して全体剛性方程式(全体剛性マトリクス)を組み立てて近似解を計算する、というようなものになります。
基礎微分方程式
簡単のために、断面積$A \ \rm{[L^{2}]}$、長さ$l \ \rm{[L]}$の物体の一方の端点が固定され、もう一方の端点に一次元的で均一な引っ張り荷重$F\ \rm{[F]}$が働く場合を考えていきます。(ここで、$ \rm{[L]} $が長さの次元、$\rm [F]$が力の次元を表すこととします。このような表記を用いれば、例えば、圧力は単位面積当たりにかかる力なので、力の次元を長さの次元の二乗で割ったもの、つまり、$\rm{[FL^{-2}]}$と書き表すことができます。パラメータの次元を明確にすることは、そのパラメータが何を表しているのかと、その他のパラメータとの関連性を分かりやすくしてくれます。)
また、この物質を弾性体(力を加えると変形をし、加えるのを止めると元の形へと戻る物質)とし、そのヤング率が$E \ \rm{[FL^{-2}]}$であるとします。ヤング率とは、弾性体において成り立つ「応力はひずみに比例する」というフックの法則の、ひずみ$\varepsilon$[無次元]にかかる比例係数のことです。(バネに力を加えた時の荷重(あるいは応力)と変位の関係と似ていますが、次元が異なることに気を付けてください。)
これらを仮定すると、この物体における力のつり合いは、物体に生じた応力を$\sigma \ \rm{[FL^{-2}]}$とすれば、以下のように書けます。
$F = \sigma A$ (力のつり合い、応力の平衡方程式) (1’)
さらに、応力を$\sigma$とすると、フックの法則はヤング率$E$とひずみ$\varepsilon$を用いれば以下のように書けます。
$\sigma = E \ \varepsilon$ (材料の性質、応力-ひずみ方程式) (2')
ひずみ$\varepsilon$は、引っ張り荷重下での物体の伸びを$u\rm{[L]}$とすれば以下と表せます。
$\varepsilon = u / L$ (形状の変化、ひずみ-変位関係式) (3')
これらの式(1')-(3')は、基礎微分方程式の特別な場合です。より一般的な状況を考慮すれば、以下の三つの基礎微分方程式が導かれます。
$d(\sigma A)/dx = 0$ (力のつり合い、応力の平衡方程式) (1)
$\sigma = E \ \varepsilon$ (材料の性質、応力-ひずみ方程式) (2)
$\varepsilon = du / dx$ (形状の変化、ひずみ-変位関係式) (3)
弾性体に生じるひずみを解析的に求める場合は、式(1)-(3)の微分方程式を解く必要があるというわけです。
仮想仕事の原理
仮想仕事の原理は以下のように説明されます。詳しい説明は長くなりますので参考文献を参考ください。
物体が変形して釣合いの状態にあるとき, 外力のなした仮想仕事は, 内力のなした仮想仕事に等しい。
仮想仕事の原理を数式で表すと、以下となります。
力が釣り合った状態で, 任意の仮想変位, 仮想ひずみについて
$\int_V \sigma \delta\varepsilon dV = F\delta u$
が成り立つ。
ここで、$\delta\varepsilon$が仮想ひずみ、$\delta u$が仮想変位です。
1次元の場合の有限要素法(要素剛性マトリクスの作製まで)
要素の分割、仮想仕事の原理
さて、有限要素法の準備が整いました。有限要素法では、まず物体を小さな領域である「要素」に分割します。それから、それぞれの要素において仮想仕事の原理を適用し、それぞれの要素の要素剛性方程式(要素剛性マトリクス)を作ります。
一次元の物体を分割した時の、ある要素一つについて仮想仕事の原理を考えていきましょう。要素は節点(要素と要素の繋ぎ目)を二つ持っているので、一方を節点1(座標$x_1$)、もう一方を節点2(座標$x_2$)とし、それぞれに働く内力を$f_1$、$f_2$、それぞれ点の変位を$u_1$、$u_2$とおきます。ここの応力が$\sigma$、各節点の仮想変位を$\delta u_1$、$\delta u_2$とすれば、この要素の仮想仕事の原理は、
$\int_{x_1}^{x_2} \sigma \delta\varepsilon A dx = f_1 \cdot \delta u_1 + f_2 \cdot \delta u_2 $ (4)
と書けます。
変位の線形近似
ここで各要素の内部で、変位$u$が座標$x$の1次関数で変化していると仮定します。
$u(x) = ax+b$ (5)
この変位は節点1および2では$u_1$、$u_2$と等しくなければならないので、
$u(x_1) = ax_1+b =u_1$ (6)
$u(x_1) = ax_2+b =u_2$ (7)
これを$a$と$b$に対してまとめると、
$a=\frac{u_2 - u_1}{x_2 - x_1} = \frac{u_2 - u_1}{l}$ (8)
$b=\frac{x_2 u_1 - x_1 u_2}{x_2 - x_1}=\frac{x_2 u_1 - x_1 u_2}{l}$ (9)
となります。ここで、$l=x_2-x_1$は要素の長さです。(8)と(9)を用いれば、(5)は、
$u(x) = ax+b = \frac{u_2 - u_1}{l}x+\frac{x_2 u_1 - x_1 u_2}{l}$ (10)
さらに、(10)式を書き換えると、
$u(x) = \frac{x_2 - x}{l}u_1+\frac{x - x_1}{l} u_2 =N_1(x)u_1 + N_2(x)u_2$ (11)
となり、節点変位$u_1$、$u_2$を未知数とした式に書き換えることができました。また、$N_1(x)$、$N_2(x)$は形状関数と言います。また、式(11)は行列表記することができ、
$u(x) = {\bf N}^{{\rm T}} {\bf U} $ (12)
と書けます。ここで、${\bf N}$と${\bf U}$は、
{\bf N} = \left(
\begin{array}{c}
N_1(x) \\
N_2(x)
\end{array}
\right),
{\bf U} = \left(
\begin{array}{c}
u_1 \\
u_2
\end{array}
\right) (13)
というベクトルです。
ひずみと応力の近似的な表示
変位の式(11)と、ひずみ-変位関係式(3)を用いれば、ひずみ$\varepsilon$は、
\varepsilon = \frac{d}{dx}(N_1(x)u_1 + N_2(x)u_2) = \frac{dN_1(x)}{dx}u_1 + \frac{dN_2(x)}{dx}u_2 (14)
と書けます。これもまた、行列表記を行えば、
$\varepsilon = {\bf B}^{{\rm T}} {\bf U} $ (15)
と書き換えることができます。ここで、${\bf B}$は、${\bf B}$マトリクス、あるいはひずみ-変位マトリクスと呼ばれ、以下のように書き表せます。
{\bf B} = \left(
\begin{array}{c}
\frac{dN_1(x)}{dx} \\
\frac{dN_2(x)}{dx}
\end{array}
\right)
=
\left(
\begin{array}{c}
\frac{-1}{l} \\
\frac{1}{l}
\end{array}
\right)
(16)
また、応力$\sigma$は応力-ひずみ関係式(2)と式(15)から、
$\sigma = E \varepsilon = E{\bf B}^{{\rm T}} {\bf U} $ (17)
と表せます。
仮想的変位とひずみ
式(12)と式(15)から、仮想変位$\delta u$と仮想ひずみ$\delta \varepsilon$は、行列表記によって、
$\delta u = {\bf N}^{{\rm T}} \delta{\bf U} $ (18)
$\delta \varepsilon = {\bf B}^{{\rm T}} \delta{\bf U} $ (19)
と書き表すことができます。
要素剛性方程式(要素剛性マトリクス)
ここで、仮想仕事の原理の式(4)について考えていきます。
$\int_{x_1}^{x_2} \sigma \delta\varepsilon A dx = f_1 \cdot \delta u_1 + f_2 \cdot \delta u_2 $ (4)
まずは、式(4)の左辺の積分内部を行列表記で書き換えると、
$\sigma \delta \varepsilon = E{\bf B}^{{\rm T}} {\bf U}{\bf B}^{{\rm T}} \delta{\bf U} = E {\bf U}^{{\rm T}}{\bf B}{\bf B}^{{\rm T}} \delta{\bf U}$ (20)
となるので、これを式(4)の左辺に適用すれば、
$\int_{x_1}^{x_2} \sigma \delta\varepsilon A dx = \int_{x_1}^{x_2} E {\bf U}^{{\rm T}}{\bf B}{\bf B}^{{\rm T}} \delta{\bf U}A dx $ (21)
(21)式で要素の中で断面積$A$、ヤング率$E$が一定であると考えれば、
$\int_{x_1}^{x_2} \sigma \delta\varepsilon A dx = E {\bf U}^{{\rm T}}{\bf B}{\bf B}^{{\rm T}} \delta{\bf U}A \int_{x_1}^{x_2} dx = E V {\bf U}^{{\rm T}}{\bf B}{\bf B}^{{\rm T}} \delta{\bf U} $ (22)
ここで、$V$は、$A \int_{x_1}^{x_2} dx = El = V$であり、要素の体積を表しています。また、${\bf K} = E V {\bf B}{\bf B}^{{\rm T}}$とおけば、以下のように書き換えることができます。
$\int_{x_1}^{x_2} \sigma \delta\varepsilon A dx = {\bf U}^{{\rm T}} {\bf K} \delta{\bf U} $ (23)
${\bf K}$は要素剛性マトリクスと呼ばれます。
次に式(4)の右辺を行列表記すれば、
$f_1 \cdot \delta u_1 + f_2 \cdot \delta u_2 = {\bf F}^{{\rm T}} \delta{\bf U}$ (24)
ここで、${\bf F}$は、以下です。
{\bf F} = \left(
\begin{array}{c}
f_1 \\
f_2
\end{array}
\right)
(25)
こうして、仮想仕事の原理の式(4)は、式(23)と式(24)を用いれば、以下のように書き換えられます。
${\bf U}^{{\rm T}} {\bf K} \delta{\bf U} = {\bf F}^{{\rm T}} \delta{\bf U}$ (26)
また、上式を$\delta{\bf U}$でまとめれば、
$({\bf U}^{{\rm T}} {\bf K} - {\bf F}^{{\rm T}} )\delta{\bf U} = {\bf 0}$ (27)
ここで、仮想仕事の原理は任意の仮想変位$\delta{\bf U}$で成り立つので、
${\bf U}^{{\rm T}} {\bf K} = {\bf F}^{{\rm T}}$ (28)
この転置をとれば、
$ {\bf K}{\bf U} = {\bf F}$ (29)
となります。式(29)は要素剛性方程式と呼ばれる方程式です。初めにも述べたように、この要素剛性方程式を各要素に対して立式したのち、全ての要素の剛性方程式を統合して全体剛性方程式(全体剛性マトリクス)を組み立てて近似解を計算することとなります。
1次元要素の要素剛性方程式(要素剛性マトリクス)の計算プログラム
さて、pythonで1次元の要素剛性方程式を計算するプログラムを作製してみました。使用したライブラリはnumpyで、コードは以下です。与えられた条件から、ある一つの要素に対して、${\bf B}$マトリクスおよび要素剛性マトリクス${\bf K}$、応力を計算することができます。
import numpy as np
class Element:
#接点の数
NumberOfNodeInElement = 2
def __init__(self):
#要素の長さ
self.length = 0
#要素の断面積
self.area = 0
#ヤング率
self.young = 0
#要素を構成する接点番号
self.ix = np.zeros(self.NumberOfNodeInElement)
#要素を構成する接点の座標値
self.xe = np.zeros(self.NumberOfNodeInElement)
#要素の接点変位
self.ue = np.zeros(self.NumberOfNodeInElement)
#Bmatrix
self.Bmatrix = np.zeros(self.NumberOfNodeInElement)
#要素剛性行列
self.Kmatrix = np.zeros((self.NumberOfNodeInElement,self.NumberOfNodeInElement))
#Bmatrixの計算
def makeBmatrix(self):
self.length = self.xe[1]-self.xe[0]
self.Bmatrix = np.array([[-1/self.length, 1/self.length]])
#要素剛性行列の計算
def makeElementStiffness(self):
self.volume = self.area*self.length
self.Kmatrix = self.young*self.volume*np.dot(self.Bmatrix.T,self.Bmatrix)
#応力の計算
def stressCalculation(self):
return self.young*np.vdot(self.Bmatrix,self.ue)
例えば、節点1・節点2の座標が$x_1 = 50$(㎜)、$x_2 = 150$(㎜)、断面積100$㎜^2$、ヤング率200 GPaである要素に対して計算を行うと
Ele = Element()
Ele.xe = np.array([50,150])
Ele.area = 100
Ele.young = 200000
Ele.ue = np.array([0.01,0.025])
Ele.makeBmatrix()
Ele.makeElementStiffness()
scal = Ele.stressCalculation()
であり、この場合の計算結果は以下のようになります。
#Bmatrix
[[-0.01 0.01]]
#Bmatrix
[[ 200000. -200000.]
[-200000. 200000.]]
#応力
30.000000000000004
まとめ
今回は有限要素法の1次元要素剛性マトリクスまでの説明と要素剛性マトリクス作製プログラムの紹介を行いました。次回は全体剛性方程式(全体剛性マトリクス)の組み立てと計算を実施します。
参考書籍
長岐滋 (2010)『Javaによるはじめての有限要素法』 コロナ社