LIBRARY

`PacejkaWheelJoint`

Joint for a wheel with Pacejka magic-formula friction rolling on the flat plane y = 0.

The magic-formula counterpart of SlipWheelJoint. The contact patch can slip and the resulting tangential force is computed from independent longitudinal and lateral Pacejka magic-formula curves evaluated by MultibodyComponents.pacejka_magic on the corresponding slip-velocity components. Each direction has its own stiffness/shape/peak/curvature factors (Bx,Cx,Dx,Ex longitudinal and By,Cy,Dy,Ey lateral). The friction curve is smooth in all derivatives, so both explicit RK and BDF-type integrators work.

Normal force

The wheel cannot leave the ground. Make sure that the normal force f_n never becomes negative.

This component extends from Renderable This component extends from CutJoint

Usage

MultibodyComponents.PacejkaWheelJoint(render=true, color=[1, 0, 0, 1], specular_coefficient=1.5, radius=0.3, width=0.035, Bx=10, Cx=1.9, Dx=1.0, Ex=0.97, By=10, Cy=1.3, Dy=0.9, Ey=0.97, v_small=1e-5)

Parameters:

Name	Description	Units	Default value
`iscut`		–	false
`residual`		–	zeros(3)
`surface`	Use an external road surface supplied via surface_frame (else flat ground y=0)	–	false
`angular_state`	Include the wheel's angular state (Cardan angles). Set false when the wheel is mounted on an already-rooted axis: the orientation is then taken from frame_a and no orientation loop/state is introduced (avoids redundant constraint forces in free-floating multi-wheel assemblies).	–	true
`sequence`	Rotation axis sequence for the Cardan angles that orient `frame_a` (default y, z, x)	–	[2, 3, 1]
`friction_ellipse`	Couple longitudinal and lateral forces through the combined friction ellipse (semi-axes `Dx`·f_n, `Dy`·f_n). When true (default), grip is shared between longitudinal and lateral demand, so the resultant force cannot exceed the friction limit under simultaneous slip. When false, the two magic-formula curves act independently (the resultant can modestly exceed the per-axis peak).	–	true
`render`		–	true
`color`		–	[1, 0, 0, 1]
`specular_coefficient`		–	1.5
`radius`	Radius of the wheel	m	0.3
`width`	Width of the wheel in animations	–	0.035
`Bx`	Longitudinal magic-formula stiffness factor (initial slope)	–	10
`Cx`	Longitudinal magic-formula shape factor	–	1.9
`Dx`	Longitudinal magic-formula peak friction coefficient	–	1.0
`Ex`	Longitudinal magic-formula curvature factor	–	0.97
`By`	Lateral magic-formula stiffness factor (initial slope)	–	10
`Cy`	Lateral magic-formula shape factor	–	1.3
`Dy`	Lateral magic-formula peak friction coefficient	–	0.9
`Ey`	Lateral magic-formula curvature factor	–	0.97
`v_small`	Small value added to `v_slip` to avoid division by zero in the slip model	–	1e-5

Connectors

frame_a - Frame3D is the fundamental 3D connector used for 6DOF motion. Most components have one or several Frame

connectors that can be connected together (Frame3D)

surface_frame - Frame3D is the fundamental 3D connector used for 6DOF motion. Most components have one or several Frame

connectors that can be connected together (Frame3D)

Variables

Name	Description	Units
`Rleaf`		–
`Rroot`		–
`x`	x-position of the wheel axis	m
`y`	y-position of the wheel axis	m
`z`	z-position of the wheel axis	m
`angles`	Angles that rotate the world frame into `frame_a` around the y-, z-, x-axis	rad
`der_angles`	Time derivatives of `angles`	rad/s
`phi_roll`	Wheel rolling angle (`angles[2]`)	rad
`w_roll`	Wheel rolling angular velocity	rad/s
`r_road_0`	Position vector from world frame to contact point on the road, resolved in world frame	m
`f_wheel_0`	Force vector on the wheel at the contact point, resolved in world frame	–
`f_n`	Contact force acting on the wheel in the normal direction	N
`f_lat`	Contact force acting on the wheel in the lateral direction	N
`f_long`	Contact force acting on the wheel in the longitudinal direction	N
`e_axis_0`	Unit vector along the wheel axis, resolved in world frame	–
`delta_0`	Distance vector from wheel center to contact point, resolved in world frame	–
`e_n_0`	Unit vector normal to the road at the contact point, resolved in world frame	–
`e_lat_0`	Unit vector in the lateral direction of the road at the contact point, resolved in world frame	–
`e_long_0`	Unit vector in the longitudinal direction of the road at the contact point, resolved in world frame	–
`e_s_0`	Road heading at (s, w), resolved in world frame (unit vector)	–
`s`	Road surface parameter 1	–
`w`	Road surface parameter 2	–
`v_0`	Velocity of wheel center, resolved in world frame	m/s
`w_0`	Angular velocity of the wheel, resolved in world frame	rad/s
`vContact_0`	Velocity of the contact point, resolved in world frame	m/s
`aux`	Auxiliary vector for building the contact frame	–
`v_slip_long`	Slip velocity in the longitudinal direction	m/s
`v_slip_lat`	Slip velocity in the lateral direction	m/s
`slip_ratio`	Slip ratio = v_slip_long / (v_0 · e_long_0)	–
`mu_x`	Longitudinal pure-slip friction coefficient from the magic formula	–
`mu_y`	Lateral pure-slip friction coefficient from the magic formula	–
`ellipse_demand`	Combined friction-ellipse demand sqrt((mu_x/Dx)^2 + (mu_y/Dy)^2)	–

Behavior

Source

dyad

"""
Joint for a wheel with Pacejka magic-formula friction rolling on the flat plane
y = 0.

The magic-formula counterpart of `SlipWheelJoint`. The contact patch can slip
and the resulting tangential force is computed from independent longitudinal and
lateral Pacejka magic-formula curves evaluated by
`MultibodyComponents.pacejka_magic` on the corresponding slip-velocity
components. Each direction has its own stiffness/shape/peak/curvature factors
(`Bx,Cx,Dx,Ex` longitudinal and `By,Cy,Dy,Ey` lateral). The friction curve is
smooth in all derivatives, so both explicit RK and BDF-type integrators work.

!!! warn "Normal force"
    The wheel cannot leave the ground. Make sure that the normal force `f_n`
    never becomes negative.
"""
component PacejkaWheelJoint
  extends Renderable(color = [1, 0, 0, 1])
  extends CutJoint()
  frame_a = Frame3D() {}
  """
   Optional road-surface frame: when surface=true the contact patch position
   equals this frame's position, so a connected component can define the road
   height y = f(x, z). When false, the wheel rolls on the flat plane y = 0.
  """
  surface_frame = Frame3D() if surface {}
  "Tire visualization: cylinder centered on frame_a with axis along z"
  tire_shape = CylinderShape(render = render, color = color, r = frame_a.r_0, R = transpose(frame_a.R), r_shape = [0, 0, -width / 2], length_direction = [0, 0, 1], length = width, width = 2 * radius, height = 2 * radius)
  "Rim cross-cylinder 1: thin spoke spanning the wheel diameter; rotates with the wheel to make the spin visible"
  rim1_shape = CylinderShape(render = render, color = [0.8, 0.8, 0.8, 1.0], r = frame_a.r_0, R = transpose(frame_a.R), r_shape = [-radius, 0, 0], length_direction = [1, 0, 0], length = 2 * radius, width = radius / 3, height = radius / 3)
  "Rim cross-cylinder 2: second spoke, perpendicular to the first"
  rim2_shape = CylinderShape(render = render, color = [0.8, 0.8, 0.8, 1.0], r = frame_a.r_0, R = transpose(frame_a.R), r_shape = [0, -radius, 0], length_direction = [0, 1, 0], length = 2 * radius, width = radius / 3, height = radius / 3)
  "Use an external road surface supplied via surface_frame (else flat ground y=0)"
  structural parameter surface::Boolean = false
  """
  Include the wheel's angular state (Cardan angles). Set false when the wheel is
   mounted on an already-rooted axis: the orientation is then taken from frame_a
   and no orientation loop/state is introduced (avoids redundant constraint forces
   in free-floating multi-wheel assemblies).
  """
  structural parameter angular_state::Boolean = true
  "Rotation axis sequence for the Cardan angles that orient `frame_a` (default y, z, x)"
  structural parameter sequence::Integer[3] = [2, 3, 1]
  """
  Couple longitudinal and lateral forces through the combined friction ellipse
   (semi-axes `Dx`·f_n, `Dy`·f_n). When true (default), grip is shared between
   longitudinal and lateral demand, so the resultant force cannot exceed the
   friction limit under simultaneous slip. When false, the two magic-formula
   curves act independently (the resultant can modestly exceed the per-axis peak).
  """
  structural parameter friction_ellipse::Boolean = true
  "Radius of the wheel"
  parameter radius::Length = 0.3
  "Width of the wheel in animations"
  parameter width::Real = 0.035
  "Longitudinal magic-formula stiffness factor (initial slope)"
  parameter Bx::Real = 10
  "Longitudinal magic-formula shape factor"
  parameter Cx::Real = 1.9
  "Longitudinal magic-formula peak friction coefficient"
  parameter Dx::Real = 1.0
  "Longitudinal magic-formula curvature factor"
  parameter Ex::Real = 0.97
  "Lateral magic-formula stiffness factor (initial slope)"
  parameter By::Real = 10
  "Lateral magic-formula shape factor"
  parameter Cy::Real = 1.3
  "Lateral magic-formula peak friction coefficient"
  parameter Dy::Real = 0.9
  "Lateral magic-formula curvature factor"
  parameter Ey::Real = 0.97
  "Small value added to `v_slip` to avoid division by zero in the slip model"
  parameter v_small::Real = 1e-5
  "x-position of the wheel axis"
  variable x::Length(statePriority = 15)
  "y-position of the wheel axis"
  variable y::Length(statePriority = 0)
  "z-position of the wheel axis"
  variable z::Length(statePriority = 15)
  "Angles that rotate the world frame into `frame_a` around the y-, z-, x-axis"
  variable angles::Angle(statePriority = 15)[3] if angular_state
  "Time derivatives of `angles`"
  variable der_angles::AngularVelocity(statePriority = 30)[3] if angular_state
  "Wheel rolling angle (`angles[2]`)"
  variable phi_roll::Angle if angular_state
  "Wheel rolling angular velocity"
  variable w_roll::AngularVelocity
  "Position vector from world frame to contact point on the road, resolved in world frame"
  variable r_road_0::Position[3]
  "Force vector on the wheel at the contact point, resolved in world frame"
  variable f_wheel_0::Real[3]
  "Contact force acting on the wheel in the normal direction"
  variable f_n::Dyad.Force
  "Contact force acting on the wheel in the lateral direction"
  variable f_lat::Dyad.Force
  "Contact force acting on the wheel in the longitudinal direction"
  variable f_long::Dyad.Force
  "Unit vector along the wheel axis, resolved in world frame"
  variable e_axis_0::Real[3]
  "Distance vector from wheel center to contact point, resolved in world frame"
  variable delta_0::Real[3]
  "Unit vector normal to the road at the contact point, resolved in world frame"
  variable e_n_0::Real[3]
  "Unit vector in the lateral direction of the road at the contact point, resolved in world frame"
  variable e_lat_0::Real[3]
  "Unit vector in the longitudinal direction of the road at the contact point, resolved in world frame"
  variable e_long_0::Real[3]
  "Road heading at (s, w), resolved in world frame (unit vector)"
  variable e_s_0::Real[3]
  "Road surface parameter 1"
  variable s::Real
  "Road surface parameter 2"
  variable w::Real
  "Velocity of wheel center, resolved in world frame"
  variable v_0::Velocity[3]
  "Angular velocity of the wheel, resolved in world frame"
  variable w_0::AngularVelocity[3]
  "Velocity of the contact point, resolved in world frame"
  variable vContact_0::Velocity[3]
  "Auxiliary vector for building the contact frame"
  variable aux::Real[3]
  "Slip velocity in the longitudinal direction"
  variable v_slip_long::Velocity
  "Slip velocity in the lateral direction"
  variable v_slip_lat::Velocity
  "Slip ratio = v_slip_long / (v_0 · e_long_0)"
  variable slip_ratio::Real
  "Longitudinal pure-slip friction coefficient from the magic formula"
  variable mu_x::Real
  "Lateral pure-slip friction coefficient from the magic formula"
  variable mu_y::Real
  "Combined friction-ellipse demand sqrt((mu_x/Dx)^2 + (mu_y/Dy)^2)"
  variable ellipse_demand::Real if friction_ellipse
relations
  # Solver guesses for algebraic intermediates (mirrors Multibody.jl defaults)
  guess x = 0
  guess y = radius
  guess z = 0
  if angular_state
    guess angles = [0, 0, 0]
    guess der_angles = [0, 0, 0]
    guess phi_roll = 0
  end
  guess w_roll = 0
  guess f_n = 1.0
  guess delta_0 = [0, -radius, 0]
  guess e_n_0 = [0, 1, 0]
  guess aux = [1, 0, 0]
  guess v_slip_long = 0
  guess v_slip_lat = 0
  # Road description
  if surface
    # The surface frame is a kinematic road reference: it exerts no force.
    surface_frame.f = [0, 0, 0]
    surface_frame.tau = [0, 0, 0]
    # Contact patch position is the surface-frame position; s, w (= x, z) remain
    # free and are solved by the contact constraints, while the connected road
    # component imposes the height (y of the frame position).
    s = surface_frame.r_0[1]
    w = surface_frame.r_0[3]
    r_road_0 = surface_frame.r_0
    # Contact tangent/normal come from the road frame's orientation: the connected
    # road component orients surface_frame so its y-axis is the surface normal and
    # its x-axis is the heading. This is exact for any surface the road describes
    # (a flat/excitation road keeps R = I → vertical normal).
    e_n_0 = resolve1(surface_frame.R, [0, 1, 0])
    e_s_0 = resolve1(surface_frame.R, [1, 0, 0])
  else
    # Flat road description
    r_road_0 = [s, 0, w]
    e_n_0 = [0, 1, 0]
    e_s_0 = [1, 0, 0]
  end
  # Pose / orientation.
  frame_a.r_0 = [x, y, z]
  Rleaf = frame_a.R
  if angular_state
    # Wheel carries its own Cardan-angle orientation state. It goes through
    # CutJoint: when iscut is true the orientation is imposed as a 3-residue cut
    # (opening the loop with an attached rooted axis) instead of full matrix
    # equality. Angular velocity is derived from der_angles (avoids state
    # selection picking D(R)-based states).
    Rroot = RR(axes_rotations(sequence, angles, der_angles))
    der_angles = der(angles)
    phi_roll = angles[2]
    w_roll = der(phi_roll)
    w_0 = resolve1(frame_a.R, angular_velocity2(axes_rotations(sequence, angles, der_angles)))
  else
    # No angular state: orientation (and its angular velocity) come from frame_a,
    # which is rooted by the axis the wheel is mounted on. No orientation loop is
    # introduced, so there is nothing to cut.
    Rroot = frame_a.R
    w_roll = dot(angular_velocity2(ori(frame_a)), [0, 0, 1])
    w_0 = angular_velocity1(ori(frame_a))
  end
  # Contact-point coordinate frame (e_long_0, e_lat_0, e_n_0)
  e_axis_0 = resolve1(frame_a.R, [0, 0, 1])
  assert(abs(dot(e_n_0, e_axis_0)) < 0.99, "Wheel lays nearly on the ground (which is a singularity)")
  aux = cross(e_n_0, e_axis_0)
  e_long_0 = aux / norm_(aux)
  e_lat_0 = -cross(e_long_0, e_n_0) # wheel rotation axis and lateral axis are opposite
  # Contact-point geometry
  delta_0 = r_road_0 - frame_a.r_0
  0 = dot(delta_0, e_axis_0)
  0 = dot(delta_0, e_long_0)
  # Holonomic constraint: wheel touches road, no penetration
  0 = radius - dot(delta_0, cross(e_long_0, e_axis_0))
  # Velocities
  v_0 = der(frame_a.r_0)
  vContact_0 = v_0 + cross(w_0, delta_0)
  # Slip kinematics
  v_slip_lat = dot(vContact_0, e_lat_0)
  v_slip_long = dot(vContact_0, e_long_0)
  # The denominator (longitudinal rolling speed) is legitimately zero for a
  # stationary wheel; regularize it (sign-preserving) so this purely-diagnostic
  # quantity cannot poison initialization, which evaluates all observed
  # equations at the initial guess.
  slip_ratio = v_slip_long / ifelse(abs(dot(v_0, e_long_0)) < v_small, v_small, dot(v_0, e_long_0))
  # Pure-slip friction coefficients from the independent magic-formula curves.
  # pacejka_magic is odd, so each carries the sign of its slip-velocity component
  # (matching SlipWheelJoint's force signs).
  mu_x = MultibodyComponents.pacejka_magic(Bx, Cx, Dx, Ex, v_slip_long)
  mu_y = MultibodyComponents.pacejka_magic(By, Cy, Dy, Ey, v_slip_lat)
  if friction_ellipse
    # Combined friction ellipse: scale both components by 1/max(1, demand) so the
    # resultant stays on/inside the ellipse with semi-axes (Dx, Dy)·f_n. The
    # scaling is inactive (demand <= 1) for pure longitudinal or pure lateral
    # slip, and only shares grip once both directions are loaded.
    ellipse_demand = sqrt((mu_x / Dx) ^ 2 + (mu_y / Dy) ^ 2)
    f_long = f_n * mu_x / max(1, ellipse_demand)
    f_lat = f_n * mu_y / max(1, ellipse_demand)
  else
    # Independent longitudinal/lateral curves (no cross-coupling).
    f_long = f_n * mu_x
    f_lat = f_n * mu_y
  end
  # Contact force in world frame
  f_wheel_0 = f_n * e_n_0 - f_lat * e_lat_0 - f_long * e_long_0
  # Force and torque balance at the wheel center
  [0, 0, 0] = frame_a.f + resolve2(frame_a.R, f_wheel_0)
  [0, 0, 0] = frame_a.tau + resolve2(frame_a.R, cross(delta_0, f_wheel_0))
metadata {
  "Dyad": {
    "icons": {"default": "dyad://MultibodyComponents/PacejkaWheelJoint.svg"},
    "labels": [
      {
        "label": "$(instance)",
        "x": 500,
        "y": 200,
        "rot": 0,
        "attrs": {"font-size": "160"}
      }
    ]
  }
}
end

Flattened Source

dyad

"""
Joint for a wheel with Pacejka magic-formula friction rolling on the flat plane
y = 0.

The magic-formula counterpart of `SlipWheelJoint`. The contact patch can slip
and the resulting tangential force is computed from independent longitudinal and
lateral Pacejka magic-formula curves evaluated by
`MultibodyComponents.pacejka_magic` on the corresponding slip-velocity
components. Each direction has its own stiffness/shape/peak/curvature factors
(`Bx,Cx,Dx,Ex` longitudinal and `By,Cy,Dy,Ey` lateral). The friction curve is
smooth in all derivatives, so both explicit RK and BDF-type integrators work.

!!! warn "Normal force"
    The wheel cannot leave the ground. Make sure that the normal force `f_n`
    never becomes negative.
"""
component PacejkaWheelJoint
  parameter render::Boolean = true
  parameter color::Real[4] = [0.5, 0.5, 0.5, 1.0]
  parameter specular_coefficient::Real = 1.5
  structural parameter iscut::Boolean = false
  final structural parameter residual::Real[3] = zeros(3) if iscut
  variable Rleaf::Real[3, 3]
  variable Rroot::Real[3, 3]
  frame_a = Frame3D() {}
  """
   Optional road-surface frame: when surface=true the contact patch position
   equals this frame's position, so a connected component can define the road
   height y = f(x, z). When false, the wheel rolls on the flat plane y = 0.
  """
  surface_frame = Frame3D() if surface {}
  "Tire visualization: cylinder centered on frame_a with axis along z"
  tire_shape = CylinderShape(render = render, color = color, r = frame_a.r_0, R = transpose(frame_a.R), r_shape = [0, 0, -width / 2], length_direction = [0, 0, 1], length = width, width = 2 * radius, height = 2 * radius)
  "Rim cross-cylinder 1: thin spoke spanning the wheel diameter; rotates with the wheel to make the spin visible"
  rim1_shape = CylinderShape(render = render, color = [0.8, 0.8, 0.8, 1.0], r = frame_a.r_0, R = transpose(frame_a.R), r_shape = [-radius, 0, 0], length_direction = [1, 0, 0], length = 2 * radius, width = radius / 3, height = radius / 3)
  "Rim cross-cylinder 2: second spoke, perpendicular to the first"
  rim2_shape = CylinderShape(render = render, color = [0.8, 0.8, 0.8, 1.0], r = frame_a.r_0, R = transpose(frame_a.R), r_shape = [0, -radius, 0], length_direction = [0, 1, 0], length = 2 * radius, width = radius / 3, height = radius / 3)
  "Use an external road surface supplied via surface_frame (else flat ground y=0)"
  structural parameter surface::Boolean = false
  """
  Include the wheel's angular state (Cardan angles). Set false when the wheel is
   mounted on an already-rooted axis: the orientation is then taken from frame_a
   and no orientation loop/state is introduced (avoids redundant constraint forces
   in free-floating multi-wheel assemblies).
  """
  structural parameter angular_state::Boolean = true
  "Rotation axis sequence for the Cardan angles that orient `frame_a` (default y, z, x)"
  structural parameter sequence::Integer[3] = [2, 3, 1]
  """
  Couple longitudinal and lateral forces through the combined friction ellipse
   (semi-axes `Dx`·f_n, `Dy`·f_n). When true (default), grip is shared between
   longitudinal and lateral demand, so the resultant force cannot exceed the
   friction limit under simultaneous slip. When false, the two magic-formula
   curves act independently (the resultant can modestly exceed the per-axis peak).
  """
  structural parameter friction_ellipse::Boolean = true
  "Radius of the wheel"
  parameter radius::Length = 0.3
  "Width of the wheel in animations"
  parameter width::Real = 0.035
  "Longitudinal magic-formula stiffness factor (initial slope)"
  parameter Bx::Real = 10
  "Longitudinal magic-formula shape factor"
  parameter Cx::Real = 1.9
  "Longitudinal magic-formula peak friction coefficient"
  parameter Dx::Real = 1.0
  "Longitudinal magic-formula curvature factor"
  parameter Ex::Real = 0.97
  "Lateral magic-formula stiffness factor (initial slope)"
  parameter By::Real = 10
  "Lateral magic-formula shape factor"
  parameter Cy::Real = 1.3
  "Lateral magic-formula peak friction coefficient"
  parameter Dy::Real = 0.9
  "Lateral magic-formula curvature factor"
  parameter Ey::Real = 0.97
  "Small value added to `v_slip` to avoid division by zero in the slip model"
  parameter v_small::Real = 1e-5
  "x-position of the wheel axis"
  variable x::Length(statePriority = 15)
  "y-position of the wheel axis"
  variable y::Length(statePriority = 0)
  "z-position of the wheel axis"
  variable z::Length(statePriority = 15)
  "Angles that rotate the world frame into `frame_a` around the y-, z-, x-axis"
  variable angles::Angle(statePriority = 15)[3] if angular_state
  "Time derivatives of `angles`"
  variable der_angles::AngularVelocity(statePriority = 30)[3] if angular_state
  "Wheel rolling angle (`angles[2]`)"
  variable phi_roll::Angle if angular_state
  "Wheel rolling angular velocity"
  variable w_roll::AngularVelocity
  "Position vector from world frame to contact point on the road, resolved in world frame"
  variable r_road_0::Position[3]
  "Force vector on the wheel at the contact point, resolved in world frame"
  variable f_wheel_0::Real[3]
  "Contact force acting on the wheel in the normal direction"
  variable f_n::Dyad.Force
  "Contact force acting on the wheel in the lateral direction"
  variable f_lat::Dyad.Force
  "Contact force acting on the wheel in the longitudinal direction"
  variable f_long::Dyad.Force
  "Unit vector along the wheel axis, resolved in world frame"
  variable e_axis_0::Real[3]
  "Distance vector from wheel center to contact point, resolved in world frame"
  variable delta_0::Real[3]
  "Unit vector normal to the road at the contact point, resolved in world frame"
  variable e_n_0::Real[3]
  "Unit vector in the lateral direction of the road at the contact point, resolved in world frame"
  variable e_lat_0::Real[3]
  "Unit vector in the longitudinal direction of the road at the contact point, resolved in world frame"
  variable e_long_0::Real[3]
  "Road heading at (s, w), resolved in world frame (unit vector)"
  variable e_s_0::Real[3]
  "Road surface parameter 1"
  variable s::Real
  "Road surface parameter 2"
  variable w::Real
  "Velocity of wheel center, resolved in world frame"
  variable v_0::Velocity[3]
  "Angular velocity of the wheel, resolved in world frame"
  variable w_0::AngularVelocity[3]
  "Velocity of the contact point, resolved in world frame"
  variable vContact_0::Velocity[3]
  "Auxiliary vector for building the contact frame"
  variable aux::Real[3]
  "Slip velocity in the longitudinal direction"
  variable v_slip_long::Velocity
  "Slip velocity in the lateral direction"
  variable v_slip_lat::Velocity
  "Slip ratio = v_slip_long / (v_0 · e_long_0)"
  variable slip_ratio::Real
  "Longitudinal pure-slip friction coefficient from the magic formula"
  variable mu_x::Real
  "Lateral pure-slip friction coefficient from the magic formula"
  variable mu_y::Real
  "Combined friction-ellipse demand sqrt((mu_x/Dx)^2 + (mu_y/Dy)^2)"
  variable ellipse_demand::Real if friction_ellipse
relations
  if iscut
    residue(Rleaf, Rroot) = residual
  else
    Rleaf = Rroot
  end
  # Solver guesses for algebraic intermediates (mirrors Multibody.jl defaults)
  guess x = 0
  guess y = radius
  guess z = 0
  if angular_state
    guess angles = [0, 0, 0]
    guess der_angles = [0, 0, 0]
    guess phi_roll = 0
  end
  guess w_roll = 0
  guess f_n = 1.0
  guess delta_0 = [0, -radius, 0]
  guess e_n_0 = [0, 1, 0]
  guess aux = [1, 0, 0]
  guess v_slip_long = 0
  guess v_slip_lat = 0
  # Road description
  if surface
    # The surface frame is a kinematic road reference: it exerts no force.
    surface_frame.f = [0, 0, 0]
    surface_frame.tau = [0, 0, 0]
    # Contact patch position is the surface-frame position; s, w (= x, z) remain
    # free and are solved by the contact constraints, while the connected road
    # component imposes the height (y of the frame position).
    s = surface_frame.r_0[1]
    w = surface_frame.r_0[3]
    r_road_0 = surface_frame.r_0
    # Contact tangent/normal come from the road frame's orientation: the connected
    # road component orients surface_frame so its y-axis is the surface normal and
    # its x-axis is the heading. This is exact for any surface the road describes
    # (a flat/excitation road keeps R = I → vertical normal).
    e_n_0 = resolve1(surface_frame.R, [0, 1, 0])
    e_s_0 = resolve1(surface_frame.R, [1, 0, 0])
  else
    # Flat road description
    r_road_0 = [s, 0, w]
    e_n_0 = [0, 1, 0]
    e_s_0 = [1, 0, 0]
  end
  # Pose / orientation.
  frame_a.r_0 = [x, y, z]
  Rleaf = frame_a.R
  if angular_state
    # Wheel carries its own Cardan-angle orientation state. It goes through
    # CutJoint: when iscut is true the orientation is imposed as a 3-residue cut
    # (opening the loop with an attached rooted axis) instead of full matrix
    # equality. Angular velocity is derived from der_angles (avoids state
    # selection picking D(R)-based states).
    Rroot = RR(axes_rotations(sequence, angles, der_angles))
    der_angles = der(angles)
    phi_roll = angles[2]
    w_roll = der(phi_roll)
    w_0 = resolve1(frame_a.R, angular_velocity2(axes_rotations(sequence, angles, der_angles)))
  else
    # No angular state: orientation (and its angular velocity) come from frame_a,
    # which is rooted by the axis the wheel is mounted on. No orientation loop is
    # introduced, so there is nothing to cut.
    Rroot = frame_a.R
    w_roll = dot(angular_velocity2(ori(frame_a)), [0, 0, 1])
    w_0 = angular_velocity1(ori(frame_a))
  end
  # Contact-point coordinate frame (e_long_0, e_lat_0, e_n_0)
  e_axis_0 = resolve1(frame_a.R, [0, 0, 1])
  assert(abs(dot(e_n_0, e_axis_0)) < 0.99, "Wheel lays nearly on the ground (which is a singularity)")
  aux = cross(e_n_0, e_axis_0)
  e_long_0 = aux / norm_(aux)
  e_lat_0 = -cross(e_long_0, e_n_0) # wheel rotation axis and lateral axis are opposite
  # Contact-point geometry
  delta_0 = r_road_0 - frame_a.r_0
  0 = dot(delta_0, e_axis_0)
  0 = dot(delta_0, e_long_0)
  # Holonomic constraint: wheel touches road, no penetration
  0 = radius - dot(delta_0, cross(e_long_0, e_axis_0))
  # Velocities
  v_0 = der(frame_a.r_0)
  vContact_0 = v_0 + cross(w_0, delta_0)
  # Slip kinematics
  v_slip_lat = dot(vContact_0, e_lat_0)
  v_slip_long = dot(vContact_0, e_long_0)
  # The denominator (longitudinal rolling speed) is legitimately zero for a
  # stationary wheel; regularize it (sign-preserving) so this purely-diagnostic
  # quantity cannot poison initialization, which evaluates all observed
  # equations at the initial guess.
  slip_ratio = v_slip_long / ifelse(abs(dot(v_0, e_long_0)) < v_small, v_small, dot(v_0, e_long_0))
  # Pure-slip friction coefficients from the independent magic-formula curves.
  # pacejka_magic is odd, so each carries the sign of its slip-velocity component
  # (matching SlipWheelJoint's force signs).
  mu_x = MultibodyComponents.pacejka_magic(Bx, Cx, Dx, Ex, v_slip_long)
  mu_y = MultibodyComponents.pacejka_magic(By, Cy, Dy, Ey, v_slip_lat)
  if friction_ellipse
    # Combined friction ellipse: scale both components by 1/max(1, demand) so the
    # resultant stays on/inside the ellipse with semi-axes (Dx, Dy)·f_n. The
    # scaling is inactive (demand <= 1) for pure longitudinal or pure lateral
    # slip, and only shares grip once both directions are loaded.
    ellipse_demand = sqrt((mu_x / Dx) ^ 2 + (mu_y / Dy) ^ 2)
    f_long = f_n * mu_x / max(1, ellipse_demand)
    f_lat = f_n * mu_y / max(1, ellipse_demand)
  else
    # Independent longitudinal/lateral curves (no cross-coupling).
    f_long = f_n * mu_x
    f_lat = f_n * mu_y
  end
  # Contact force in world frame
  f_wheel_0 = f_n * e_n_0 - f_lat * e_lat_0 - f_long * e_long_0
  # Force and torque balance at the wheel center
  [0, 0, 0] = frame_a.f + resolve2(frame_a.R, f_wheel_0)
  [0, 0, 0] = frame_a.tau + resolve2(frame_a.R, cross(delta_0, f_wheel_0))
metadata {
  "Dyad": {
    "icons": {"default": "dyad://MultibodyComponents/PacejkaWheelJoint.svg"},
    "labels": [
      {
        "label": "$(instance)",
        "x": 500,
        "y": 200,
        "rot": 0,
        "attrs": {"font-size": "160"}
      }
    ]
  }
}
end

Test Cases

No test cases defined.

Examples
Experiments
Analyses

Scalar Types

Partial Types

Enumerated Types

Function Types

Partial Types

Enumerated Types

Partial Types

Enumerated Types

Partial Types

Scalar Types

Partial Types

Partial Types

Types

Partial Types

Enumerated Types

Struct Types

Connectors

Components

Analyses

Examples

Tests

Scalar Types

Partial Types

Enumerated Types

PacejkaWheelJoint ​

Usage ​

Parameters: ​

Connectors ​

Variables ​

Behavior ​

Source ​

Test Cases ​

Related ​

`PacejkaWheelJoint`

Usage

Parameters:

Connectors

Variables

Behavior

Source

Test Cases

Related