improve slip wheel docs

Author: baggepinnen

docs/src/examples/wheel.md

@@ -72,7 +72,35 @@ nothing # hide
 ![wheel animation](worldwheel.gif)
 
 ### Add slip
-The example below is similar to that above, but models a wheel with slip properties instead of ideal rolling.
+The example below is similar to that above, but models a wheel with slip properties instead of ideal rolling. We start by showing the slip model
+```@example WHEEL
+using Plots
+vAdhesion = 0.2
+vSlide = 0.4
+mu_A = 0.95
+mu_S = 0.7
+v = range(0, stop=1, length=500) # Simulating the slip velocity
+μ = Multibody.PlanarMechanics.limit_S_triple.(vAdhesion, vSlide, mu_A, mu_S, v)
+plot(v, μ, label=nothing, lw=2, color=:black, xlabel = "\$v_{Slip}\$", ylabel = "\$\\mu\$")
+scatter!([vAdhesion, vSlide], [mu_A, mu_S], color=:white, markerstrokecolor=:black)
+hline!([mu_A, mu_S], linestyle=:dash, color=:black, alpha=0.5)
+vline!([vAdhesion, vSlide], linestyle=:dash, color=:black, alpha=0.5)
+plot!(
+    xticks = ((vAdhesion, vSlide), ["\$v_{Adhesion}\$", "\$v_{Slide}\$"]),
+    yticks = ((mu_A, mu_S), ["\$\\mu_{adhesion}\$", "\$\\mu_{slide}\$"]),
+    framestyle = :zerolines,
+    legend = false,
+)
+```
+The longitudinal force on the tire is given by
+```math
+f_{long} = - f_n \dfrac{\mu(v_{Slip})}{v_{Slip}} v_{SlipLong}
+```
+where `f_n` is the normal force on the tire, `μ` is the friction coefficient from the slip model, and `v_{Slip}` is the magnitude of the slip velocity.
+
+The slip velocity is defined such that when the wheel is moving with positive velocity and increasing in speed (accelerating), the slip velocity is negative, i.e., the contact patch is moving slightly backwards. When the wheel is moving with positive velocity and decreasing in speed (braking), the slip velocity is positive, i.e., the contact patch is moving slightly forwards.
+
+
 ```@example WHEEL
 @mtkmodel SlipWheelInWorld begin
     @components begin
@@ -85,7 +113,6 @@ The example below is similar to that above, but models a wheel with slip propert
             x0 = 0.2,
             z0 = 0.2,
             der_angles = [0, 25, 0.1],
-
             mu_A = 0.95,             # Friction coefficient at adhesion
             mu_S = 0.5,             # Friction coefficient at sliding
             sAdhesion = 0.04,       # Adhesion slippage
@@ -334,13 +361,22 @@ prob = ODEProblem(ssys, [
     D(model.revolute.frame_b.phi) => 0,
     D(model.prismatic.r0[2]) => 0,
 ], (0.0, 15.0))
-sol = solve(prob, Rodas5Pr())
+sol = solve(prob, Rodas5P())
 render(model, sol, show_axis=false, x=0, y=0, z=4, traces=[model.slipBasedWheelJoint.frame_a], filename="slipwheel.gif", cache=false)
 nothing # hide
 ```
 
 ![slipwheel animation](slipwheel.gif)
 
+```@example WHEEL
+plot(sol, idxs=[
+    model.slipBasedWheelJoint.w_roll
+    model.slipBasedWheelJoint.v_long
+    model.slipBasedWheelJoint.v_slip_long
+    model.slipBasedWheelJoint.f_long
+], layout=4)
+```
+
 
 ## Planar two-track model
 A more elaborate example with 4 wheels.

src/wheels.jl

@@ -267,7 +267,7 @@ end
 """
     SlipWheelJoint(; name, radius, angles = zeros(3), der_angles = zeros(3), x0 = 0, y0 = radius, z0 = 0, sequence, iscut = false, surface = nothing, vAdhesion_min = 0.1, vSlide_min = 0.1, sAdhesion = 0.04, sSlide = 0.12, mu_A = 0.8, mu_S = 0.6, phi_roll = 0, w_roll = 0)
 
-Joint for a wheel with slip rolling on a surface.
+Joint for a wheel with slip rolling on a surface. See https://people.inf.ethz.ch/fcellier/MS/andres_ms.pdf for details.
 
 !!! tip "Integrator choice"
     The slip model contains a discontinuity in the second derivative at the transitions between adhesion and sliding. This can cause problems for integrators, in particular BDF-type integrators.
@@ -277,17 +277,39 @@ Joint for a wheel with slip rolling on a surface.
 
 # Parameters
 - `radius`: Radius of the wheel
-- `vAdhesion_min`: Minimum adhesion velocity
-- `vSlide_min`: Minimum sliding velocity
-- `sAdhesion`: Adhesion slippage
-- `sSlide`: Sliding slippage
+- `vAdhesion_min`: Minimum velocity for the peak of the adhesion curve (regularization close to 0)
+- `vSlide_min`: Minimum velocity for the start of the flat region of the slip curve (regularization close to 0)
+- `sAdhesion`: Adhesion slippage. The peak of the adhesion curve appears when the wheel slip is equal to `sAdhesion`.
+- `sSlide`: Sliding slippage. The flat region of the adhesion curve appears when the wheel slip is greater than `sSlide`.
 - `mu_A`: Friction coefficient at adhesion
 - `mu_S`: Friction coefficient at sliding
 - `surface`: By default, the wheel is rolling on a flat xz plane. A function `surface = (x, z)->y` may be provided to define a road surface. The function should return the height of the road at `(x, z)`. Note: if a function that depends on parameters is provided, make sure the parameters are scoped appropriately using, e.g., `ParentScope`.
 - `state`: (structural) whether or not the component has angular state variables. Default is `true`.
 
 # State and iscut
 When the wheel is mounted on an axis that is rooted, one may either supply `state=false` or `iscut = true`. With `state = false`, the angular state variables are not included in the wheel and there is thus no kinematic chain introduced. This reduces the total number of variables in the system. if the angular variables are required, one may instead pass `iscut=true` to cut the kinematic loop that is introduced when coupling the angles of the wheel to the orientation of the `frame_a`, unless this is cut elsewhere.
+
+# Understaning the slip model
+The following Julia code draws the slip model with descriptive labels
+```
+using Plots
+vAdhesion = 0.2
+vSlide = 0.4
+mu_A = 0.95
+mu_S = 0.7
+v = range(0, stop=1, length=500) # Simulating the slip velocity
+μ = Multibody.PlanarMechanics.limit_S_triple.(vAdhesion, vSlide, mu_A, mu_S, v)
+plot(v, μ, label=nothing, lw=2, color=:black, xlabel = "\$v_{Slip}\$", ylabel = "\$\\mu\$")
+scatter!([vAdhesion, vSlide], [mu_A, mu_S], color=:white, markerstrokecolor=:black)
+hline!([mu_A, mu_S], linestyle=:dash, color=:black, alpha=0.5)
+vline!([vAdhesion, vSlide], linestyle=:dash, color=:black, alpha=0.5)
+plot!(
+    xticks = ((vAdhesion, vSlide), ["\$v_{Adhesion}\$", "\$v_{Slide}\$"]),
+    yticks = ((mu_A, mu_S), ["\$\\mu_{adhesion}\$", "\$\\mu_{slide}\$"]),
+    framestyle = :zerolines,
+    legend = false,
+)
+```
 """
 @component function SlipWheelJoint(; name, radius, angles = zeros(3), der_angles=zeros(3), x0=0, y0 = radius, z0=0, sequence = [2, 3, 1], iscut=false, surface = nothing, vAdhesion_min = 0.05, vSlide_min = 0.15, sAdhesion = 0.04, sSlide = 0.12, mu_A = 0.8, mu_S = 0.6, phi_roll = 0, w_roll = 0, v_small = 1e-5, state=true)
     @parameters begin
@@ -296,7 +318,7 @@ When the wheel is mounted on an axis that is rooted, one may either supply `stat
         vSlide_min = vSlide_min, [description = "Minimum sliding velocity"]
         sAdhesion = sAdhesion, [description = "Adhesion slippage"]
         sSlide = sSlide, [description = "Sliding slippage"]
-        mu_A = mu_A, [description = "Friction coefficient at adhesion"]
+        mu_A = mu_A, [description = "Friction coefficient at adhesion peak"]
         mu_S = mu_S, [description = "Friction coefficient at sliding"]
         v_small = v_small, [description = "Small value added to v_slip to avoid division by zero in slip model."]
     end
@@ -362,9 +384,10 @@ When the wheel is mounted on an axis that is rooted, one may either supply `stat
         v_slip_long(t), [guess=0, description="Slip velocity in longitudinal direction"]
         v_slip_lat(t), [guess=0, description="Slip velocity in lateral direction"]
         v_slip(t), [description="Slip velocity, norm of component slip velocities"]
+        slip_ratio(t)
         f(t), [description="Total traction force"]
-        vAdhesion(t), [description="Adhesion velocity"]
-        vSlide(t), [description="Sliding velocity"]
+        vAdhesion(t), [description="Slip velocity at which adhesion is maximized"]
+        vSlide(t), [description="Slip velocity at which the flat region of the slip model starts"]
     end
 
     angles,der_angles,r_road_0,f_wheel_0,e_axis_0,delta_0,e_n_0,e_lat_0,e_long_0,e_s_0,v_0,w_0,vContact_0,aux = collect.((angles,der_angles,r_road_0,f_wheel_0,e_axis_0,delta_0,e_n_0,e_lat_0,e_long_0,e_s_0,v_0,w_0,vContact_0,aux))
@@ -441,15 +464,16 @@ When the wheel is mounted on an axis that is rooted, one may either supply `stat
                 v_slip ~ sqrt(v_slip_long^2 + v_slip_lat^2) + v_small
                 # -f_long * radius ~ flange_a.tau # No longer needed?
                 # frame_a.tau ~ 0
+                slip_ratio ~ v_slip_long / (v_0'e_long_0)
                 vAdhesion ~ max(vAdhesion_min, sAdhesion * abs(radius * w_roll))
                 vSlide ~ max(vSlide_min, sSlide * abs(radius * w_roll))
 
                 f ~ f_n * PlanarMechanics.limit_S_triple(vAdhesion, vSlide, mu_A, mu_S, v_slip) # limit_S_triple(x_max, x_sat, y_max, y_sat, x)
-                f_long ~ -f * v_slip_long / v_slip
-                f_lat ~ -f * v_slip_lat / v_slip
+                f_long ~ f * v_slip_long / v_slip
+                f_lat ~ f * v_slip_lat / v_slip
 
                 # Contact force (world frame)
-                f_wheel_0 .~ f_n * e_n_0 + f_lat * e_lat_0 + f_long * e_long_0
+                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
                 zeros(3) .~ collect(frame_a.f) + resolve2(Ra, f_wheel_0)