法線方向の作用力

要素が接触している場合には、お互いに作用する力を計算します。 力の方向として、要素同士の法線方向とせん断方向（法線方向の垂直方向）の2方向を分けて考えます。
要素間の法線方向の力は線形バネによる反力（弾性力）とダンパーによる減衰を合わせたモデル化をします。

バネの弾性力は要素間の相対変位に弾性係数を乗じた値で、 ダンパーの粘性抵抗力は要素間の相対速度に粘性係数を乗じた値になります。
ただし、DEM（というより差分法）は、微小時間の間で生じた作用力の増分を 積算していく手法ですので、弾性力の増分を前のステップまでの弾性力に加算していきます。
粘性係数はその時の速度に応じて生じる力のため、積算しません。

バネの弾性力　　　　：⊿en = Kn×⊿un
ダンパーの粘性抵抗力：⊿dn = Dn×⊿un/⊿t
　　　　　　　　　　（⊿unは法線方向の相対変位量） 
法線方向の力　　　　：fn(t) = en(t-⊿t)+⊿en + ⊿dn

せん断方向の作用力

せん断方向のモデルは、線形バネとダンパーは法線方向と同じですが、 せん断方向の力は摩擦力以上は働かないので、上限値を制限するためのスライダーを追加で設けます。

バネの弾性力　　　　：⊿es = Ks×⊿us
ダンパーの粘性抵抗力：⊿ds = Ds×⊿us/⊿t
　　　　　　　　　　（⊿usはせん断方向の相対変位量） 
せん断方向の力　　　：fs(t) = es(t-⊿t)+⊿es + ⊿ds
　　　　　　　　　　　ただし、fs(t) ≦ u×fn(t)

ソースコード

 
class DEM:
 
    def _force_par2par(self,p1,p2,cos_a,sin_a):
 
        #相対的変位増分
        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)
 
        inf = self.interfaces.get(p1.etype,p2.etype)
 
        #合力
        p1.en[p2.n] += inf.kn*un
        p1.es[p2.n] += inf.ks*us
        hn = p1.en[p2.n] + inf.etan*un/self.dt
        hs = p1.es[p2.n] + inf.etas*us/self.dt
 
        if hn <= 0.0:
            #法線力がなければ、せん断力は０
            hs = 0.0
        elif abs(hs) >= inf.frc*hn:
            #摩擦力以上のせん断力は働かない
            hs = inf.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
        p2.fx +=  hn*cos_a - hs*sin_a
        p2.fy +=  hn*sin_a + hs*cos_a
        p2.fm -=  p2.r*hs