ソースコード

Pythonで個別要素法プログラムの作成



ソースコード

早速ですが、全体のソースコードをば。
動作確認：Python3.8 

  1. # -*- coding: utf-8 -*-
  2.  
  3. import math
  4. import time
  5. import tkinter
  6.  
  7. G = 9.80665 #重力加速度
  8.  
  9. class Element:
  10.     def __init__(self):
  11.         self.n = 0 #要素No.
  12.         self.r = 0.0 #半径
  13.         self.x = 0.0 #X座標
  14.         self.y = 0.0 #Y座標
  15.         self.a = 0.0 #角度
  16.         self.dx = 0.0 #X方向増加量
  17.         self.dy = 0.0 #Y方向増加量
  18.         self.da = 0.0 #角度増加量
  19.         self.vx = 0.0 #X方向速度
  20.         self.vy = 0.0 #Y方向速度
  21.         self.va = 0.0 #角速度
  22.         self.fx = 0.0 #Ｘ方向節点力
  23.         self.fy = 0.0 #Ｙ方向節点力
  24.         self.fm = 0.0 #モーメント
  25.         self.en = [] #弾性力（直方向）
  26.         self.es = [] #弾性力（せん断方向）
  27.  
  28. class Particle(Element):
  29.     def __init__(self,x,y,r=5.0):
  30.         super().__init__()
  31.         self.etype = 1 #要素タイプ
  32.         self.r = r #半径
  33.         self.x = x #X座標
  34.         self.y = y #Y座標
  35.         #粒子専用
  36.         self.m = 0.0 #質量
  37.         self.Ir = 0.0 #慣性モーメント
  38.  
  39. class Line(Element):
  40.     def __init__(self,x1=0,y1=0,x2=0,y2=0):
  41.         super().__init__()
  42.         self.etype = 2 #要素タイプ
  43.  
  44.         #線要素専用
  45.         self.x1 = x1
  46.         self.y1 = y1
  47.         self.x2 = x2
  48.         self.y2 = y2
  49.  
  50. class Interface:
  51.     def __init__(self,kn=0,ks=0,etan=0,etas=0,frc=0):
  52.         self.kn = kn     #弾性係数（法線方向）
  53.         self.etan = etan #粘性係数（法線方向）
  54.         self.ks = ks     #弾性係数（せん断方向）
  55.         self.etas = etas #粘性係数（せん断方向）
  56.         self.frc = frc   #摩擦係数
  57.  
  58. class InterfaceContainer:
  59.     def __init__(self):
  60.         self.data = {}
  61.  
  62.     def set(self,interface,etype1,etype2):
  63.         if not etype1 in self.data:
  64.             self.data[etype1] = {}
  65.         if not etype2 in self.data:
  66.             self.data[etype2] = {}
  67.         self.data[etype1][etype2] = interface
  68.         self.data[etype2][etype1] = interface
  69.  
  70.     def get(self,etype1,etype2):
  71.         return self.data[etype1][etype2]
  72.  
  73. class DEM:
  74.  
  75.     def __init__(self):
  76.         self.particles = []
  77.         self.lines = []
  78.         self.interfaces = InterfaceContainer()
  79.  
  80.         self.rho = 10 #粒子間密度
  81.         self.dt = 0.001 #⊿ｔ
  82.         self.step = 0 #ステップ
  83.  
  84.     def setup(self):
  85.         #境界（暫定処理）
  86.         l1 = Line(self.area[0],self.area[2],self.area[0],self.area[3])
  87.         l2 = Line(self.area[1],self.area[2],self.area[1],self.area[3])
  88.         l3 = Line(self.area[0],self.area[2],self.area[1],self.area[2])
  89.         l4 = Line(self.area[0],self.area[3],self.area[1],self.area[3])
  90.         self.lines += [l1,l2,l3,l4]
  91.  
  92.         self.par_count = len(self.particles)
  93.         self.line_count = len(self.lines)
  94.         self.elem_count = self.par_count + self.line_count
  95.  
  96.         #粒子要素
  97.         pe = self.particles
  98.         for i in range(self.par_count):
  99.             pe[i].n = i
  100.             pe[i].m = 4.0/3.0*math.pi*self.rho*pe[i].r**3.0 # 質量
  101.             pe[i].Ir = math.pi*self.rho*pe[i].r**4.0/2.0 #慣性モーメント
  102.             pe[i].en = [0.0 for j in range(self.elem_count)]
  103.             pe[i].es = [0.0 for j in range(self.elem_count)]
  104.  
  105.         #線要素
  106.         for i in range(self.line_count):
  107.             self.lines[i].n = self.par_count + i
  108.  
  109.     def cycle(self,step=1):
  110.         for i in range(step):
  111.             self.next_step()
  112.  
  113.     def next_step(self):
  114.         self.calc_force()
  115.         self.update_coord()
  116.         self.step += 1
  117.  
  118.     def calc_force(self):
  119.         self._reset_force()
  120.         self._particles_collision()
  121.         self._par_and_line_collision()
  122.         self._apply_gravitiy()
  123.  
  124.     def _reset_force(self):
  125.         for p in self.particles:
  126.             p.fx = 0
  127.             p.fy = 0
  128.             p.fm = 0
  129.  
  130.     def _particles_collision(self):
  131.         for i in range(self.par_count):
  132.             for j in range(i+1,self.par_count):
  133.                 p1 = self.particles[i]
  134.                 p2 = self.particles[j]
  135.                 lx = p1.x - p2.x
  136.                 ly = p1.y - p2.y
  137.                 ld = (lx**2+ly**2)**0.5
  138.                 if (p1.r+p2.r)>ld: #接触
  139.                     cos_a = lx/ld
  140.                     sin_a = ly/ld
  141.                     self._force_par2par(p1,p2,cos_a,sin_a)
  142.                 else:
  143.                     p1.en[p2.n] = 0.0
  144.                     p1.es[p2.n] = 0.0
  145.  
  146.     def _par_and_line_collision(self):
  147.         #粒子と線の接触判定
  148.         for i in range(self.par_count):
  149.             for j in range(self.line_count):
  150.                 hit = False
  151.                 p = self.particles[i]
  152.                 l = self.lines[j]
  153.                 th0 = math.atan2(l.y2-l.y1, l.x2-l.x1)
  154.                 th1 = math.atan2(p.y-l.y1, p.x-l.x1)
  155.                 a = math.sqrt((p.x-l.x1)**2+(p.y-l.y1)**2)
  156.                 d = abs(a*math.sin(th1-th0))
  157.                 if d < p.r:
  158.                     b = math.sqrt((p.x-l.x2)**2+(p.y-l.y2)**2)
  159.                     s = math.sqrt((l.x2-l.x1)**2+(l.y2-l.y1)**2)
  160.                     s = math.sqrt(s**2 + p.r**2)
  161.                     if a < s and b < s:
  162.                         c = math.sqrt(a**2-d**2)
  163.                         x = l.x1 + c*math.cos(th0)
  164.                         y = l.y1 + c*math.sin(th0)
  165.                         hit = True
  166.                     elif a < b and a < p.r:
  167.                         x = l.x1
  168.                         y = l.y1
  169.                         hit = True
  170.                     elif b < p.r:
  171.                         x = l.x2
  172.                         y = l.y2
  173.                         hit = True
  174.                 if hit:
  175.                     lx = x - p.x
  176.                     ly = y - p.y
  177.                     ld = math.sqrt(lx**2+ly**2)
  178.                     cos_a = lx/ld
  179.                     sin_a = ly/ld
  180.                     self._force_par2par(p,l,cos_a,sin_a)
  181.                 else:
  182.                     p.en[l.n] = 0.0
  183.                     p.es[l.n] = 0.0
  184.  
  185.     def _apply_gravitiy(self):
  186.         for p in self.particles:
  187.             p.fy += -G*p.m #重力
  188.  
  189.     def _force_par2par(self,p1,p2,cos_a,sin_a):
  190.  
  191.         #相対的変位増分
  192.         un = +(p1.dx-p2.dx)*cos_a+(p1.dy-p2.dy)*sin_a
  193.         us = -(p1.dx-p2.dx)*sin_a+(p1.dy-p2.dy)*cos_a+(p1.r*p1.da+p2.r*p2.da)
  194.  
  195.         inf = self.interfaces.get(p1.etype,p2.etype)
  196.  
  197.         #合力
  198.         p1.en[p2.n] += inf.kn*un
  199.         p1.es[p2.n] += inf.ks*us
  200.         hn = p1.en[p2.n] + inf.etan*un/self.dt
  201.         hs = p1.es[p2.n] + inf.etas*us/self.dt
  202.  
  203.         if hn <= 0.0:
  204.             #法線力がなければ、せん断力は０
  205.             hs = 0.0
  206.         elif abs(hs) > inf.frc*hn:
  207.             #摩擦力以上のせん断力は働かない
  208.             hs = inf.frc*abs(hn)*hs/abs(hs)
  209.         #合力
  210.         p1.fx += -hn*cos_a + hs*sin_a
  211.         p1.fy += -hn*sin_a - hs*cos_a
  212.         p1.fm -=  p1.r*hs
  213.         p2.fx +=  hn*cos_a - hs*sin_a
  214.         p2.fy +=  hn*sin_a + hs*cos_a
  215.         p2.fm -=  p2.r*hs
  216.  
  217.     def update_coord(self):
  218.         pe = self.particles
  219.         for i in range(self.par_count):
  220.             #位置更新（オイラー差分）
  221.             ax = pe[i].fx/pe[i].m
  222.             ay = pe[i].fy/pe[i].m
  223.             aa = pe[i].fm/pe[i].Ir
  224.             pe[i].vx += ax*self.dt
  225.             pe[i].vy += ay*self.dt
  226.             pe[i].va += aa*self.dt
  227.             pe[i].dx = pe[i].vx*self.dt
  228.             pe[i].dy = pe[i].vy*self.dt
  229.             pe[i].da = pe[i].va*self.dt
  230.             pe[i].x += pe[i].dx
  231.             pe[i].y += pe[i].dy
  232.             pe[i].a += pe[i].da
  233.  
  234.  
  235. class Window(tkinter.Tk):
  236.     dem = DEM
  237.  
  238.     def __init__(self,dem):
  239.         super().__init__()
  240.         self.dem = dem
  241.         self.loop_n = 0
  242.  
  243.         self.scale = 1.0
  244.         dem_width = self.dem.area[1] - self.dem.area[0]
  245.         dem_height= self.dem.area[3] - self.dem.area[2]
  246.         self.width = dem_width + 10
  247.         self.height = dem_height + 10
  248.         self.offset_x = (dem_width - self.width)/2
  249.         self.offset_y = (dem_height - self.height)/2
  250.  
  251.         self.canvas = tkinter.Canvas(self, bg="white")
  252.         self.canvas.pack(fill=tkinter.BOTH,expand=True)
  253.         self.geometry('%dx%d' % (self.width,self.height))
  254.         self.title('DEM')
  255.  
  256.         for l in self.dem.lines:
  257.             xy = self.view_pos([l.x1,l.y1,l.x2,l.y2])
  258.             self.canvas.create_line(xy,width=1)
  259.  
  260.         self.redraw()
  261.         self.update_idletasks()
  262.  
  263.     def calcloop(self):
  264.         if self.loop_n % 5 == 0:
  265.             print('Step %d' % self.dem.step)
  266.         if self.loop_n % 1 == 0:
  267.             self.redraw()
  268.             time.sleep(0.001) #ゆっくり見たい場合
  269.         if self.loop_n >= 1000:
  270.             print('解析終了.設定最大ループに達しました')
  271.         else:
  272.             self.after(0,self.calcloop)
  273.         self.dem.cycle(10)
  274.         self.loop_n += 1
  275.         self.update_idletasks()
  276.  
  277.     def redraw(self):
  278.         self.canvas.delete('par')
  279.         h = self.height
  280.         for p in self.dem.particles:
  281.             x1,y1 = self.view_pos([p.x-p.r,p.y-p.r])
  282.             x2,y2 = self.view_pos([p.x+p.r,p.y+p.r])
  283.             self.canvas.create_oval(x1,y1,x2,y2,tags='par')
  284.  
  285.             x1,y1 = self.view_pos([p.x,p.y])
  286.             x2,y2 = self.view_pos([p.x+p.r*math.cos(p.a),
  287.                                     p.y+p.r*math.sin(p.a)])
  288.             self.canvas.create_line(x1,y1,x2,y2,tags='par')
  289.  
  290.     def view_pos(self,coords,offset=(0,0)):
  291.         s = self.scale # 表示倍率
  292.         h = self.height #表示画面高さ
  293.         w = self.width #表示画面幅
  294.         xy_list = []
  295.         for i in range(0,len(coords),2):
  296.             x = round(s*coords[i]) - self.offset_x
  297.             y = round(h-s*coords[i+1]) + self.offset_y
  298.             x = x + offset[0]
  299.             y = y + offset[1]
  300.             xy_list.append(x)
  301.             xy_list.append(y)
  302.         return xy_list
  303.  
  304. def main():
  305.     print('初期設定中...',end='')
  306.     dem = DEM()
  307.     dem.area = [0,300,0,200]
  308.     dem.dt = 0.01
  309.  
  310.     p1 = Particle(250,190,10)
  311.     p2 = Particle(190,190,10)
  312.     p3 = Particle(120,190,10)
  313.     p4 = Particle(60,190,10)
  314.     dem.particles = [p1,p2,p3,p4]
  315.  
  316.     l1 = Line(50,140,300,170)
  317.     l2 = Line(0,110,170,80)
  318.     l3 = Line(30,30,300,50)
  319.     dem.lines = [l1,l2,l3]
  320.  
  321.     #粒子同士
  322.     inf1 = Interface()
  323.     inf1.kn  = 100000 #弾性係数（法線方向）
  324.     inf1.etan= 50000 #粘性係数（法線方向）
  325.     inf1.ks  = 5000 #弾性係数（せん断方向）
  326.     inf1.etas= 5000 #粘性係数（せん断方向）
  327.     inf1.frc = 100 #摩擦係数
  328.     #粒子と線要素
  329.     inf2 = Interface()
  330.     inf2.kn  = 1000000
  331.     inf2.etan= 50000
  332.     inf2.ks  = 5000
  333.     inf2.etas= 5000
  334.     inf2.frc = 100
  335.  
  336.     dem.interfaces.set(inf1,1,1)
  337.     dem.interfaces.set(inf2,1,2)
  338.  
  339.     dem.setup()
  340.  
  341.     print('完了')
  342.     print('粒子要素数： %d ' % dem.par_count)
  343.     print('解析開始')
  344.  
  345.     w = Window(dem)
  346.     w.after(0,w.calcloop)
  347.     w.mainloop()
  348.  
  349. if __name__ == '__main__':
  350.     main()

実行するとウィドウが立ち上がり、 ４つの球がコロコロ転がる描写が実行させるはずです。

目次にもどる
ページトップへ