ObjectConnectorRollingDiscPenalty

A (flexible) connector representing a rolling rigid disc (marker 1) on a flat surface (marker 0, ground body, not moving) in global \(x\)-\(y\) plane. The connector is based on a penalty formulation and adds friction and slipping. The contraints works for discs as long as the disc axis and the plane normal vector are not parallel. Parameters may need to be adjusted for better convergence (e.g., dryFrictionProportionalZone). The formulation for the arbitrary disc axis is still under development and needs further testing. Note that the rolling body must have the reference point at the center of the disc.

Simpler: mbs.CreateRollingDiscPenalty adds this item, with what it needs, in one call.

Interface

Python names: ObjectConnectorRollingDiscPenalty or RollingDiscPenalty, and VRollingDiscPenalty for its visualization
Nodes it takes: NodeGenericData
Markers it acts on: those providing Position and Orientation: MarkerBodyRigid, MarkerNodeRigid, MarkerSuperElementRigid, MarkerKinematicTreeRigid

Parameters

The parameters of the item; in a dictionary, its type is ‘ConnectorRollingDiscPenalty’:

name	type	size	default	description
name	String		‘’	constraints’s unique name
markerNumbers	ArrayMarkerIndex	2	[ invalid (-1), invalid (-1) ]	(symbol: \([m0,m1]\tp\)) list of markers used in connector; \(m0\) represents a point at the plane surface (normal of surface plane defined by planeNormal); the ground can also be a moving rigid body; \(m1\) represents the rolling body, which has its reference point (=local position [0,0,0]) at the disc center point
nodeNumber	NodeIndex		invalid (-1)	(symbol: \(n_d\)) node number of a NodeGenericData (size=3) for 3 dataCoordinates, needed for discontinuous iteration (friction and contact)
discRadius	PReal		0.	defines the disc radius; must be given: the default is only a placeholder
discAxis	Vector3D	3	[1,0,0]	(symbol: \(\LU{m1}{\wv_{1}}, \;\; \vert \LU{m1}{\wv_{1}}\vert = 1\)) axis of disc defined in marker \(m1\) frame
planeNormal	Vector3D	3	[0,0,1]	(symbol: \(\LU{m0}{\vv_{PN}}, \;\; \vert \LU{m0}{\vv_{PN}}\vert = 1\)) normal to the contact / rolling plane (ground); note that the plane reference point can be arbitrarily chosen by the location of the marker \(m0\)
dryFrictionAngle	Real		0.	(symbol: \(\alpha_t\)) angle [SI:1 (rad)] which defines a rotation of the local tangential coordinates dry friction; this allows to model Mecanum wheels with specified roll angle
contactStiffness	UReal		0.	(symbol: \(k_c\)) normal contact stiffness [SI:N/m]
contactDamping	UReal		0.	(symbol: \(d_c\)) normal contact damping [SI:N/(m s)]
dryFriction	Vector2D	2	[0,0]	(symbol: \([\mu_x,\mu_y]\tp\)) dry friction coefficients [SI:1] in local marker 1 joint \(J1\) coordinates; if \(\alpha_t==0\), lateral direction \(l=x\) and forward direction \(f=y\); assuming a normal force \(f_n\), the local friction force can be computed as \(\LU{J1}{\vp{f_{t,x}}{f_{t,y}}} = \vp{\mu_x f_n}{\mu_y f_n}\)
dryFrictionProportionalZone	Real		0.	(symbol: \(v_\mu\)) limit velocity [m/s] up to which the friction is proportional to velocity (for regularization / avoid numerical oscillations)
viscousFriction	Vector2D	2	[0,0]	(symbol: \([d_x, d_y]\tp\)) viscous friction coefficients [SI:1/(m/s)] in local marker 1 joint \(J1\) coordinates; proportional to slipping velocity, leading to increasing slipping friction force for increasing slipping velocity
rollingViscousFriction	Real		0.	(symbol: \(\mu_r\)) viscous rolling friction [SI:s/m]: the force acts against the velocity of the trail on ground and is proportional to this velocity and to the contact normal force; currently, only implemented for disc axis parallel to ground!
useLinearProportionalZone	Bool		False	if True, a linear proportional zone is used; the linear zone performs better in implicit time integration as the Jacobian has a constant tangent in the sticking case
activeConnector	Bool		True	flag, which determines, if the connector is active; used to deactivate (temporarily) a connector or constraint
visualization	VObjectConnectorRollingDiscPenalty			parameters for visualization of item

Renamed parameters, still taken with a DeprecationWarning: rollingFrictionViscous (deprecated since 1.12.258, removed in 2031): use rollingViscousFriction.

Visualization parameters

The parameters of VObjectConnectorRollingDiscPenalty, given as visualization:

name	type	size	default	description
show	Bool		True	set true, if item is shown in visualization and false if it is not shown
discWidth	float		0.1	width of disc for drawing
color	Float4	4	[-1.,-1.,-1.,-1.]	RGBA connector color; if R==-1, use default color

Drawing

Drawn as all connectors; the disc of width discWidth, with 4 times general.cylinderTiling segments, and its axis.

Settings beyond those of all connectors: general.cylinderTiling.

Output variables

Available as OutputVariableType in sensors, Get...Output() and other functions:

output variable	symbol	description
Position	\(\LU{0}{\pv}_{G}\)	current global position of contact point between rolling disc and ground
Velocity	\(\LU{0}{\vv}_{trail}\)	current velocity of the trail (according to motion of the contact point along the trail!) in global coordinates; this is not the velocity of the contact point!
VelocityLocal	\(\LU{J1}{\vv}\)	relative slip velocity at contact point in special \(J1\) joint coordinates
ForceLocal	\(\LU{J1}{\fv} = \LU{0}{[f_{t,x},\, f_{t,y},\, f_{n}]\tp}\)	contact forces acting on disc, in special \(J1\) joint coordinates, see section Connector Forces, \(f_{t,x}\) being the lateral force (parallel to ground plane), \(f_{t,y}\) being the longitudinal force and \(f_{n}\) being the contact normal force
RotationMatrix	\(\LU{0,J1}{\Am} = [\LU{0}{\wv_{lat}},\, \LU{0}{\wv}_2,\, \LU{0}{\vv_{PN}}]\)	transformation matrix of special joint coordinates \(J1\) to global coordinates

Detailed description

Definition of quantities

intermediate variables	symbol	description
marker m0 position	\(\LU{0}{\pv}_{m0}\)	current global position which is provided by marker m0, any ground reference point; currently unused
marker m0 orientation	\(\LU{0,m0}{\Rot}\)	current rotation matrix provided by marker m0; currently unused
marker m1 position	\(\LU{0}{\pv}_{m1}\)	center of disc
marker m1 orientation	\(\LU{0,m1}{\Rot}\)	current rotation matrix provided by marker m1
data coordinates	\(\xv=[x_0,\,x_1,\,x_2]\tp\)	data coordinates for \([x_0,\,x_1]\): hold the sliding velocity in lateral and longitudinal direction of last discontinuous iteration; \(x_2\): represents gap of last discontinuous iteration (in contact normal direction)
marker m1 velocity	\(\LU{0}{\vv}_{m1}\)	accordingly
marker m1 angular velocity	\(\LU{0}{\tomega}_{m1}\)	current angular velocity vector provided by marker m1
ground normal vector	\(\LU{0}{\vv_{PN}} = \LU{0,m0}{\Am} \LU{m0}{\vv_{PN}}\)	normalized normal vector to the ground body (rotates with marker \(m0\) if not fixed to ground)
ground position B	\(\LU{0}{\pv}_{B}\)	disc center point projected on ground (normal projection)
ground position C	\(\LU{0}{\pv}_{C}\)	contact point of disc with ground
ground velocity C	\(\LU{0}{\vv}_{C}\)	velocity of disc at ground contact point (must be zero at end of iteration)
wheel axis vector	\(\LU{0}{\wv_1} =\LU{0,m1}{\Rot} \LU{m1}{\wv_{1}} \)	normalized disc axis vector in global coordinates
longitudinal vector	\(\LU{0}{\wv_2}\)	vector in longitudinal (motion) direction
contact point vector	\(\LU{0}{\wv_3}\)	normalized vector from disc center point in direction of contact point C
lateral vector	\(\LU{0}{\wv_{lat}} = \LU{0}{\vv_{PN}} \times \LU{0}{\wv}_2\)	vector in lateral direction, parallel to ground plane
\(D1\) transformation matrix	\(\LU{0,D1}{\Am} = [\LU{0}{\wv_1},\, \LU{0}{\wv_2},\, \LU{0}{\wv_3}]\)	transformation of special disc coordinates \(D1\) to global coordinates
connector forces	\(\LU{J1}{\fv}=[f_{t,x},\,f_{t,y},\,f_n]\tp\)	joint force vector at contact point in joint 1 coordinates: x=lateral direction, y=longitudinal direction, z=plane normal (contact normal)

Geometric relations

The main geometrical setup is shown in the following figure:

../../../_images/ObjectJointRollingDiscSketch.png

First, the contact point \(\LU{0}{\pv}_{C}\) must be computed. With the helper vector,

\[ \LU{0}{\xv} = \LU{0}{\wv}_1 \times \LU{0}{\vv_{PN}} \]

we create a disc coordinate system \(D1\) (\(\LU{0}{\wv}_1, \; \LU{0}{\wv}_2, \; \LU{0}{\wv}_3\)), with the longitudinal direction,

\[ \LU{0}{\wv}_2 = \frac{1}{|\LU{0}{\xv}|} \LU{0}{\xv} \]

and the vector to the contact point,

\[ \LU{0}{\wv}_3 = \LU{0}{\wv}_1 \times \LU{0}{\wv}_2 \]

The vector from marker \(m0\) position to the contact point can be computed from

\[ \LU{0}{\pv}_{C} = \LU{0}{\pv}_{m1} + r \cdot \LU{0}{\wv}_3 - \LU{0}{\pv}_{m0} \]

The velocity of the contact point at the disc is computed from,

\[ \LU{0}{\vv}_{C} = \LU{0}{\vv}_{m1} + \LU{0}{\tomega}_{m1} \times (r\cdot \LU{0}{\wv}_3) - \left(\LU{0}{\vv}_{m0} + \LU{0}{\tomega}_{m0} \times \LU{0}{\pv}_{C} \right) \]

A second coordinate system, denoted as \(J1\), is defined by vectors (\(\LU{0}{\wv}_{lat}, \; \LU{0}{\wv}_2, \; \LU{0}{\vv}_{PN}\)), using

\[ \LU{0}{\wv}_{lat} = \LU{0}{\vv_{PN}} \times \LU{0}{\wv}_2 \]

Note that in the case that the rolling axis \(\LU{0}{\wv}_1\) lies in the rolling plane, we obtain the special case \(\LU{0}{\wv}_{lat} = \LU{0}{\wv}_1\) and \(\LU{0}{\wv}_3 = -\LU{0}{\vv}_{PN}\).

Computation of normal and tangential forces

The connector forces at the contact point \(C\) are computed as follows. The normal contact force reads

\[ f_n = \left(k_c \cdot \LU{0}{\pv}_{C} + d_c \cdot \LU{0}{\vv}_{C} \right)\tp \LU{0}{\vv_{PN}} \, . \]

Note that due to the projection onto \(\LU{0}{\vv_{PN}}\), this equation also works for inclined planes and reference points, that are not at \([0,0,0]\tp\).

The inplane velocity in joint coordinates,

\[ \LU{J1}{\vv_t} = [\LU{0}{\vv}_{C}\tp \LU{0}{\wv}_{lat}, \; \LU{0}{\vv}_{C}\tp \LU{0}{\wv}_2 ]\tp \, , \]

is used for the computation of tangential forces,

\[ \LU{J1}{\fv_t} = [f_{t,x} ,\; f_{t,y}]\tp = \LU{J1}{\tmu} \cdot \left( \phi(|\vv_t|,v_\mu) \cdot f_n \cdot \LU{J1}{\ev_t} \right) \, , \]

with the regularization function, see Geradin and Cardona GeradinCardona2001 (Sec. 7.9.3), if useLinearProportionalZone=False,

\[\begin{split} \phi(v, v_\mu) = \left\{ \begin{array}{ccl} \displaystyle \left( 2-\frac{v}{v_\mu} \right)\frac{v}{v_\mu} & \mathrm{if} & v \le v_\mu \\ 1 & \mathrm{if} & v > v_\mu \\ \end{array} \right. \end{split}\]

and the linear regularization function, if useLinearProportionalZone=True,

\[\begin{split} \phi(v, v_\mu) = \left\{ \begin{array}{ccl} \displaystyle \frac{v}{v_\mu} & \mathrm{if} & v \le v_\mu \\ 1 & \mathrm{if} & v > v_\mu \\ \end{array} \right. \end{split}\]

The direction of tangential slip is given as

\[\begin{split} \LU{J1}{\ev_t} = \left\{ \begin{array}{ccl} \displaystyle \frac{\LU{J1}{\vv_t}}{|\vv_t|} &\mathrm{if}& |\vv_t|>0 \\ \vp{0}{0} &\mathrm{else}& \\ \end{array} \right. \end{split}\]

The friction coefficient matrix \(\LU{J1}{\tmu}\) is given in joint coordinates and computed from

\[ \LU{J1}{\tmu} = \mp{\mu_x + d_x \cdot |\vv_t|}{0}{0}{\mu_y + d_y \cdot |\vv_t|} \]

where for isotropic behaviour of surface and wheel, it will give a diagonal matrix with the friction coefficient in the diagonal. In case that the dry friction angle \(\alpha_t\) is not zero, the \(\tmu\) changes to

\[ \LU{J1}{\tmu} = \mp{\cos(\alpha_t)}{\sin(\alpha_t)}{-\sin(\alpha_t)}{\cos(\alpha_t)} \mp{\mu_x + d_x \cdot |\vv_t|}{0}{0}{\mu_y + d_y \cdot |\vv_t|} \mp{\cos(\alpha_t)}{-\sin(\alpha_t)}{\sin(\alpha_t)}{\cos(\alpha_t)} \]

Connector forces

Finally, the connector forces read in joint coordinates

(77)\[ \LU{J1}{\fv} = \vr{f_{t,x}}{f_{t,y}}{f_n} \]

and in global coordinates, they are computed from

\[ \LU{0}{\fv} = f_{t,x}\LU{0}{\wv}_{lat} + f_{t,y} \LU{0}{\wv}_2 + f_n \LU{0}{\vv}_{PN} \]

Due to the fact that the marker positions are not collocated with the contact point, there are additional torques that need to be considered in the action on the body. The torque onto the disc (marker \(m1\)) is computed as

\[ \LU{0}{\ttau_{m1}} = (r\cdot \LU{0}{\wv}_3) \times \LU{0}{\fv} \]

The torque onto the ground (marker \(m0\)) is computed as

\[ \LU{0}{\ttau_{m0}} = \LU{0}{\pv}_{C} \times \LU{0}{\fv} \]

Note that if activeConnector = False, we replace (77) with

\[ \LU{J1}{\fv} = \Null \]

Mini example

#a disc rolling on the ground plane with penalty contact and friction
r = 0.2
inertia = InertiaCylinder(density=1000, length=0.05, outerRadius=r, axis=0)
node = mbs.AddNode(NodeRigidBodyEP(referenceCoordinates=[0,0,r]+eulerParameters0,
                   initialVelocities=[0,-2,0]+list(AngularVelocity2EulerParameters_t([2/r,0,0], eulerParameters0))))
disc = mbs.AddObject(ObjectRigidBody(nodeNumber=node, mass=inertia.Mass(), inertia=inertia.GetInertia6D()))
mbs.AddLoad(LoadMassProportional(markerNumber=mbs.AddMarker(MarkerBodyMass(bodyNumber=disc)), loadVector=[0,0,-9.81]))
mGround = mbs.AddMarker(MarkerBodyRigid(bodyNumber=oGround, localPosition=[0,0,0]))
mDisc = mbs.AddMarker(MarkerBodyRigid(bodyNumber=disc, localPosition=[0,0,0]))
nData = mbs.AddNode(NodeGenericData(numberOfDataCoordinates=3, initialCoordinates=[0,0,0]))
mbs.AddObject(ObjectConnectorRollingDiscPenalty(markerNumbers=[mGround, mDisc], nodeNumber=nData, discRadius=r,
                                                discAxis=[1,0,0], planeNormal=[0,0,1], dryFriction=[0.5,0.5],
                                                dryFrictionProportionalZone=1e-2, contactStiffness=1e5, contactDamping=1e3))

mbs.Assemble()
mbs.SolveDynamic()

#it rolls on with the initial velocity: y = -2*t at t=1
exu.sys['testResult'] = mbs.GetNodeOutput(node, exu.OutputVariableType.Position)[1] #-2

Relevant Examples (Ex) and TestModels (TM) with weblink to github: bicycleIftommBenchmark.py (Ex), leggedRobot.py (Ex), mobileMecanumWheelRobotWithLidar.py (Ex), reinforcementLearningRobot.py (Ex), carRollingDiscTest.py (TM), laserScannerTest.py (TM), mecanumWheelRollingDiscTest.py (TM), rollingCoinPenaltyTest.py (TM), rollingDiscTangentialForces.py (TM), rotatingTableTest.py (TM), createFunctionsTest.py (TM), createRollingDiscPenaltyTest.py (TM)