ソースコード

Pythonで有限要素法プログラムの作成

更新日：2020/10/18



ダウンロード

サンプルデータ込み

ソースコード

シンプルな三角形１次元要素の２次元FEM

  1. # -*- coding: utf-8 -*-
  2.  
  3. import math
  4.  
  5. class Material:
  6.     def __init__(self,no):
  7.         self.no = no
  8.         self.young = .0
  9.         self.poisson = .0
  10.         self.thickness = .0
  11.  
  12. class Node:
  13.     def __init__(self,no):
  14.         self.no = no  # 節点番号
  15.         self.x = .0    # X座標値
  16.         self.y = .0    #
  17.         self.disp = [.0,.0] # 強制変位
  18.         self.cond = [0,0] # 拘束条件の有無
  19.         self.force = [.0,0] # 節点荷重
  20.  
  21. class Element:
  22.     def __init__(self,no):
  23.         self.no = no
  24.         self.node = [0,0,0] # 節点番号
  25.         self.x = [.0,.0,.0] # 節点番号のX座標
  26.         self.y = [.0,.0,.0]
  27.         self.u = [.0 for i in range(6)] # 節点番号の変位
  28.         self.elNode = [0,0,0,0,0,0]
  29.         self.matNo = 0
  30.         self.area = .0
  31.         self.Bmat = [[.0 for i in range(6)] for j in range(3)]
  32.         self.Smat = [[.0 for i in range(6)] for j in range(3)]
  33.         self.Dmat = [[.0 for i in range(3)] for j in range(3)]
  34.         self.Kmat = [[.0 for i in range(6)] for j in range(6)]
  35.         self.stress = [.0 for i in range(9)]   # 要素の応力
  36.  
  37.     def setDmatrix(self,mat):
  38.         a = mat.young / (1+mat.poisson) / (1-mat.poisson)
  39.         b = a * mat.poisson
  40.         c = 0.5 * mat.young / (1+mat.poisson)
  41.         self.Dmat = [[a,b,0,],
  42.                      [b,a,0,],
  43.                      [0,0,c]]
  44.  
  45.     def setBmatrix(self):
  46.         self.elNode = [2*self.node[0],2*self.node[0]+1,
  47.                        2*self.node[1],2*self.node[1]+1,
  48.                        2*self.node[2],2*self.node[2]+1]
  49.         Ai = self.y[1] - self.y[2]
  50.         Aj = self.y[2] - self.y[0]
  51.         Ak = self.y[0] - self.y[1]
  52.         Bi = self.x[2] - self.x[1]
  53.         Bj = self.x[0] - self.x[2]
  54.         Bk = self.x[1] - self.x[0]
  55.         dl = Bk * Aj - Ak * Bj
  56.         self.area = 0.5 * dl
  57.         self.Bmat = [[Ai/dl,    0,Aj/dl,    0,Ak/dl,    0],
  58.                      [    0,Bi/dl,    0,Bj/dl,    0,Bk/dl],
  59.                      [Bi/dl,Ai/dl,Bj/dl,Aj/dl,Bk/dl,Ak/dl]]
  60.  
  61.     def setSmatrix(self):
  62.         for i in range(3):
  63.             for j in range(6):
  64.                 s = 0
  65.                 for k in range(3):
  66.                     s += self.Dmat[i][k]*self.Bmat[k][j]
  67.                 self.Smat[i][j] = s
  68.  
  69.     def setKmatrix(self,mat): #要素剛性行列の設定
  70.         v = mat.thickness * self.area
  71.         for i in range(6):
  72.             for j in range(6):
  73.                 s = 0
  74.                 for k in range(3):
  75.                     s += self.Bmat[k][i]*self.Smat[k][j]
  76.                 self.Kmat[i][j] = v * s
  77.  
  78.     def calcStress(self):
  79.         for i in range(3):
  80.             s = 0
  81.             for j in range(6):
  82.                 s += self.Smat[i][j] * self.u[j]
  83.             self.stress[i] = s
  84.         cc = 0.5 * (self.stress[0] + self.stress[1])
  85.         ss = 0.5 * (self.stress[0] - self.stress[1])
  86.         cr = math.sqrt(ss*ss+self.stress[2]*self.stress[2])
  87.         self.stress[3] = cc + cr
  88.         self.stress[4] = cc - cr
  89.         self.stress[5] = cr
  90.         if cr <= ss:
  91.             self.stress[6]  = 0
  92.         else:
  93.             wk = abs((ss+cr)/(ss-cr))
  94.             self.stress[6] = 57.29577951 * math.atan(math.sqrt(wk))
  95.             if self.stress[6] < 0:
  96.                 self.stress[6] = -self.stress[6]
  97.             self.stress[6] = 90 - self.stress[6]
  98.         self.stress[7] = math.sqrt(self.stress[3]*self.stress[3]
  99.                                   - self.stress[3] * self.stress[4]
  100.                                   + self.stress[4] * self.stress[4])
  101.         if self.stress[4] > 0:
  102.             self.stress[8] = self.stress[3]
  103.         elif self.stress[3] > 0:
  104.             self.stress[8] = self.stress[3] - self.stress[4]
  105.         else:
  106.             self.stress[8] = - self.stress[4]
  107.  
  108. class Fem2D():
  109.     def __init__(self):
  110.         self.nodes = []
  111.         self.mats = []
  112.         self.elems = []
  113.         self.numberOfNode = 0
  114.         self.numberOfElement = 0
  115.         self.bandWidth = 0
  116.         self.totalMat = []
  117.  
  118.     def readNodeData(self,fp='input_node.txt'):
  119.         '''節点データの読込'''
  120.         f = open(fp,'r')
  121.         lines = f.readlines()
  122.         f.close()
  123.         lines = lines[1:]
  124.         self.nodes = []
  125.         for line in lines:
  126.             line = line.strip('\n')
  127.             d = line.split('\t')
  128.             d = [x if x is not '' else 0 for x in d] + [0,0,0,0]
  129.             n = Node(int(d[0])-1)
  130.             n.x = float(d[1])
  131.             n.y = float(d[2])
  132.             n.cond = [int(d[3]),int(d[4])]
  133.             n.disp = [float(d[5]), float(d[6])]
  134.             n.force = [float(d[7]),float(d[8])]
  135.             self.nodes.append(n)
  136.  
  137.         self.numberOfNode = len(self.nodes)
  138.         print('Loaded Nodes:%d' % self.numberOfNode)
  139.  
  140.     def readMaterialData(self,fp='input_material.txt'):
  141.         '''材料データの読込'''
  142.         f = open(fp)
  143.         lines = f.readlines()
  144.         lines = lines[1:]
  145.         f.close()
  146.  
  147.         self.mats = []
  148.         for line in lines:
  149.             line = line.strip('\n')
  150.             d = line.split('\t')
  151.             m = Material(int(d[0])-1)
  152.             m.young = float(d[1])
  153.             m.poisson = float(d[2])
  154.             m.thickness = float(d[3])
  155.             self.mats.append(m)
  156.         print('Loaded Materials:%d' % len(self.mats))
  157.  
  158.     def readElementData(self,fp='input_element.txt'):
  159.         '''要素データの読込'''
  160.         f = open(fp)
  161.         lines = f.readlines()
  162.         lines = lines[1:]
  163.         f.close()
  164.  
  165.         self.elems = []
  166.         for line in lines:
  167.             d = line.split('\t')
  168.             e = Element(int(d[0])-1)
  169.             n0 = int(d[1])-1
  170.             n1 = int(d[2])-1
  171.             n2 = int(d[3])-1
  172.             e.node = [n0,n1,n2]
  173.             e.matNo = int(d[4])-1
  174.             e.x[0] = self.nodes[n0].x
  175.             e.x[1] = self.nodes[n1].x
  176.             e.x[2] = self.nodes[n2].x
  177.             e.y[0] = self.nodes[n0].y
  178.             e.y[1] = self.nodes[n1].y
  179.             e.y[2] = self.nodes[n2].y
  180.             self.elems.append(e)
  181.         self.numberOfElement = len(self.elems)
  182.         print('Loaded Elements:%d' % self.numberOfElement)
  183.  
  184.     def setElementMatrix(self):
  185.         '''要素剛性行列'''
  186.         for i in range(self.numberOfElement):
  187.             n = self.elems[i].matNo
  188.             self.elems[i].setDmatrix(self.mats[n])
  189.             self.elems[i].setBmatrix()
  190.             self.elems[i].setSmatrix()
  191.             self.elems[i].setKmatrix(self.mats[n])
  192.  
  193.     def setBandWidth(self):
  194.         n = 0
  195.         for k in range(self.numberOfElement):
  196.             for i in range(6):
  197.                 for j in range(6):
  198.                     ii = self.elems[k].elNode[i]
  199.                     jj = self.elems[k].elNode[j]-ii
  200.                     if jj > n:
  201.                         n = jj
  202.         self.bandWidth = n
  203.         print('BandWidth:%d' % n)
  204.  
  205.     def setTotalMatrix(self):
  206.         '''全体剛性行列'''
  207.         i = self.numberOfNode*2
  208.         j = self.bandWidth+1
  209.         self.totalMat = [[0 for a in range(j)] for b in range(i)]
  210.         for k in range(self.numberOfElement):
  211.             for i in range(6):
  212.                 for j in range(6):
  213.                     ii = self.elems[k].elNode[i]
  214.                     jj = self.elems[k].elNode[j] - ii
  215.                     if jj >= 0:
  216.                         self.totalMat[ii][jj] += self.elems[k].Kmat[i][j]
  217.  
  218.     def setBoundaryCondition(self):
  219.         '''境界条件による全体剛性行列，荷重ベクトルの変更'''
  220.         for i in range(self.numberOfNode):
  221.             node = self.nodes[i]
  222.             for j in range(2):
  223.                 if node.cond[j] != 0:
  224.                     self.totalMat[2*i+j][0] = 1e30
  225.                     node.force[j] = self.totalMat[2*i+j][0]*node.disp[j]
  226.  
  227.     def solve(self):
  228.         '''連立方程式の解（ガウスの消去法）'''
  229.         n1 = self.bandWidth + 1
  230.         n2 = self.numberOfNode*2
  231.         akj = [0 for i in range(n2)]
  232.         aik = [0 for i in range(n2)]
  233.         for k in range(n2):
  234.             m = self.bandWidth + 1
  235.             if k+self.bandWidth+1 > n2:
  236.                 m = n2 - k
  237.             for j in range(1,m):
  238.                 akj[j] = self.totalMat[k][j]/self.totalMat[k][0]
  239.                 aik[j] = self.totalMat[k][j]
  240.             for i in range(1,m):
  241.                 it = i + k
  242.                 im1 = i
  243.                 for j in range(i,m):
  244.                     jt = j - im1
  245.                     self.totalMat[it][jt] -= aik[i]*akj[j]
  246.         for n in range(self.numberOfNode):
  247.             for j in range(2):
  248.                 k = 2*n+j
  249.                 self.nodes[n].force[j] /= self.totalMat[k][0]
  250.                 m = self.bandWidth + 1
  251.                 if (k+self.bandWidth+1)>n2:
  252.                     m = n2 - k
  253.                 for i in range(1,m):
  254.                     nt,jt = divmod(i+k,2)
  255.                     diff = self.totalMat[k][i]*self.nodes[n].force[j]
  256.                     self.nodes[nt].force[jt] -= diff
  257.                 if n == self.numberOfNode:
  258.                     break
  259.         self.nodes[-1].disp[1] = self.nodes[-1].force[1]/self.totalMat[-1][0]
  260.         for n in range(self.numberOfNode-1,-1,-1):
  261.             for j in [1,0]:
  262.                 k = 2*n+j
  263.                 s = 0
  264.                 m = self.bandWidth + 1
  265.                 if (k+self.bandWidth+1)> n2:
  266.                     m = n2-k
  267.                 for i in range(1,m):
  268.                     nt,jt = divmod(i+k,2)
  269.                     s += self.totalMat[k][i]*self.nodes[nt].disp[jt]
  270.                 disp = self.nodes[n].force[j] - s / self.totalMat[k][0]
  271.                 self.nodes[n].disp[j] = disp
  272.  
  273.     def setStress(self):
  274.         '''応力の設定'''
  275.         for i in range(self.numberOfElement):
  276.             e = self.elems[i]
  277.             for j in range(3):
  278.                 e.u[2*j] = self.nodes[e.node[j]].disp[0]
  279.                 e.u[2*j+1] = self.nodes[e.node[j]].disp[1]
  280.             e.calcStress()
  281.  
  282.     def saveTotalMatrix(self,fp='output_matrix.txt'):
  283.         '''全体剛性行列の保存'''
  284.         f = open(fp,'w')
  285.         for i in range(self.numberOfNode*2):
  286.             for j in range(self.bandWidth+1):
  287.                 f.write(str(self.totalMat[i][j])+'\t')
  288.             f.write('\n')
  289.         f.close()
  290.         print('Saved TotalMatrix: %s' % fp)
  291.  
  292.     def saveDisp(self,fp='output_disp.txt'):
  293.         '''節点変位の保存'''
  294.         f = open(fp,'w')
  295.         f.write('節点番号\tX変位\tY変位\n')
  296.         for i in range(self.numberOfNode):
  297.             f.write('%d\t' % (i+1))
  298.             f.write('%f\t' % self.nodes[i].disp[0])
  299.             f.write('%f\n' % self.nodes[i].disp[1])
  300.         f.close()
  301.         print('Saved NodeDisp: %s' % fp)
  302.  
  303.     def saveForce(self,fp='output_force.txt'):
  304.         '''節点力の保存'''
  305.         f = open(fp,'w')
  306.         f.write('節点番号\tX\tY\n')
  307.         for i in range(self.numberOfNode):
  308.             f.write('%d\t' % (i+1))
  309.             f.write('%f\t' % self.nodes[i].force[0])
  310.             f.write('%f\n' % self.nodes[i].force[1])
  311.         f.close()
  312.         print('Saved NodeDisp: %s' % fp)
  313.  
  314.     def saveStress(self,fp='output_stress.txt'):
  315.         '''要素応力の保存'''
  316.         f = open(fp,'w')
  317.         f.write('Element\tX-Stress\tY\tXY\tMAX\tMin\tMax-T\tTheta\tEq-S\tTau_Max\n')
  318.         for i in range(self.numberOfElement):
  319.             f.write(str(i+1)+'\t')
  320.             for j in range(9):
  321.                 f.write(str(self.elems[i].stress[j])+'\t')
  322.             f.write('\n')
  323.         f.close()
  324.         print('Saved Stress: %s' % fp)
  325.  
  326.     def main(self):
  327.         self.readNodeData()
  328.         self.readMaterialData()
  329.         self.readElementData()
  330.         self.setElementMatrix()
  331.         self.setBandWidth()
  332.         self.setTotalMatrix()
  333.         self.setBoundaryCondition()
  334.         self.saveForce()
  335.         self.saveTotalMatrix()
  336.         self.solve()
  337.         self.saveDisp()
  338.         self.setStress()
  339.         self.saveStress()
  340.  
  341. def main():
  342.     fem = Fem2D()
  343.     fem.main()
  344.  
  345. if __name__ == '__main__':
  346.     main()
  347.  
目次にもどる
ページトップへ