前回は下記の論文に参考に転がり抵抗を導入してみました。今回は、要素を土にするためにもう一つ必要な、粒子間の吸着力（土の世界でいう粘着力）を考慮したモールクーロンの破壊基準の適用にチャレンジします。
本来は、転がり抵抗モデルが意図した通りに組めているかの検証をしなくてはいけないところですが、やりたい・やれそうなことは、順序は関係なくやっていくスタンスです＾＾
地盤材料の破壊基準を表現するためのシンプルな個別要素モデル
何をするかと言うと、せん断に対する抵抗力を一定値ではなく、拘束力（直方向の力）に係数を掛けた値を上限とするモールクーロンの破壊基準を組み込みます。加えて、通常のDEMでは扱わない引張抵抗力（粘着力）を考慮するようなモデルにします。この引張力が転がり抵抗にも関与するので、それに応じて、転がり抵抗部分も変更します。
とにかくやってみよー。
ボンドモデル部分ソース（全ソースは最後に）
# 粒子間ボンドモデル
cdef double bs,tu=0
if i < peCount and j < peCount: #粒子同士
tu = pe[i].et[j]
if hn < -tu and tu > 0: #ボンド破壊条件
pe[i].et[j] = 0
tu = 0
print('bond break',i,j,hn,tu)
#モールクーロン破壊基準の適用
bs = (hn + tu) * tan(pe[i].ph)
if hn <= -tu:
hs = 0.0
elif fabs(hs) > fabs(bs):
hs = bs * hs/fabs(hs)
#hn = hn - tu
else: #粒子同士以外
if hn <= 0.0:
#法線力がなければ、せん断力は０
hs = 0.0
elif fabs(hs) >= inf.frc*hn:
#摩擦力以上のせん断力は働かない
hs = inf.frc*fabs(hn) *hs/fabs(hs)
上記論文ではDF条件と称して、接触状態にある粒子間では粘着力が発揮されていて、一定の引張力が加わったときに粘着がなくなり（ボンドが破壊）、その後の接触は通常の摩擦として扱うモデルになっています。
普通は離れている要素の粘着力は既にボンド破壊状態であり、粘着力は働かない（通常の摩擦であらわす）のですが、確かめるための便宜上、初期値では粘着力を有効にしておいて、破壊されるまでの間は、接触する要素間に粘着力が働くようにしておきます。
さて、どうなるか？
（左：粘着なし 右：粘着有り）
お、粘着力は働いてるっぽい。しかし、2つ問題が発生。動きが止まってきたあたりにバネのようにスライドする謎の動きと、後半、なぜか要素が離れていく・・・。
前者は、今回のボンドモデルに関係なく、線と円要素の接触処理によるところの部分の不具合ですかね。そういえば、前の転がり抵抗のときもそうでした。摩擦の上限を超えることなく、せん断バネでスライドしているせいかな？検証が必要ですなー。パラメータが非適正な値になっているだけだといいんだけど。
後者について、実はイメージで言うと、今回の動きは、僕の想像と違って、もっと粒子同士が重なって、重なった状態で安定するイメージでした。上のソースだと、粒子が反発して離れる過程で、引張力（粘着力）が働くけども、発揮される直方向の力自体の値には影響を与えません。収束状態は力の釣り合ったところなので、ボンドが破壊するまで反発が働き続けると、結果、離れてしまうのかなと。
ともすれば、最終的な力は、弾性力（反発力）に引張力を差し引いた値にするのでは？と思い、上のソースの最後に
hn = hn - tu
を加えてみました。すると、
イメージ通り、粒子同士は結着したものの、転がり抵抗力に余計な影響を与えるらしく、結着したまま回り続ける謎の動きにorz
転がり抵抗をOFFってみると、こんな感じ。
ありゃー、結着しながら転がっていくー。
転がり抵抗のモーメント算定式を修正するのと他に、粒子同士が結着した接触点が固定されるように回転を考慮しなくはいけないことが判明。
上記論文、モデルについて丁寧に書かれてくれている部類の論文だと思いますが、いざコード化すると、色々と不明な点が出てきます。さて、どうしようかな・・・。
# -*- coding: utf-8 -*-
from libc.math cimport sqrt,sin,cos,tan,atan2,fabs,ceil,floor
from cpython cimport bool
from cpython.mem cimport PyMem_Malloc, PyMem_Realloc, PyMem_Free
cdef double PI = 3.1415926535 #円周率
cdef double G = 9.80665 #重力加速度
class _Particle:
pass
cdef struct Particle:
#要素共通
int etype #要素タイプ
int n #要素No.
double r#半径(m)
double x#X座標(m)
double y#Y座標(m)
double a#角度(rad)
double dx#X方向増加量(m)
double dy#Y方向増加量(m)
double da#角度増加量(m)
double vx#X方向速度(m/s)
double vy#Y方向速度(m/s)
double va#角速度(rad/s)
double fx#X方向力(N)
double fy#Y方向力(N)
double fm#回転モーメント(Nm)
double *en #直方向弾性力(N)
double *es #せん断方向弾性力(N)
#粒子専用
double m #質量 (kg)
double Ir #慣性モーメント(Nm)
double ph #せん断抵抗角(rad)
double t #ボンド強度(N/m)
int *nearList #近傍リスト
int nearCount #近傍リスト数
int *co #ボンド判定
double *et #粒子間ボンド力(N)
class _Line:
pass
cdef struct Line:
#要素共通
int etype #要素タイプ
int n #要素No.
## double r#半径
## double x#X座標
## double y#Y座標
## double a#角度
## double dx#X方向増加量
## double dy#Y方向増加量
## double da#角度増加量
## double vx#X方向速度
## double vy#Y方向速度
## double va#角速度
## double fy
## double fx
## double fm
## #double *en #弾性力（直方向）
## #double *es #弾性力（せん断方向）
#線要素専用
double x1
double y1
double x2
double y2
class _Interface:
pass
cdef struct Interface:
double kn #弾性係数（法線方向）
double etan#粘性係数（法線方向）
double ks#弾性係数（せん断方向）
double etas#弾性係数（せん断方向）
double frc#摩擦係数
#グローバル変数
cdef Particle *pe
cdef Line *le
cdef Interface infs[9]
cdef int infNo[9][9]
cdef int area[4] # 解析範囲
cdef int elCount #全要素数
cdef int peCount #粒子の数
cdef int leCount #線要素の数
cdef double dt = 0.001 #計算間隔
cdef int st = 0 #ステップ
cdef double rho = 10 #粒子間密度
def _nextStep():
pass
cdef void nextStep():
global st,nearUpdate
if nearUpdate == True:
cellRegist()
resetForce()
calcForce()
if nearUpdate == True:
nearUpdate = False
updateCoord()
st += 1
def _resetForce():
pass
cdef void resetForce():
cdef int i
for i in range(peCount):
pe[i].fx = 0
pe[i].fy = 0
pe[i].fm = 0
def _interface():
pass
cdef Interface interface(int etype1,int etype2):
cdef Interface inf
cdef int n
n = infNo[etype1][etype2]
inf.kn = infs[n].kn
inf.etan = infs[n].etan
inf.ks = infs[n].ks
inf.etas = infs[n].etas
inf.frc = infs[n].frc
return inf
def _calcForce():
pass
cdef void calcForce():
#2粒子間の接触判定
cdef double lx,ly,ld,cos_a,sin_a
cdef int i,j,n,m
for i in range(peCount):
#粒子登録法による近傍リストの絞込み
if nearUpdate == True:
m = 0
for n in range(pe[i].nearCount):
j = pe[i].nearList[n]
rc = pe[i].r + pe[j].r + nearAlpha * r_min*2
lx = pe[j].x - pe[i].x
ly = pe[j].y - pe[i].y
ld = (lx**2+ly**2)**0.5
if ld < rc:
nearList[m] = j
m += 1
pe[i].nearCount = m
for n in range(pe[i].nearCount):
pe[i].nearList[n] = nearList[n]
#for j in range(peCount):
for n in range(pe[i].nearCount):
j = pe[i].nearList[n]
if i == j:
continue
lx = pe[j].x - pe[i].x
ly = pe[j].y - pe[i].y
ld = (lx**2+ly**2)**0.5
if ld<(pe[i].r+pe[j].r):
cos_a = lx/ld
sin_a = ly/ld
force2par(i,j,cos_a,sin_a,ld*pe[i].r/(pe[i].r+pe[j].r))
else:
pe[i].en[j] = 0.0
pe[i].es[j] = 0.0
#粒子と壁の接触判定
cdef double cond[4]
cdef double xs[4]
cdef double ys[4]
for i in range(peCount):
cond = [pe[i].x-pe[i].r-area[0],-(pe[i].x+pe[i].r-area[1]),
pe[i].y-pe[i].r-area[2],-(pe[i].y+pe[i].r-area[3])]
xs = [area[0],area[1],pe[i].x,pe[i].x]
ys = [pe[i].y,pe[i].y,area[2],area[3]]
for j in range(4):
n = peCount + leCount + j
if cond[j] < 0:
lx = xs[j] - pe[i].x
ly = ys[j] - pe[i].y
ld = sqrt(lx**2+ly**2)
cos_a = lx/ld
sin_a = ly/ld
force2par(i,n,cos_a,sin_a,ld)
else:
pe[i].en[n] = 0.0
pe[i].es[n] = 0.0
#粒子と線の接触判定
cdef bool hit
cdef double x,y,a,d,b,s
for i in range(peCount):
for j in range(leCount):
hit = False
th0 = atan2(le[j].y2-le[j].y1, le[j].x2-le[j].x1)
th1 = atan2(pe[i].y -le[j].y1, pe[i].x -le[j].x1)
a = sqrt((pe[i].x-le[j].x1)**2+(pe[i].y-le[j].y1)**2)
d = fabs(a*sin(th1-th0))
if d < pe[i].r:
b = sqrt((pe[i].x -le[j].x2)**2+(pe[i].y -le[j].y2)**2)
s = sqrt((le[j].x2-le[j].x1)**2+(le[j].y2-le[j].y1)**2)
if a < s and b < s:
s = sqrt(a**2-d**2)
x = le[j].x1 + s*cos(th0)
y = le[j].y1 + s*sin(th0)
hit = True
elif a < b and a < pe[i].r:
x = le[j].x1
y = le[j].y1
hit = True
elif b < pe[i].r:
x = le[j].x2
y = le[j].y2
hit = True
if hit:
lx = x - pe[i].x
ly = y - pe[i].y
ld = sqrt(lx**2+ly**2)
cos_a = lx/ld
sin_a = ly/ld
force2par(i,le[j].n,cos_a,sin_a,ld)
else:
pe[i].en[le[j].n] = 0.0
pe[i].es[le[j].n] = 0.0
#外力
for i in range(peCount):
pe[i].fy += -G*pe[i].m #重力
def _force2par():
pass
cdef void force2par(int i,int j,double cos_a,double sin_a,double rd):
cdef int m,n
cdef int mat[2]
cdef double dx[2]
cdef double dy[2]
cdef double da[2]
cdef double r[2]
cdef double un,us,vn,vs,hn,hs
cdef Interface inf
for m in range(2):
if m == 0 : n = i
else: n = j
if n < peCount:
dx[m] = pe[n].dx
dy[m] = pe[n].dy
da[m] = pe[n].da
r[m] = pe[n].r
mat[m] = pe[n].etype
else:
dx[m] = 0
dy[m] = 0
da[m] = 0
r[m] = 0
mat[m] = 2
#相対的変位増分
un = +(dx[0]-dx[1])*cos_a+(dy[0]-dy[1])*sin_a
us = -(dx[0]-dx[1])*sin_a+(dy[0]-dy[1])*cos_a+(r[0]*da[0]+r[1]*da[1])
#相対的速度増分
#vn = +(pe[0].vx-pe[1].vx)*cos_a+(pe[0].vy-pe[1].vy)*sin_a
#vs = -(pe[0].vx-pe[1].vx)*sin_a+(pe[0].vy-pe[1].vy)*cos_a+(pe[0].r*pe[0].va+pe[1].r*pe[1].va)
inf = interface(mat[0],mat[1])
#合力（局所座標系）
pe[i].en[j] += inf.kn*un
pe[i].es[j] += inf.ks*us
hn = pe[i].en[j] + inf.etan*un/dt #圧縮が正
hs = pe[i].es[j] + inf.etas*us/dt #時計回りが正
# 粒子間ボンドモデル
cdef double bs,tu=0
if i < peCount and j < peCount: #粒子同士
tu = pe[i].et[j]
if hn < -tu and tu > 0: #ボンド破壊条件
pe[i].et[j] = 0
tu = 0
print('bond break',i,j,hn,tu)
#モールクーロン破壊基準の適用
bs = (hn + tu) * tan(pe[i].ph)
if hn <= -tu:
hs = 0.0
elif fabs(hs) > fabs(bs):
hs = bs * hs/fabs(hs)
#hn = hn - tu
else: #粒子同士以外
if hn <= 0.0:
#法線力がなければ、せん断力は０
hs = 0.0
elif fabs(hs) >= inf.frc*hn:
#摩擦力以上のせん断力は働かない
hs = inf.frc*fabs(hn) *hs/fabs(hs)
#転がり摩擦抵抗
cdef double mr = 0
cdef double b = 0.2
#if hn > -tu and hs != 0:
if hn > 0 and hs != 0:
mr = (hn + tu) * b * rd
mr = fabs(mr) * hs/fabs(hs)#回転を打ち消す方向に作用
if fabs(hs*pe[i].r) < fabs(mr):
#上限設定（逆回転するようなMrは発生しない）
mr = hs*pe[i].r
#print('fm',i,j,hs*pe[i].r-mr)
#test
#hn = hn - tu
#全体合力（全体座標系）
pe[i].fx += -hn*cos_a + hs*sin_a
pe[i].fy += -hn*sin_a - hs*cos_a
pe[i].fm -= hs*pe[i].r - mr
def _updateCoord():
pass
cdef void updateCoord():
global xbook
cdef double ax,ay,aa,v,vmax
vmax = 0
cdef int i
for i in range(peCount):
#位置更新（オイラー差分）
ax = pe[i].fx/pe[i].m
ay = pe[i].fy/pe[i].m
aa = pe[i].fm/pe[i].Ir
pe[i].vx += ax*dt
pe[i].vy += ay*dt
pe[i].va += aa*dt
pe[i].dx = pe[i].vx*dt
pe[i].dy = pe[i].vy*dt
pe[i].da = pe[i].va*dt
pe[i].x += pe[i].dx
pe[i].y += pe[i].dy
pe[i].a += pe[i].da
#以下近傍リスト用処理
v = sqrt(pe[i].vx**2+pe[i].vy**2)
if vmax < v:
vmax = v
xbook += vmax * dt
if xbook > nearAlpha*r_min:
global nearUpdate
nearUpdate = True
xbook = 0
#セル分割法と粒子登録法
##cell = []
##cellCount = []
cdef int ***cell
cdef int **cellCount
cdef int cellInfo[3]
cdef bool nearUpdate = False
cdef double nearAlpha = 0.1
cdef double xbook
cdef int *nearList
def _cellInit():
pass
cdef void cellInit():
global cell,cellCount,nearList
cell_width = r_max*2 + nearAlpha * r_min*2
xn = int(ceil((area[1]-area[0])/cell_width))
yn = int(ceil((area[3]-area[2])/cell_width))
cn = int(ceil(cell_width/r_min)+1)**2 #+1予備（重なりを考慮）
## cell = [[[-1 for i in range(cn)] for j in range(yn)] for k in range(xn)]
## cellCount = [[0 for i in range(yn)] for j in range(xn)]
cell = PyMem_Malloc(xn * sizeof(int*))
cellCount = PyMem_Malloc(xn * sizeof(int*))
for i in range(xn):
cell[i] = PyMem_Malloc(yn * sizeof(int*))
cellCount[i] = PyMem_Malloc(yn * sizeof(int))
for j in range(yn):
cell[i][j] = PyMem_Malloc(cn * sizeof(int))
cellInfo[0] = xn #列数
cellInfo[1] = yn #行数
cellInfo[2] = cn #セル格納最大数
cellReset()
nearList = PyMem_Malloc(cn * sizeof(int))
def _cellRegist():
pass
cdef void cellRegist():
cdef double cell_width
cdef int xn,yn,cn,n,i,k,j,l
cellReset()
cell_width = r_max*2 + nearAlpha * r_min*2
#セルに要素No格納
for i in range(peCount):
xn = int(floor(pe[i].x/cell_width))
yn = int(floor(pe[i].y/cell_width))
if xn<0 or cellInfo[0]<=xn or yn<0 or cellInfo[1]<=yn:
#print('Error! particle coord is out of range.')
#print('Particle n=%d x=%0.3f y=%0.3f' % (i,pe[i].x,pe[i].y))
continue
cn = cellCount[xn][yn]
cell[xn][yn][cn] = i
cellCount[xn][yn] += 1
#近傍リスト作成
for i in range(peCount):
xn = int(floor(pe[i].x/cell_width))
yn = int(floor(pe[i].y/cell_width))
n = 0
for j in range(xn-1,xn+2):
for k in range(yn-1,yn+2):
if j < 0 or cellInfo[0] <= j:
continue
if k < 0 or cellInfo[1] <= k:
continue
for l in range(cellCount[j][k]):
pe[i].nearList[n] = cell[j][k][l]
n += 1
pe[i].nearCount = n
cdef void cellReset():
cdef int i,j,k
for i in range(cellInfo[0]):
for j in range(cellInfo[1]):
cellCount[i][j] = 0
for k in range(cellInfo[2]):
cell[i][j][k] = -1
# -------------------------
# Pythonからの設定用
# -------------------------
def setNumberOfParticle(n):
global peCount,pe
peCount = n
pe = PyMem_Malloc(n * sizeof(Particle))
if not pe:
raise MemoryError()
def setNumberOfLine(n):
global leCount,le
leCount = n
le = PyMem_Malloc(n * sizeof(Line))
if not le:
raise MemoryError()
def initialize():
global elCount
elCount = peCount + leCount + 4
cdef int i,j
for i in range(peCount):
pe[i].etype = 1
pe[i].n = i
pe[i].r = 0
pe[i].x = 0; pe[i].y = 0; pe[i].a = 0
pe[i].dx = 0
pe[i].dy = 0
pe[i].da = 0
pe[i].vx = 0
pe[i].vy = 0
pe[i].va = 0
pe[i].fx = 0
pe[i].fy = 0
pe[i].fm = 0
pe[i].m = 0
pe[i].Ir = 0
pe[i].en = PyMem_Malloc(elCount * sizeof(double))
pe[i].es = PyMem_Malloc(elCount * sizeof(double))
pe[i].co = PyMem_Malloc(elCount * sizeof(int))
pe[i].et = PyMem_Malloc(elCount * sizeof(double))
for j in range(elCount):
pe[i].en[j] = 0
pe[i].es[j] = 0
pe[i].et[j] = 0
if j < peCount: pe[i].co[j] = True
else: pe[i].co[j] = False
for i in range(leCount):
le[i].etype = 2
le[i].n = peCount+i
def setDeltaTime(sec):
global dt
dt = sec
def setParticle(pe_no,pe_obj):
pe[pe_no].r = pe_obj.r
pe[pe_no].x = pe_obj.x
pe[pe_no].y = pe_obj.y
pe[pe_no].a = pe_obj.a
pe[pe_no].dx = pe_obj.dx
pe[pe_no].dy = pe_obj.dy
pe[pe_no].da = pe_obj.da
pe[pe_no].vx = pe_obj.vx
pe[pe_no].vy = pe_obj.vy
pe[pe_no].va = pe_obj.va
pe[pe_no].fx = pe_obj.fx
pe[pe_no].fy = pe_obj.fy
pe[pe_no].fm = pe_obj.fm
pe[pe_no].m = pe_obj.m
pe[pe_no].Ir = pe_obj.Ir
pe[pe_no].ph = pe_obj.ph
pe[pe_no].t = pe_obj.t
cdef int i
for i in range(elCount):
if i < len(pe_obj.en):
pe[pe_no].en[i] = pe_obj.en[i]
pe[pe_no].es[i] = pe_obj.es[i]
def particle(pe_no,pe_obj):
pe_obj.n = pe[pe_no].n
pe_obj.r = pe[pe_no].r
pe_obj.x = pe[pe_no].x
pe_obj.y = pe[pe_no].y
pe_obj.a = pe[pe_no].a
pe_obj.dx = pe[pe_no].dx
pe_obj.dy = pe[pe_no].dy
pe_obj.da = pe[pe_no].da
pe_obj.vx = pe[pe_no].vx
pe_obj.vy = pe[pe_no].vy
pe_obj.va = pe[pe_no].va
pe_obj.fx = pe[pe_no].fx
pe_obj.fy = pe[pe_no].fy
pe_obj.fm = pe[pe_no].fm
pe_obj.m = pe[pe_no].m
pe_obj.Ir = pe[pe_no].Ir
pe_obj.ph = pe[pe_no].ph
pe_obj.t = pe[pe_no].t
pe_obj.en = [0.0 for i in range(elCount)]
pe_obj.es = [0.0 for i in range(elCount)]
for i in range(elCount):
pe_obj.en[i] = pe[pe_no].en[i]
pe_obj.es[i] = pe[pe_no].es[i]
return pe_obj
def setLine(l_no,l_obj):
le[l_no].x1 = l_obj.x1
le[l_no].y1 = l_obj.y1
le[l_no].x2 = l_obj.x2
le[l_no].y2 = l_obj.y2
return l_obj
def line(l_no,l_obj):
l_obj.x1 = le[l_no].x1
l_obj.y1 = le[l_no].y1
l_obj.x2 = le[l_no].x2
l_obj.y2 = le[l_no].y2
return l_no
def setArea(x_min,x_max,y_min,y_max):
global area
area[0] = x_min
area[1] = x_max
area[2] = y_min
area[3] = y_max
def setInterface(mat1,mat2,inf_no,inf_obj):
infNo[mat1][mat2] = inf_no
infs[inf_no].kn = inf_obj.kn
infs[inf_no].etan = inf_obj.etan
infs[inf_no].ks = inf_obj.ks
infs[inf_no].etas = inf_obj.etas
infs[inf_no].frc = inf_obj.frc
cdef double r_max #粒子の最大半径
cdef double r_min #粒子の最小半径
def setup():
global r_max,r_min
cdef int i
r_max = 0
r_min = float('inf')
for i in range(peCount):
if pe[i].r > r_max:
r_max = pe[i].r
if pe[i].r < r_min:
r_min = pe[i].r
#セル分割法設定
cellInit()
pn = cellInfo[2] * 9
for i in range(peCount):
pe[i].nearList = PyMem_Malloc(pn * sizeof(int))
cellRegist()
#粒子間ボンド力設定
for i in range(peCount):
for j in range(peCount):
pe[i].et[j] = pe[i].t * (pe[i].r + pe[j].r)
def step():
return st
def calcStep(int n=1):
cdef int i
for i in range(n):
nextStep()