introduce stateless option in slip wheels (#150)

* introduce stateless option in slip wheels

* add damping to wheel rotation

* add some tests

* tweak tol

* move from self-hosted to ubuntu-latest since we made repo public

* update docs manifest

* rm explicit dev

* not ssh

Author: baggepinnen

.github/workflows/CI.yml

@@ -13,32 +13,30 @@ env:
 jobs:
   test:
     name: Julia ${{ matrix.version }} - ${{ matrix.os }} - ${{ matrix.arch }} - ${{ github.event_name }}
-    runs-on: self-hosted
+    runs-on: ubuntu-latest
     strategy:
       fail-fast: false
       matrix:
         version:
           - '1'
         os:
-          - self-hosted
+          - ubuntu-latest
         arch:
           - x64
     steps:
       - name: Set debug env
         run: export JULIA_DEBUG="loading"
       - uses: actions/[email protected]
+      - uses: julia-actions/cache@v1
       - uses: julia-actions/setup-julia@v1
         with:
           version: ${{ matrix.version }}
           arch: ${{ matrix.arch }}
-      - uses: webfactory/[email protected] 
-        with: 
-          ssh-private-key: ${{ secrets.JULIASIM_REGISTRY_SSH_KEY }}
       - uses: PumasAI/add-private-registry@main
         with:
           juliahub_token_encoded: ${{ secrets.JULIAHUB_TOKEN_ENCODED }}
-          private_registry_name: JuliaSimRegistry
-          private_registry_uuid: 309a7822-a73e-4490-9504-7d1983f27685
+          private_registry_name: JuliaHubRegistry
+          private_registry_uuid: de52bcdf-fcb2-40cf-a397-3d64b64f4d9c
       - uses: julia-actions/julia-buildpkg@v1
       - uses: julia-actions/julia-runtest@v1
       - uses: julia-actions/julia-processcoverage@v1

.github/workflows/CompatHelper.yml

@@ -5,7 +5,7 @@ on:
   workflow_dispatch:
 jobs:
   CompatHelper:
-    runs-on: self-hosted
+    runs-on: ubuntu-latest
     steps:
       - name: Pkg.add("CompatHelper")
         run: julia -e 'using Pkg; Pkg.add("CompatHelper")'

.github/workflows/Documentation.yml

@@ -17,24 +17,21 @@ jobs:
     permissions:
       contents: write
       statuses: write
-    runs-on: self-hosted
+    runs-on: ubuntu-latest
     steps:
       - uses: actions/[email protected]
+      - uses: julia-actions/cache@v1
       - uses: julia-actions/setup-julia@latest
         with:
           version: '1'
-      - uses: webfactory/[email protected]
-        with:
-          ssh-private-key: ${{ secrets.JULIASIM_REGISTRY_SSH_KEY }}
       - uses: PumasAI/add-private-registry@main
         with:
           juliahub_token_encoded: ${{ secrets.JULIAHUB_TOKEN_ENCODED }}
-          private_registry_name: JuliaSimRegistry
-          private_registry_uuid: 309a7822-a73e-4490-9504-7d1983f27685
+          private_registry_name: JuliaHubRegistry
+          private_registry_uuid: de52bcdf-fcb2-40cf-a397-3d64b64f4d9c
       - name: Install dependencies
         run: DISPLAY=:0 xvfb-run --auto-servernum -s '-screen 0 1024x768x24' julia --project=docs -e 'using Pkg;
-                                      Pkg.setprotocol!(; domain = "github.com", protocol = "ssh"); 
-                                      Pkg.add(PackageSpec(name = "JuliaSimCompiler", url = "https://github.com/JuliaComputing/JuliaSimCompiler.jl", rev = "master"));
+                                      Pkg.add(PackageSpec(name = "ModelingToolkitStandardLibrary", url = "https://github.com/SciML/ModelingToolkitStandardLibrary.jl", rev = "main"));
                                       Pkg.develop(PackageSpec(path=pwd()));
                                       Pkg.instantiate()'
       - name: Build and deploy

.github/workflows/TagBot.yml

@@ -12,7 +12,7 @@ jobs:
         github.actor == 'ChrisRackauckas' &&
         startsWith(github.event.comment.body, 'Triggering TagBot')
       )
-    runs-on: self-hosted
+    runs-on: ubuntu-latest
     steps:
       - uses: JuliaRegistries/TagBot@v1
         with:

docs/Manifest.toml

@@ -183,6 +183,8 @@ nothing # hide
 
 ![suspension](suspension.gif)
 
+Due to the high excitation frequency, we make use of the argument `timescale = 3` when we render the 3D animation to slow down the animation by a factor of 3.
+
 
 # Half-car suspension
 
@@ -235,13 +237,13 @@ defs = [
     model.excited_suspension_r.freq => 10
     model.excited_suspension_r.suspension.ks => 30*44000
     model.excited_suspension_r.suspension.cs => 30*4000
-    model.excited_suspension_r.suspension.r2.phi => -0.6031*(1)
+    model.excited_suspension_r.suspension.r2.phi => -0.6031
 
     model.excited_suspension_l.amplitude => 0.05
     model.excited_suspension_l.freq => 9.5
     model.excited_suspension_l.suspension.ks => 30*44000
     model.excited_suspension_l.suspension.cs => 30*4000
-    model.excited_suspension_l.suspension.r2.phi => -0.6031*(+1)
+    model.excited_suspension_l.suspension.r2.phi => -0.6031
 
     model.ms => 1500/2
 
@@ -269,10 +271,15 @@ nothing # hide
 
 ## Adding wheels
 The example below further extends the example from above by adding wheels to the suspension system. The excitation is not modeled as a time-varying surface profile, provided through the `surface` argument to the [`SlippingWheel`](@ref) component.
-The connection between the wheels and the ground form two kinematic loops together with the `body_upright` joint, we thus set all wheels to be cut joints using `iscut=true`. We start by adding a wheel to the quarter-car setup and then to the half-car setup.
+The connection between the wheels and the ground form two kinematic loops together with the `body_upright` joint, which we can handle in one of two ways.
+1. include a cut joint in the loop, e.g., by setting all wheels to be cut joints using `iscut=true`.
+2. Remove angular state variables from the wheels by passing `state = false`. Since we do not need the angular state variables in this case, we choose this option which reduces the total number of unknowns to solve for.
+
+We start by adding a wheel to the quarter-car setup and then to the half-car setup.
 
 ### Quarter car
 ```@example suspension
+using ModelingToolkitStandardLibrary.Mechanical.Rotational
 @mtkmodel ExcitedWheelAssembly begin
     @structural_parameters begin
         mirror = false
@@ -289,24 +296,28 @@ The connection between the wheels and the ground form two kinematic loops togeth
     @components begin
         chassis_frame = Frame()
         suspension = QuarterCarSuspension(; spring=true, mirror, rod_radius)
-
+        wheel_rotation = Revolute(n = [0, 0, 1], axisflange=true) # Wheel rotation axis
+        rotational_losses = Rotational.Damper(d = 0.1)
         wheel = SlippingWheel(
             radius = 0.2,
             m = 15,
-            I_axis = 0.06,
-            I_long = 0.12,
+            I_axis = 0.8,
+            I_long = 0.8,
             x0 = 0.0,
             z0 = 0.0,
-            der_angles = [0, 0, 0],
-            iscut = iscut,
+            state = false,
             # Note the ParentScope qualifier, without this, the parameters are treated as belonging to the wheel.wheel_joint component instead of the ExcitedWheelAssembly
             surface = (x,z)->ParentScope(ParentScope(amplitude))*(sin(2pi*ParentScope(ParentScope(freq))*t)), # Excitation from a time-varying surface profile
         )
 
     end
     @equations begin
-        connect(wheel.frame_a, suspension.r123.frame_ib)
+        connect(wheel.frame_a, wheel_rotation.frame_b)
+        connect(wheel_rotation.frame_a, suspension.r123.frame_ib)
         connect(chassis_frame, suspension.chassis_frame)
+
+        connect(rotational_losses.flange_a, wheel_rotation.axis)
+        connect(rotational_losses.flange_b, wheel_rotation.support)
     end
 end
 
@@ -339,7 +350,7 @@ defs = [
     model.excited_suspension.freq => 10
     model.excited_suspension.suspension.ks => 30*44000
     model.excited_suspension.suspension.cs => 30*4000
-    model.excited_suspension.suspension.r2.phi => -0.6031*(1)
+    model.excited_suspension.suspension.r2.phi => -0.6031
     model.body_upright.s => 0.17
     model.body_upright.v => 0.14
 ]
@@ -392,13 +403,13 @@ defs = [
     model.excited_suspension_r.freq => 10
     model.excited_suspension_r.suspension.ks => 30*44000
     model.excited_suspension_r.suspension.cs => 30*4000
-    model.excited_suspension_r.suspension.r2.phi => -0.6031*(1)
+    model.excited_suspension_r.suspension.r2.phi => -0.6031
 
     model.excited_suspension_l.amplitude => 0.05
     model.excited_suspension_l.freq => 9.5
     model.excited_suspension_l.suspension.ks => 30*44000
     model.excited_suspension_l.suspension.cs => 30*4000
-    model.excited_suspension_l.suspension.r2.phi => -0.6031*(+1)
+    model.excited_suspension_l.suspension.r2.phi => -0.6031
 
     model.ms => 1500
 
@@ -478,29 +489,33 @@ model = complete(model)
 @time "simplification" ssys = structural_simplify(IRSystem(model))
 
 defs = [
+    model.excited_suspension_br.wheel.wheeljoint.v_small => 10
     model.excited_suspension_br.amplitude => 0.02
     model.excited_suspension_br.freq => 10
     model.excited_suspension_br.suspension.ks => 30*44000
     model.excited_suspension_br.suspension.cs => 30*4000
-    model.excited_suspension_br.suspension.r2.phi => -0.6031*(1)
+    model.excited_suspension_br.suspension.r2.phi => -0.6031
 
+    model.excited_suspension_bl.wheel.wheeljoint.v_small => 10
     model.excited_suspension_bl.amplitude => 0.02
     model.excited_suspension_bl.freq => 10.5
     model.excited_suspension_bl.suspension.ks => 30*44000
     model.excited_suspension_bl.suspension.cs => 30*4000
-    model.excited_suspension_bl.suspension.r2.phi => -0.6031*(+1)
+    model.excited_suspension_bl.suspension.r2.phi => -0.6031
 
+    model.excited_suspension_fr.wheel.wheeljoint.v_small => 10
     model.excited_suspension_fr.amplitude => 0.02
     model.excited_suspension_fr.freq => 10
     model.excited_suspension_fr.suspension.ks => 30*44000
     model.excited_suspension_fr.suspension.cs => 30*4000
-    model.excited_suspension_fr.suspension.r2.phi => -0.6031*(1)
+    model.excited_suspension_fr.suspension.r2.phi => -0.6031
 
+    model.excited_suspension_fl.wheel.wheeljoint.v_small => 10
     model.excited_suspension_fl.amplitude => 0.02
     model.excited_suspension_fl.freq => 9.7
     model.excited_suspension_fl.suspension.ks => 30*44000
     model.excited_suspension_fl.suspension.cs => 30*4000
-    model.excited_suspension_fl.suspension.r2.phi => -0.6031*(+1)
+    model.excited_suspension_fl.suspension.r2.phi => -0.6031
 
     model.ms => 100
 
@@ -516,12 +531,12 @@ defs = [
 display(sort(unknowns(ssys), by=string))
 
 prob = ODEProblem(ssys, defs, (0, 3))
-sol = solve(prob, Rodas5P(autodiff=false), initializealg = BrownFullBasicInit())
+sol = solve(prob, Rodas5P(autodiff=false), initializealg = BrownFullBasicInit(), u0 = prob.u0.+1e-6.*randn.())
 @test SciMLBase.successful_retcode(sol)
 import GLMakie
-Multibody.render(model, sol, show_axis=false, x=-3.5, y=0.5, z=0.15, lookat=[0,0.1,0.0], timescale=2, filename="suspension_fullcar_wheels.gif") # Video
+Multibody.render(model, sol, show_axis=false, x=-3.5, y=0.5, z=0.15, lookat=[0,0.1,0.0], timescale=3, filename="suspension_fullcar_wheels.gif") # Video
+nothing # hide
 ```
 
 ![suspension with 4 wheels](suspension_fullcar_wheels.gif)
-```
 

docs/src/examples/suspension.md

@@ -281,8 +281,12 @@ Joint for a wheel with slip rolling on a surface.
 - `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.
 """
-@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.1, vSlide_min = 0.1, sAdhesion = 0.1, sSlide = 0.1, mu_A = 0.8, mu_S = 0.6, phi_roll = 0, w_roll = 0)
+@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.1, vSlide_min = 0.1, sAdhesion = 0.1, sSlide = 0.1, mu_A = 0.8, mu_S = 0.6, phi_roll = 0, w_roll = 0, v_small = 1e-5, state=true)
     @parameters begin
         radius = radius, [description = "Radius of the wheel"]
         vAdhesion_min = vAdhesion_min, [description = "Minimum adhesion velocity"]
@@ -291,6 +295,7 @@ Joint for a wheel with slip rolling on a surface.
         sSlide = sSlide, [description = "Sliding slippage"]
         mu_A = mu_A, [description = "Friction coefficient at adhesion"]
         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
     @variables begin
         (x(t) = x0), [state_priority = 15, description = "x-position of the wheel axis"]
@@ -361,8 +366,8 @@ Joint for a wheel with slip rolling on a surface.
 
     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))
 
-    @named frame_a = Frame(varw=true)
-    Ra = ori(frame_a, true)
+    @named frame_a = Frame(varw=state)
+    Ra = ori(frame_a, state)
 
     Rarot = axes_rotations(sequence, angles, -der_angles) # The - is the neg_w change
 
@@ -386,18 +391,25 @@ Joint for a wheel with slip rolling on a surface.
 
     equations = [
                 equations;
-                connect_orientation(Ra, Rarot; iscut)   # Ra ~ Rarot
-                Ra.w ~ Rarot.w
 
-                phi_roll ~ angles[2]
-                w_roll ~ D(phi_roll)
+                if state
+                    [
+                        connect_orientation(Ra, Rarot; iscut)
+                        Ra.w ~ Rarot.w
+                        phi_roll ~ angles[2]
+                        w_roll ~ D(phi_roll)
+                        der_angles .~ D.(angles)
+
+                    ]
+                else
+                    w_roll ~ -Ra.w[3] # w: Absolute angular velocity of local frame, resolved in local frame
+                end
 
                 # frame_a.R is computed from generalized coordinates
                 collect(frame_a.r_0) .~ [x, y, z]
-                der_angles .~ D.(angles)
 
 
-                # Coordinate system at contact point (e_long_0, e_lat_0, e_n_0)
+                # Coordinate system at contact point (e_long_0, e_lat_0, e_n_0), resolved on world frame
                 e_axis_0 .~ resolve1(Ra, [0, 0, 1])
                 aux .~ (cross(e_n_0, e_axis_0))
                 e_long_0 .~ (aux ./ _norm(aux))
@@ -411,19 +423,19 @@ Joint for a wheel with slip rolling on a surface.
                 # One holonomic positional constraint equation (no penetration in to the ground)
                 0 ~ radius - delta_0'cross(e_long_0, e_axis_0)
 
-                # Slip velocities
+                # Slip velocities (world frame)
                 v_0 .~ D.(frame_a.r_0)
                 w_0 .~ angular_velocity1(Ra)
                 vContact_0 .~ v_0 + cross(w_0, delta_0)
 
-                # Contact dynamics =============================================
+                # Contact dynamics (world frame) ===============================
 
                 v_slip_lat ~ vContact_0' * e_lat_0
                 v_slip_long ~ vContact_0' * e_long_0
                 # v_slip_lat ~ v_lat - 0
                 # v_slip_long ~ v_long - radius * w_roll
 
-                v_slip ~ sqrt(v_slip_long^2 + v_slip_lat^2) + 0.00001
+                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
                 vAdhesion ~ max(vAdhesion_min, sAdhesion * abs(radius * w_roll))
@@ -433,14 +445,16 @@ Joint for a wheel with slip rolling on a surface.
                 f_long ~ -f * v_slip_long / v_slip
                 f_lat ~ -f * v_slip_lat / v_slip
 
-                # Contact force
+                # Contact force (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
                 zeros(3) .~ collect(frame_a.f) + resolve2(Ra, f_wheel_0)
                 zeros(3) .~ collect(frame_a.tau) +
                             resolve2(Ra, cross(delta_0, f_wheel_0))]
 
+                    
+
                 # continuous_events = [
                 #     v_slip~vAdhesion
                 #     v_slip~vSlide
@@ -465,6 +479,7 @@ Wheel with slip rolling on a surface.
 - `width`: Width of the wheel (for rendering)
 - `x0`: Initial x-position of the wheel axis
 - `z0`: Initial z-position of the wheel axis
+- `state`: (structural) whether or not the component has angular state variables. 
 
 # Variables
 - `x`: x-position of the wheel axis
@@ -479,8 +494,8 @@ Wheel with slip rolling on a surface.
 See [Docs: Wheels](https://help.juliahub.com/multibody/dev/examples/wheel/)
 """
 @component function SlippingWheel(; name, radius, m, I_axis, I_long, width = 0.035, x0=0, z0=0,
-                      angles = zeros(3), der_angles = zeros(3), kwargs...)
-    @named wheeljoint = SlipWheelJoint(; radius, angles, x0, z0, der_angles, kwargs...)
+                      angles = zeros(3), der_angles = zeros(3), state = true, kwargs...)
+    @named wheeljoint = SlipWheelJoint(; radius, angles, x0, z0, der_angles, state, kwargs...)
     @named begin
         frame_a = Frame()
         body = Body(r_cm = [0, 0, 0],
@@ -512,10 +527,13 @@ See [Docs: Wheels](https://help.juliahub.com/multibody/dev/examples/wheel/)
 
     equations = Equation[wheeljoint.x ~ x
                          wheeljoint.z ~ z
-                         collect(wheeljoint.angles) .~ collect(angles)
-                         collect(wheeljoint.der_angles) .~ collect(der_angles)
                          connect(body.frame_a, frame_a)
                          connect(wheeljoint.frame_a, frame_a)]
+    if state
+        append!(equations, [collect(wheeljoint.angles) .~ collect(angles)
+                            collect(wheeljoint.der_angles) .~ collect(der_angles)]
+        )
+    end
     compose(ODESystem(equations, t; name), frame_a, wheeljoint, body)
 end