有限要素法 事始め2 (1次元全体剛性マトリクスの作製)

#今回の記事について
前回の記事は要素剛性マトリクスの作製プログラムまででしたので、今回は要素剛性マトリクスから全体剛性マトリクスを作製するプログラムの説明をします。前回同様、コロナ社発行のJavaによるはじめての有限要素法を参考としています。

#1次元の有限要素法(全体剛性マトリクスの作製)
##要素剛性マトリクスを結合し全体剛性マトリクスを作製

前回は物体の一つの要素に関して要素剛性マトリクス及び要素剛性方程式を作製しました。これを物体中のそれぞれの要素に対して立式していきます。ある要素iが節点j,kを持つとすれば、要素剛性マトリクス及び要素剛性方程式は、

  {\bf K}^{(i)} \left(
    \begin{array}{c}
      u_j \\
      u_k
    \end{array}
  \right)
  = 
  \left(
    \begin{array}{c}
      K_{11}^{(i)}\ K_{12}^{(i)} \\
      K_{21}^{(i)}\ K_{22}^{(i)}
    \end{array}
  \right)
   \left(
    \begin{array}{c}
      u_j \\
      u_k
    \end{array}
  \right)
  =
  \left(
    \begin{array}{c}
      f^{(i)}_j \\
      f^{(i)}_k
    \end{array}
  \right)
  (1)

と書くことができます。例えば、要素が3つあり、i番目の要素に2つの節点(i, i+1)があるとすれば要素剛性方程式は、以下のようになります。

  \left(
    \begin{array}{c}
      K_{11}^{(1)}\ K_{12}^{(1)} \\
      K_{21}^{(1)}\ K_{22}^{(1)}
    \end{array}
  \right)
   \left(
    \begin{array}{c}
      u_1 \\
      u_2
    \end{array}
  \right)
  =
  \left(
    \begin{array}{c}
      f^{(1)}_1 \\
      f^{(1)}_2
    \end{array}
  \right)
  (2)

  \left(
    \begin{array}{c}
      K_{11}^{(2)}\ K_{12}^{(2)} \\
      K_{21}^{(2)}\ K_{22}^{(2)}
    \end{array}
  \right)
   \left(
    \begin{array}{c}
      u_2 \\
      u_3
    \end{array}
  \right)
  =
  \left(
    \begin{array}{c}
      f^{(2)}_2 \\
      f^{(2)}_3
    \end{array}
  \right)
  (3)

  \left(
    \begin{array}{c}
      K_{11}^{(3)}\ K_{12}^{(3)} \\
      K_{21}^{(3)}\ K_{22}^{(3)}
    \end{array}
  \right)
   \left(
    \begin{array}{c}
      u_3 \\
      u_4
    \end{array}
  \right)
  =
  \left(
    \begin{array}{c}
      f^{(3)}_3 \\
      f^{(3)}_4
    \end{array}
  \right)
  (4)

それぞれの要素が各節点を共有していることを考慮すれば、これらの方程式は一つにまとめることができ、

  \left(
    \begin{array}{cccc}
      K_{11}^{(1)}& K_{12}^{(1)} & 0 & 0 \\
      K_{21}^{(1)}& K_{22}^{(1)} + K_{11}^{(2)} & K_{12}^{(2)} & 0 \\
     0 & K_{21}^{(2)} & K_{22}^{(2)} + K_{11}^{(3)} & K_{12}^{(3)} \\
     0 & 0  & K_{21}^{(3)} & K_{22}^{(3)} \\
    \end{array}
  \right)
   \left(
    \begin{array}{c}
      u_1 \\
      u_2 \\
      u_3 \\
      u_4 
    \end{array}
  \right)
  =
  \left(
    \begin{array}{c}
      f^{(1)}_1\\
      f^{(1)}_2+f^{(2)}_2 \\
      f^{(2)}_3+f^{(3)}_3 \\
      f^{(3)}_4
    \end{array}
  \right)
  (5)

となります。このように要素剛性方程式を一つにまとめたものを、全体剛性方程式といい、左辺の係数マトリクスを全体剛性マトリクスと言います。外力や境界条件を考慮して全体剛性方程式を解くことで全体の応力、ひずみを計算します。

##1次元要素の全体剛性マトリクスの計算プログラム
節点座標、要素の面積・ヤング率、節点にかかる力、各節点の節点変位拘束条件を入力すると、対応した全体剛性方程式を出力するプログラムは以下のようになりました。(節点にかかる力、各節点の節点変位拘束条件は今回は使用しませんが、次回の全体剛性方程式を解く際に必要なパラメータとなります。)プログラム中では、前回作製したElementクラスを使用しています。

import numpy as np

    #データの入力及び計算準備
    #list_x : 節点座標
    #list_area : 各要素の面積
    #list_young : 各要素のヤング率
    #list_force : 各節点にかかる力
    #list_nodeCondition : 各節点の節点変位拘束条件
    def __init__(self,list_x,list_area,list_young,list_force,list_nodeCondition):
        #節点数の設定 NumberOfNode, nnode
        self.NumberOfNode = int(list_x.shape[0])
        nnode = self.NumberOfNode
        
        #要素数の設定 NumberOfElement, nelm
        self.NumberOfElement = nnode-1
        nelm = self.NumberOfElement
        
        #データ格納・計算・出力用のリスト・配列を作製
        self.x = np.zeros(nnode)                     #節点の座標
        self.force = np.zeros(nnode)                 #節点にかかる力
        self.nodeCondition = np.zeros((nnode,2))     #節点変位拘束条件
        self.totalStiffness = np.zeros((nnode,nnode))#全体剛性行列
        self.disp = np.zeros(nnode)                  #節点変位
        self.stress = np.zeros(nelm)                 #応力
        
        #Elementクラス格納用のリスト
        self.elem = np.array([])
        
        #節点座標、節点にかかる力、節点変位条件の初期条件格納
        self.x = list_x
        self.force = list_force
        self.nodeCondition = list_nodeCondition
        
        #Elementクラスを作製し、elemリストに格納
        #その後、要素に関するデータを格納。要素数分実施。
        for i in range(nelm):
            self.elem = np.append(self.elem, Element())
            self.elem[i].young=list_young[i]
            self.elem[i].area=list_area[i]
            for j in range(2):
                #各要素を構成する接点の番号ixと座標xeを格納
                self.elem[i].ix[j]=i+j
                self.elem[i].xe[j]=list_x[i+j]
                
    #全体剛性行列の作製
    def assembleMatrix(self):
        #各要素のBmatrix・要素剛性行列作製
        print("\nAssemble Stiffness Matrix : \n")
        nelm = self.NumberOfElement
        for i in range(nelm):
            self.elem[i].makeBmatrix()
            self.elem[i].makeElementStiffness()
        
        #各要素剛性行列から全体剛性行列を作製
        for k in range(nelm):
            print("  Element",k)
            for i in range(2):
                for j in range(2):
                    ii = int(self.elem[k].ix[i])
                    jj = int(self.elem[k].ix[j])
                    self.totalStiffness[ii,jj] += self.elem[k].Kmatrix[i,j] 

上記の動作例は以下になります。ここでは、節点が4つで各座標座標が(50, 150, 250, 350)(㎜)、各要素の断面積が100㎜$^2$、ヤング率が200 GPaであるとしています。また、今回は使用していませんが、list_nodeConditionで節点座標50㎜での固定とlist_forceで節点座標350㎜に1000 Nの力の印可を指定しています。

入力

#nodeのパラメータ
list_x = np.array([50,150,250,350])
list_force = np.array([0.0,0.0,0.0,1000.])
list_nodeCondition = np.array([[1,0],[0,0],[0,0],[0,0]])

#Elementのパラメータ
list_area = np.array([100,100,100])
list_young = np.array([200000,200000,200000])
 
Fem = Fem1dim(list_x,list_area,list_young,list_force,list_nodeCondition)
Fem.assembleMatrix()
print(Fem.totalStiffness)

出力

Assemble Stiffness Matrix : 

  Element 0
  Element 1
  Element 2
[[ 200000. -200000.       0.       0.]
 [-200000.  400000. -200000.       0.]
 [      0. -200000.  400000. -200000.]
 [      0.       0. -200000.  200000.]]

#まとめ
　今回は要素剛性方程式から全体剛性方程式の立式を説明し、全体剛性マトリクス作製プログラムの紹介を行いました。次回は、変位境界条件の適用と全体剛性方程式(全体剛性マトリクス)の計算するプログラムを紹介します。

#参考書籍

長岐滋 (2010)『Javaによるはじめての有限要素法』 コロナ社