MV-9003: LuGre Friction Tutorial

In this tutorial, you will learn how to implement and use a new MotionSolve object.

In this tutorial, you will learn how to:

Create a friction force on a translational joint using the LuGre friction model, and using Object Oriented Design with classes.

A class is a blueprint, or template, that allows you to create a higher level entity and treat it just like any other MotionSolve object.

In this example, you will combine a series of MotionSolve entities into a class called LuGre to define a new modeling element that exhibits new behavior. This process is called aggregation and MotionSolve provides a special composite class to represent such objects.

Composite is a special class that allows you to aggregate entities, such as parts, markers, forces, differential equations, and other MotionSolve primitives. A composite class object behaves just like any other solver object except it is a collection of entities conveniently packaged together. They behave as atomic entities and can be instantiated multiple times.

The advantage of using the composite class is that the LuGre object inherits behaviors from the base class. In particular, you will be able to modify the LuGre object at run-time, just like you would issue a modify command on any other MotionSolve primitive object, such as a part or a marker. The example below illustrates how you can create your own composite class. First, you'll examine the LuGre friction model.

LuGre Friction Model

The LuGre friction uses a bristle to model the friction force. Friction is modeled as the average deflection force of elastic springs between two mating surfaces. When a tangential motion is applied, the bristles deflect like springs. If the deflection is sufficiently large, the bristles start to slip. The slip velocity between the mating surfaces determines the average bristle deflection for a steady state motion. It is lower at low velocities, which implies that the steady state deflection decreases with increasing velocity.

The LuGre model is defined in terms of the following input parameters:

$V_{s}$ = Stiction to dynamic friction transition velocity
$μ_{s}$ = Coefficient of static friction
$μ_{d}$ = Coefficient of dynamic friction
$k_{0}$ = Bristle stiffness
$k_{1}$ = Bristle damping
$k_{2}$ = Viscous coefficient

These variables are used in the formulation:

Velocity in joint v = $V_{z} (I, J, J, J)$
Stribeck factor $P = - {(\frac{v}{V_{s}})}^{2}$
The normal force $N = \sqrt{F_{X} {(I, J, J)}^{2} + F_{Y} {(I, J, J)}^{2}}$
Coulomb friction $F c = μ d * N$
Static Friction $F s = μ d * N$
$G = \frac{F c + (F s - F c) * e^{P}}{k_{0}}$
ODE defining bristle deflection $\dot{z} = v - \frac{| v | * z}{G}$
The friction force $F = - (k_{0} * z + k_{1} * \dot{z} + k_{2} * v)$

The LuGre friction force consists of four atomic elements:

DIFF defining the bristle deflection
MARKER defining the point of action of the friction force
FORCE defining the friction element
REQUEST capturing the friction force on the block

Load the Msolve Module

Load the msolve module, which is done by issuing this command:

In [1]: from msolve import*

The above command requires the msolve module to be in your computer path. Assuming the above is successful, you have now imported the msolve names in the current namespace. This means that you have access to all the classes and functions defined in msolve and you can start creating the LuGre class and supporting model.

Below is the entire implementation of the class. Each section will be reviewed in the paragraphs following the code cell.

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)

Now that the class has been defined, you need to to fill the implementation details. The LuGre friction uses a single component force between two bodies, and a first-order differential equation that models the average deflection of the bristles as a function of slip velocity, coefficients of friction, and normal force. First you'll work on the Single Component Force. The goal of the Sforce function (normally called SfoSub) is to compute and return the scalar value of the friction force. The three component depicted here are computed separately in the force function.
```
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  
```
The composite class LuGre is inherited from the built-in composite class. You have to write two special methods, createChildren and setChildDataValue, to define the behavior of this class. You can also optionally write a validate method, validateSelf.
Each composite class must have two methods.
1. The first method is called createChildren. It is used to create all the MotionSolve entities that are needed by the composite class. When the composite class object is created, the createChildren method is invoked. This happens only once, at the creation of the composite object. In createChildren, you will instantiate various MotionSolve objects and define their immutable properties.
2. The updateChildren method is used to define all the mutable properties of the class and updates the class object with the new values. This facilitates modifying those values at runtime, when the solver is running and it can be done easily just like any other primitive MotionSolve objects. This is explained in more detail later in the tutorial.
An optional method called validateSelf can be defined. This method checks the input to ensure that the data provided to the LuGre class interface is physically meaningful. For example, by verifying that the static friction coefficient is larger than the dynamic friction coefficient.

After the methods described above, there are two additional functions. These are the Sforce and Diff user written subroutines. Note that these functions are created directly within the class definition. The advantage of doing this is to make the class properties available to the implementation of the subroutines.

This obviates the need to pass parameters to the user subroutines.

The Sforce subroutine computes the three components of the friction Sforce, as described in the image below.

Figure 2.

Note that the function's first argument is self because Sforce and Diff user written subroutines are Python callable class methods. Also, note that the signature of the Sforce function is in the familiar form (id, par, npar, and so on) but you are not actually unpacking the parameters, as there is no need to pass anything to the subroutine.

You are in fact accessing the properties of the LuGre instance by using the dot operator on the instance itself.

Similarly, you can define the differential equation. This also accesses the LuGre properties directly.

Create a Model

Now that you have completely defined the LuGre class and implemented the two methods (Sforce and Diff), you have to create a suitable model for testing.

You can use a simple friction element defined between two bodies on a sliding block. The sliding block is constrained by a translational joint and a forcing function acts upon it.

You can study the effect of changing the Static Friction Coefficient μ_s on the computed friction force by modifying it between simulations.

The model is created inside a function, called _slidingblock, in the code cell below. This is done purely for convenience.

Issue this command:

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

Note that the creation of the LuGre friction is done by calling one single line:

model.lugre = LuGre (joint=joint)

All you need to do is pass the Translational Joint, joint, to the LuGre class instantiation.

All of the properties that are used in the calculation of the Friction Force default to the values defined in the class interface. You can, of course, pass different values to the LuGre class but for now, you will use default values.

Test the model and plot some results.
1. Create a model by calling the following function:
```
model = sliding_block ()
```
2. Run a 4 second simulation and store the results of the run in an appropriate container:
```
run = model.simulate (type="DYNAMICS", returnResults=True, end=4, dtout=.01)
```
Now that you have a run container with a series of requests, you can selectively plot some of the relevant channels, such as the magnitude of the friction force, the block velocity,and its displacement as a function of time.

Import a plotting package:

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

Then, extract some of the data to plot. At this point, you can create a plot object and display three separate curves.

The friction force.
The velocity of the block.
The block displacement.

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)

To increase the friction static coefficient, set the LuGre $μ_{s t a t i c}$ to a different value and continue the simulation for an additional four seconds.
```
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)
```
Figure 4.

With the increased $μ_{s}$ , you can see that the forcing function cannot overcome the bristle resistance and the block is not able to move.