FEMの勉強

プログラミング
Comments

4月に転職、引越と怒涛の2週間が過ぎました。まだカーテンも一部しかついていない状況だけど、ＰＣ開くくらいもういいよね。

結局、転職前の職場もギリギリまで勤め、1週間だけ有給消化週間をもらい、1歳児の相手をしながらだと全然進まない引越準備を妻と大喧嘩しながらも、なんとか引越は無事完了。最初の1週間は家の事に意識を持ってきていたけど、転職者なので、そろそろ戦力として業務することを求められるだろうし、こっちとしてもアピールしていきたいところ。

とりあえず、採用の目論見として、これまでの定常業務＋数値解析を内製化したい意向があるらしい。前の職場ではホントかじった程度なので、新卒のような甘えが効かない分、本腰を入れて勉強しなくてはー。

現在、数値解析として最もベーシックに使われる有限要素法（FEM）の勉強を進めていこうと思い、本を買いました。

 

[商品価格に関しましては、リンクが作成された時点と現時点で情報が変更されている場合がございます。]

ExcelとVBAによる実用有限要素法入門 [ 白井豊 ]

価格：2592円（税込、送料無料) (2017/4/16時点)

 

いかんせん、数学が苦手で（なぜ解析を見込まれたし）、数式を見てもちんぷんかんぷんなので、プログラムで学ぶことにしました。

 

本編の一部とExcel(VBA)はこちらからダウンロードできます。

http://souzousha.iinaa.net/www/FEM/Download.htm

 

とりあえず、この本で出てくる１次元のFEMプログラムを、何も考えずPython化してみました。

 

# -*- coding: utf-8 -*-

#節点番号,座標値,拘束条件,変位,節点力
node_data = [[0, 0,1,0,   0],
             [1, 4,0,0,   0],
             [2, 8,0,0,   0],
             [3,10,0,0,   0],
             [4,12,0,0,2000]]

#要素番号,構成節点0,構成節点2,断面積,ヤング率
elm_data = [[0,0,1,50,400000],
            [1,1,2,50,300000],
            [2,2,3,50,200000],
            [3,3,4,50,100000]]

class Node:
    def __init__(self,no):
        self.no = no  # 節点番号
        self.x = 0    # X座標値
        self.disp = 0 # 強制変位
        self.cond = 0 # 拘束条件の有無
        self.force = 0 # 節点荷重

class Element:
    def __init__(self,no):
        self.no = no
        self.length = .0   # 長さ
        self.area = .0     # 断面積
        self.young = .0    # ヤング率
        self.node = [0,0]  # 節点番号0，1
        self.x = [.0,.0]    # 節点番号0，1のX座標
        self.u = [.0,.0]    # 節点番号0，1の変位
        self.Bmat = [.0,.0] # B行列（要素の形状を示す）
        self.Kmat = [[.0,.0],[.0,.0]] # 要素剛性行列
        self.stress = .0   # 要素の応力

    def setBmatrix(self): #B行列の設定
        self.Bmat[0] = -1/self.length
        self.Bmat[1] =  1/self.length
        
    def setKmatrix(self): #要素剛性行列の設定
        yv = self.young * self.area * self.length
        self.Kmat[0][0] = yv * self.Bmat[0] * self.Bmat[0]
        self.Kmat[0][1] = yv * self.Bmat[0] * self.Bmat[1]
        self.Kmat[1][0] = yv * self.Bmat[1] * self.Bmat[0]
        self.Kmat[1][1] = yv * self.Bmat[1] * self.Bmat[1]
        
class Fem1D():
    
    def setNodeData(self):
        '''節点データ設定'''
        self.numberOfNode = len(node_data)
        self.nodes = []
        for i in range(self.numberOfNode):
            n = Node(i)
            n.x = node_data[i][1]
            n.cond = node_data[i][2]
            n.disp = node_data[i][3]
            n.force = node_data[i][4]
            self.nodes.append(n)
    
    def setElementData(self):
        '''要素データ設定'''
        self.numberOfElement = len(elm_data)
        self.elements = []
        for i in range(self.numberOfElement):
            e = Element(i)
            n0 = elm_data[i][1]
            n1 = elm_data[i][2]
            e.node[0] = n0
            e.node[1] = n1
            e.area = float(elm_data[i][3])
            e.young = float(elm_data[i][4])
            e.length = self.nodes[n1].x - self.nodes[n0].x
            self.elements.append(e)
    
    def setElementStiffness(self):
        '''要素剛性行列設定'''
        for e in self.elements:
            e.setBmatrix()
            e.setKmatrix()
    
    def setTotalMatrix(self):
        '''全体剛性行列設定'''
        #初期設定
        x = y = self.numberOfNode
        self.totalMat = [[0.0 for i in range(x)] for j in range(y)]
        #要素剛性行列の加算
        for k in range(self.numberOfElement):
            for i in range(2):
                for j in range(2):
                    ii = self.elements[k].node[i]
                    jj = self.elements[k].node[j]
                    self.totalMat[ii][jj] += self.elements[k].Kmat[i][j]

    def setBoundaryCondition(self):
        '''境界条件による全体剛性行列，荷重ベクトルの変更'''
        for i in range(self.numberOfNode):
            n1 = self.nodes[i]
            disp = n1.disp
            if n1.cond != 0:
                #非対角要素の設定
                for j in range(self.numberOfNode):
                    n2 = self.nodes[j]
                    n1.force = n2.force - self.totalMat[i][j]*disp
                    self.totalMat[j][i] = 0
                    self.totalMat[i][j] = 0
                #対角要素の設定
                n1.force = disp
                self.totalMat[i][i] = 1
    
    def saveTotalMatrix(self):
        '''全体剛性行列の保存'''
        print('Total Matrix')
        s = ''
        for i in range(self.numberOfNode):
            for j in range(self.numberOfNode):
                s += str(self.totalMat[i][j])+'\t'
            s = s[:-1]+ '\n'
        print(s[:-1])
    
    def solver(self):
        '''連立方程式の解（ガウスの消去法）'''
        n1 = self.numberOfNode - 1
        n2 = n1 - 1
        for k in range(n1+1):
            self.nodes[k].disp = self.nodes[k].force
        #前進消去
        for k in range(n2+1):
            p = self.totalMat[k][k]
            for j in range(k,n1+1):
                self.totalMat[k][j] = self.totalMat[k][j]/p
            self.nodes[k].disp = self.nodes[k].disp/p
            for i in range(k+1,n1+1):
                t = self.totalMat[i][k]
                for j in range(k,n1+1):
                    self.totalMat[i][j] -= t * self.totalMat[k][j]
                self.nodes[i].disp -= t * self.nodes[k].disp
        #後退代入
        self.nodes[n1].disp = self.nodes[n1].disp / self.totalMat[n1][n1]
        for k in range(n1,-1,-1):
            t = self.nodes[k].disp
            for j in range(k+1,n1+1):
                t = t - self.totalMat[k][j] * self.nodes[j].disp
            self.nodes[k].disp = t

    def saveNodeDisp(self):
        '''求まった節点変位データの保存'''
        print('Disp')
        for n in self.nodes:
            print(n.no,n.disp)
    
    def setStress(self):
        '''応力の設定'''
        for i in range(self.numberOfElement):
            e = self.elements[i]
            n1 = self.nodes[e.node[0]]
            n2 = self.nodes[e.node[1]]
            e.u[0] = n1.disp
            e.u[1] = n2.disp
            #cal stress
            s = 0
            for j in [0,1]:
                s += e.young * e.Bmat[j] * e.u[j]
            e.stress = s

    def saveStress(self):
        '''求まった応力データの保存'''
        print('Stress')
        for e in self.elements:
            print(e.no,e.stress)
    
    def main(self):    
        self.setNodeData()
        self.setElementData()
        self.setElementStiffness()
        self.setTotalMatrix()
        self.setBoundaryCondition()
        self.saveTotalMatrix()
        self.solver()
        self.saveNodeDisp()
        self.setStress()
        self.saveStress()

def main():
    fem = Fem1D()
    fem.main()

if __name__ == '__main__':
    main()

 

・・・後でじっくる見るよ。。。

コメントを残す

メールアドレスが公開されることはありません。 ※ が付いている欄は必須項目です