MV-9003：LuGre摩擦のチュートリアル

In [1]: from msolve import*

上記のコマンドを実行するには、msolveモジュールがコンピュータのパスになければなりません。上記が成功したと仮定すると、msolveの名称は現在のネームスペースにインポートされています。これにより、msolveで定義されているすべてのクラスと関数にアクセスできるようになり、LuGreクラスの作成とモデルのサポートを開始することができます。

In [2]: # Create a LuGre Composite Class
class LuGre(Composite):
  joint  = Reference (Joint)
  vs     = Double (1.E-3)
  mus    = Double (0.3)
  mud    = Double (0.2)
  k0     = Double (1e5)
  k1     = Double (math.sqrt(1e5))
  k2     = Double (0.4)
  def createChildren (self):
    """This is called when the object is created so the children objects
    """
    # The DIFF defining bristle deflection
    self.diff = Diff(routine=self.lugre_diff, ic=[0,0])

    # The MARKER on which the friction force acts
    self.im = Marker()

    # The FORCE defining the friction force
    self.friction = Sforce  (type="TRANSLATION",
    actiononly = True,
    routine = self.lugre_force,
    )

    # The REQUEST capturing the friction force
    self.request = Request (type="FORCE", comment="Friction force")
def updateChildren(self):

    # Set topology
    self.friction.i =self.joint.iself.friction.j =self.joint.jself.im.setValues (body=self.joint.i.body, qp=self.joint.i.qp, zp=self.joint.i.zp)
    self.request.setValues (i=self.im, j=self.joint.j, rm=self.joint.j)
def validateSelf(self, validator):
    validator.checkGe0 ("vs")
    validator.checkGe0 ("mus")
    validator.checkGe0 ("mud")
    validator.checkGe0 ("k0")
    validator.checkGe0 ("k1")
    validator.checkGe0 ("k2")
    if self.mud > self.mus:
      msg = tr("Mu dynamic({0}) needs to be smaller than Mu static ({1})",
                self.mud,self.mus)
      validator.sendError(msg)

def lugre_diff(self, id, time, par, npar, dflag, iflag):
"Diff user function"
    i = self.joint.i
    j   = self.joint.j
    vs  = self.vs
    mus = self.mus
    mud = self.mud
    k0  = self.k0

    z   = DIF(self)
    v   = VZ(i,j,j,j)

    N   = math.sqrt (FX(i,j,j)**2 + FY(i,j,j)**2)
    fs  = mus*N
    fc  = mud*N

    p   = -(v/vs)**2
    g   = (fc + (fs - fc) * math.exp(p))/k0
    def smooth_fabs(x):
       returnmath.fabs(x) # temp solution
    if iflag or math.fabs(g) <1e-8:
       return v
    else:
       return v - smooth_fabs(v) * z / g
def lugre_force (self,id, time, par, npar, dflag, iflag):
    "Friction Force user function"
    i    = self.joint.i
    j    = self.joint.j
    diff = self.diff
    k0   = self.k0
    k1   = self.k1
    k2   = self.k2

    v    = VZ(i,j,j,j)
    z    = DIF(diff)
    zdot = DIF1(diff)

    F    = k0*z + k1*zdot + k2*v
#print "Time=", time, "  V=", v, "  Z=", z, "  ZDOT=", zdot, "  F=", F
#print "k0*z=", k0*z, "  k1*zdot=", k1*zdot, "  k2*v=", k2*v
return -F  

複合クラスLuGreは、組み込みの複合クラスから継承されます。2つの特別なメソッドcreateChildrenおよびsetChildDataValueを記述して、このクラスの挙動を定義する必要があります。オプションとして検証メソッドvalidateSelfを記述することもできます。

validateSelfというオプションのメソッドを定義することができます。このメソッドは、入力をチェックして、LuGreクラスに渡されたデータが物理的に意味があることを確認します。例えば、静摩擦係数が動摩擦係数よりも大きいことを検証することで、それを確認することができます。

上記の方法の後には、2つの機能が追加されます。それらは、SforceおよびDiffユーザーサブルーチンです。これらの関数は、クラス定義内で直接作成されることにご留意ください。これを行うメリットは、クラスのプロパティをサブルーチンの実装で利用できるようにすることです。

これにより、ユーザのサブルーチンにパラメータを渡す必要がなくなります。

Sforceサブルーチンは、下の画像で示すように、摩擦Sforceの3つの成分を計算します。

SforceおよびDiffユーザーサブルーチンは、Python の呼び出し可能なクラスメソッドなので、関数の第1引数はselfです。また、Sforce関数の署名はおなじみの形式（id、par、nparなど）になっていますが、サブルーチンに何も渡す必要がないので、実際にはパラメータを解凍していないことにも注意が必要です。

実際には、インスタンス自体にドット演算子を使用してLuGreインスタンスのプロパティにアクセスしています。

同様に、微分方程式を定義することもできます。これはまた、LuGreプロパティに直接アクセスします。

In [2]: 
################################################################################ 
# Model definition                                                             #
################################################################################
def sliding_block ():
   """Test case for the LuGre friction model
        Slide a block across a surface.
        The block is attached to ground with a translational joint that has
        LuGre friction
     """
  model = Model ()
  
  Units      (mass="KILOGRAM", length="METER", time="SECOND", force="NEWTON")
  Accgrav    (jgrav=-9.800)
  Integrator (error=1e-5)
  Output     (reqsave=True)

# Location for Part cm and markers used for the joint markers
  qp = Point (10,0,0)
  xp = Point (10,1,0)
  zp = Point (11,0,0)

  ground    = Part   (ground=True)
  ground.jm = Marker (part=ground, qp=qp, zp=zp, xp=xp, label="Joint Marker")

  block    = Part   (mass=1, ip=[1e3,1e3,1e3], label="Block")
  block.cm = Marker (part=block, qp=qp, zp=zp, xp=xp, label="Block CM")

  # Attach the block to ground
  joint = Joint (  type="TRANSLATIONAL",
    i    = block.cm,
    j    = ground.jm,
  )

   # Apply the LuGre friction to the joint
  model.lugre = LuGre (joint=joint)

# Push the block
  Sforce (
    type="TRANSLATION",
    actiononly =True,
    i          = block.cm,
    j          = ground.jm,
    function   ="3*sin(2*pi*time)",
    label      ="Actuation Force",
  )

# Request the DISPLACEMENT, VELOCITY and FORCE of the Joint
  model.r1 = Request (
  type ="DISPLACEMENT",
    i       = block.cm,
    j       = ground.jm,
    rm      = ground.jm,
    comment ="Joint displacement",
  )
  model.r2 = Request (
   type="VELOCITY",
    i       = block.cm,
    j       = ground.jm,
    rm      = ground.jm,
    comment ="Joint velocity",
  )
  model.r3 = Request (
    type="FORCE",
    i       = block.cm,
    j       = ground.jm,
    rm      = ground.jm,
    comment ="Joint force",
  )
  return model

model.lugre = LuGre (joint=joint)

必要なのは、並進ジョイント jointをLuGreクラスのインスタンスに渡すことだけです。

摩擦力の計算に使用されるすべてのプロパティは、クラスインターフェースで定義された値がデフォルトです。もちろん、LuGreクラスに異なる値を渡すこともできますが、今のところはデフォルト値を使用します。

# this line causes the plots to be inlined on this web page
%matplotlib inline
import matplotlib.pyplot as plt

def plot(model):# get the channels from the model...
  disp = run.getObject (model.r1)
  velo = run.getObject (model.r3)
  force = run.getObject (model.r2)
  
plt.subplot(3, 1, 1)
  plt.plot(force.times, force.getComponent (3))
  plt.title('Friction force')
  plt.ylabel('force [N]')
  
plt.subplot(3, 1, 2)
  plt.plot(velo.times, velo.getComponent (3),'r')
  plt.title('Block Velocity')
  plt.ylabel('velocity [m/s]')
  plt.subplot(3, 1, 3)
  plt.plot(disp.times, disp.getComponent (3),'g')
  plt.title('Block Displacement')
  plt.xlabel('time (s)')
  plt.ylabel('travel [m]')

  plt.grid(True)
  plt.tight_layout()

# Show the plot
  plt.show()

In [7]: plot(model)

In [8]:  model.lugre.mus= 0.4
run = model.simulate (type="DYNAMICS", returnResults=True, duration=4, dtout=.01)

# replot the quantities and observe the changes to the Friction Force
plot(model)

${\mu }_s$ を大きくしたことで、強制関数がブリッスルの抵抗に打ち勝つことができず、ブロックを動かすことができないことがわかります。