FEMやったりDEMやったり節操がないですが、手っ取り早くできそうなことには手をつけていくスタンスでさ。
個別要素法に関して、コードを公開してくれている講義資料を見つけましたので、とりあえずPython化。ありがたい。
大学の講義資料なので、一部、意図的に書かれていないこととかあって、あってることやら。
http://www.fml.t.u-tokyo.ac.jp/lecture/lecture_6/text_DEM2014.pdf
粒子間の動作だけだけど、とりあえず動いた。
おー
# -*- coding: utf-8 -*-
import sys
import math
from PySide import QtGui,QtCore
G = 9.80665 # 重力加速度
rho = 10 # 粒子間密度
pn = 0 #粒子数
class Particle():
def __init__(self,pos):
global pn
self.n = pn
pn += 1
self.r = 5 #半径
self.x = pos[0] #X座標
self.y = pos[1] #Y座標
self.a = 0 #角度
self.dx = 0 #X方向増加量
self.dy = 0 #Y方向増加量
self.da = 0 #角度増加量
self.vx = 0 #X方向速度
self.vy = 0 #Y方向速度
self.va = 0 #角速度
self.fy = 0
self.fx = 0
self.fm = 0
self.m = 4.0/3.0*math.pi*rho*self.r**3 # 質量
self.Ir = math.pi*rho*self.r**4/2.0 #慣性モーメント
self.en = [] #弾性力（直方向）
self.es = [] #弾性力（せん断方向）
def config(self):
self.en = [0 for i in range(pn)]
self.es = [0 for i in range(pn)]
def nextStep(self,dt):
#位置更新（オイラー差分）
ax = self.fx/self.m
ay = self.fy/self.m
aa = self.fm/self.Ir
self.vx += ax*dt
self.vy += ay*dt
self.va += aa*dt
self.dx = self.vx*dt
self.dy = self.vy*dt
self.da = self.va*dt
self.x += self.dx
self.y += self.dy
self.a += self.da
class DEM():
def __init__(self):
self.step = 0
self.dt = 0.005
p1 = Particle([-100,0])
p1.vy = 30
p1.vx = 60
p2 = Particle([100,0])
p2.vy = 30
p2.vx = -60
p3 = Particle([10,-20])
p3.vy = 30
p3.vx = -20
p4 = Particle([-80,60])
p4.vx = 20
p4.vy = 5
self.pe = []
self.pe.append(p1)
self.pe.append(p2)
self.pe.append(p3)
self.pe.append(p4)
for p in self.pe:
p.config()
def calcForce(self):
#2粒子間の力
combs = [(p1, p2) for p1 in self.pe for p2 in self.pe]
for p1,p2 in combs:
if p1.n == p2.n:
continue
#接触判定
lx = p1.x - p2.x
ly = p1.y - p2.y
ld = math.sqrt(lx**2+ly**2)
if (p1.r+p2.r)>ld:
#接触
cos_a = lx/ld
sin_a = ly/ld
self.force2par(p1,p2,cos_a,sin_a)
else:
#接触しない
p1.en[p2.n] = 0.0
p1.es[p2.n] = 0.0
#外力
for p in self.pe:
p.fy += -G * p.m #重力
def force2par(self,p1,p2,cos_a,sin_a):
kn = 10**7 #弾性係数（法線方向）
etan = 1000 #粘性係数（法線方向）
ks = 5000 #弾性係数（せん断方向）
etas = 1000 #粘性係数（せん断方向）
frc = 10 #摩擦係数
#相対的変位増分
un = +(p1.dx-p2.dx)*cos_a+(p1.dy-p2.dy)*sin_a
us = -(p1.dx-p2.dx)*sin_a+(p1.dy-p2.dy)*cos_a+(p1.r*p1.da+p2.r*p2.da)
#相対的速度増分
vn = +(p1.vx-p2.vx)*cos_a+(p1.vy-p2.vy)*sin_a
vs = -(p1.vx-p2.vx)*sin_a+(p1.vy-p2.vy)*cos_a+(p1.r*p1.va+p2.r*p2.va)
#合力
p1.en[p2.n] += kn*un
p1.es[p2.n] += ks*us
hn = p1.en[p2.n] + etan*vn
hs = p1.es[p2.n] + etas*vs
if hn <= 0.0:
#法線力がなければ、せん断力は０
hs = 0.0
elif abs(hs)-frc*hn >= 0.0:
#摩擦力以上のせん断力は働かない
hs = frc*abs(hn)*hs/abs(hs)
#合力
p1.fx += -hn*cos_a + hs*sin_a
p1.fy += -hn*sin_a - hs*cos_a
p1.fm -= p1.r*hs
def clacStep(self):
for p in self.pe:
p.fx = 0
p.fy = 0
p.fm = 0
self.calcForce()
for p in self.pe:
p.nextStep(self.dt)
self.step += 1
class MainWindow(QtGui.QMainWindow):
def __init__(self, *argv, **keywords ):
super(MainWindow,self).__init__(*argv,**keywords)
self.resize(300, 200)
self.view = QtGui.QGraphicsView(self)
self.scene = QtGui.QGraphicsScene(self.view)
self.view.setScene(self.scene)
self.setCentralWidget(self.view)
self.pen = QtGui.QPen()
self.pen.setColor(QtGui.QColor('#000000'))
self.dem = DEM()
self.redraw()
self.timer = QtCore.QTimer(self)
self.timer.setInterval(int(1000*self.dem.dt))
self.timer.timeout.connect(self.drawLoop)
self.timer.start()
def redraw(self):
for item in self.scene.items():
self.scene.removeItem(item)
pass
for p in self.dem.pe:
self.scene.addEllipse(p.x,-p.y,p.r*2,p.r*2,self.pen)
def drawLoop(self):
self.dem.clacStep()
self.redraw()
def main():
app = QtGui.QApplication(sys.argv)
w = MainWindow()
w.show()
sys.exit(app.exec_())
if __name__ == '__main__':
main()