Skip to content

Commit 8fc7bfa

Browse files
baggepinnenbaggepinnen
committed
add suspension example
1 parent 0ffe0e8 commit 8fc7bfa

File tree

2 files changed

+378
-0
lines changed

2 files changed

+378
-0
lines changed

‎docs/make.jl

Lines changed: 1 addition & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -39,6 +39,7 @@ makedocs(;
3939
"Bodies in space" => "examples/space.md",
4040
"Gyroscopic effects" => "examples/gyroscopic_effects.md",
4141
"Wheels" => "examples/wheel.md",
42+
"Suspension systems" => "examples/suspension.md",
4243
"Quadrotor with cable-suspended load" => "examples/quad.md",
4344
],
4445
"Rotations and orientation" => "rotations.md",

‎docs/src/examples/suspension.md

Lines changed: 377 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,377 @@
1+
# Quarter-car suspension
2+
3+
A simple suspension system for one wheel, roughly modeled after "Interaction between asymmetrical damping and geometrical nonlinearity in vehicle suspension systems improves comfort", Fernandes et al. 2020.
4+
5+
![suspension system](https://www.researchgate.net/publication/337953665/figure/fig1/AS:941829033828436@1601560948959/Schematic-representation-of-a-vehicle-suspension-system-with-Double-Wishbone-geometry-a_W640.jpg)
6+
7+
We model only the suspension, and excite the system with a force from the ground through a rod connected where the wheel would be. The actuation rod is modeled as a [`SphericalSpherical`](@ref) joint to avoid over constraining the system. We further add a prismatic joint called `body_upright` to hold the "car body" to the world frame. This allows the body to move up and down, but not sideways or rotate.
8+
9+
```@example suspension
10+
using Multibody
11+
using ModelingToolkit
12+
import ModelingToolkitStandardLibrary.Mechanical.TranslationalModelica as Translational
13+
using Plots
14+
using OrdinaryDiffEq
15+
using LinearAlgebra
16+
using JuliaSimCompiler
17+
using Test
18+
19+
t = Multibody.t
20+
D = Differential(t)
21+
W(args...; kwargs...) = Multibody.world
22+
23+
n = [1, 0, 0]
24+
AB = 146.5 / 1000
25+
BC = 233.84 / 1000
26+
CD = 228.60 / 1000
27+
DA = 221.43 / 1000
28+
BP = 129.03 / 1000
29+
DE = 310.31 / 1000
30+
t5 = 19.84 |> deg2rad
31+
32+
@mtkmodel QuarterCarSuspension begin
33+
@structural_parameters begin
34+
spring = true
35+
(jc = [0.5, 0.5, 0.5, 0.7])#, [description = "Joint color"]
36+
mirror = false
37+
end
38+
@parameters begin
39+
cs = 4000, [description = "Damping constant [Ns/m]"]
40+
ks = 44000, [description = "Spring constant [N/m]"]
41+
rod_radius = 0.02
42+
jr = 0.03, [description = "Radius of revolute joint"]
43+
end
44+
begin
45+
dir = mirror ? -1 : 1
46+
rRod1_ia = AB*normalize([0, -0.1, 0.2dir])
47+
rRod2_ib = BC*normalize([0, 0.2, 0dir])
48+
end
49+
@components begin
50+
r123 = JointRRR(n_a = n*dir, n_b = n*dir, rRod1_ia, rRod2_ib, rod_radius=0.018, rod_color=jc)
51+
r2 = Revolute(; n=n*dir, radius=jr, color=jc)
52+
b1 = FixedTranslation(radius = rod_radius, r = CD*normalize([0, -0.1, 0.3dir])) # CD
53+
chassis = FixedTranslation(r = DA*normalize([0, 0.2, 0.2*sin(t5)*dir]), render=false)
54+
chassis_frame = Frame()
55+
56+
if spring
57+
springdamper = SpringDamperParallel(c = ks, d = cs, s_unstretched = 1.3*BC, radius=rod_radius)
58+
end
59+
if spring
60+
spring_mount_F = FixedTranslation(r = 0.7*CD*normalize([0, -0.1, 0.3dir]), render=false)
61+
end
62+
if spring
63+
spring_mount_E = FixedTranslation(r = 1.3DA*normalize([0, 0.2, 0.2*sin(t5)*dir]), render=true)
64+
end
65+
end
66+
begin
67+
A = chassis.frame_b
68+
D = chassis.frame_a
69+
end
70+
@equations begin
71+
# Main loop
72+
connect(A, r123.frame_a)
73+
connect(r123.frame_b, b1.frame_b)
74+
connect(b1.frame_a, r2.frame_b)
75+
connect(r2.frame_a, D)
76+
77+
# Spring damper
78+
if spring
79+
connect(springdamper.frame_b, spring_mount_E.frame_b)
80+
connect(b1.frame_a, spring_mount_F.frame_a)
81+
connect(D, spring_mount_E.frame_a)
82+
connect(springdamper.frame_a, spring_mount_F.frame_b)
83+
end
84+
85+
connect(chassis_frame, chassis.frame_a)
86+
end
87+
end
88+
89+
mirror = false
90+
@mtkmodel SuspensionWithExcitation begin
91+
@structural_parameters begin
92+
mirror = false
93+
end
94+
@parameters begin
95+
ms = 1500, [description = "Mass of the car [kg]"]
96+
rod_radius = 0.02
97+
amplitude = 0.1, [description = "Amplitude of wheel displacement"]
98+
freq = 2, [description = "Frequency of wheel displacement"]
99+
end
100+
begin
101+
dir = mirror ? -1 : 1
102+
end
103+
@components begin
104+
world = W()
105+
chassis_frame = Frame()
106+
suspension = QuarterCarSuspension(; spring=true, mirror, rod_radius)
107+
108+
wheel_prismatic = Prismatic(n = [0,1,0], axisflange=true, state_priority=100, iscut=false)
109+
actuation_rod = SphericalSpherical(radius=rod_radius, r_0 = [0, BC, 0])
110+
actuation_position = FixedTranslation(r = [0, 0, CD*dir], render=false)
111+
wheel_position = Translational.Position(exact=true)
112+
113+
end
114+
@equations begin
115+
wheel_position.s_ref.u ~ amplitude*(sin(2pi*freq*t)) # Displacement of wheel
116+
connect(wheel_position.flange, wheel_prismatic.axis)
117+
118+
connect(world.frame_b, actuation_position.frame_a)
119+
connect(actuation_position.frame_b, wheel_prismatic.frame_a)
120+
connect(wheel_prismatic.frame_b, actuation_rod.frame_a,)
121+
connect(actuation_rod.frame_b, suspension.r123.frame_ib)
122+
123+
connect(chassis_frame, suspension.chassis_frame)
124+
end
125+
126+
end
127+
128+
@mtkmodel SuspensionWithExcitationAndMass begin
129+
@structural_parameters begin
130+
mirror = false
131+
end
132+
@parameters begin
133+
ms = 1500, [description = "Mass of the car [kg]"]
134+
rod_radius = 0.02
135+
end
136+
begin
137+
dir = mirror ? -1 : 1
138+
end
139+
@components begin
140+
world = W()
141+
mass = Body(m=ms, r_cm = 0.5DA*normalize([0, 0.2, 0.2*sin(t5)*dir]))
142+
excited_suspension = SuspensionWithExcitation(; suspension.spring=true, mirror, rod_radius)
143+
body_upright = Prismatic(n = [0, 1, 0], render = false, state_priority=1000)
144+
end
145+
@equations begin
146+
connect(world.frame_b, body_upright.frame_a)
147+
connect(body_upright.frame_b, excited_suspension.chassis_frame, mass.frame_a)
148+
end
149+
150+
end
151+
152+
@named model = SuspensionWithExcitationAndMass()
153+
model = complete(model)
154+
ssys = structural_simplify(IRSystem(model))
155+
156+
defs = [
157+
model.body_upright.s => 0.17
158+
model.excited_suspension.amplitude => 0.05
159+
model.excited_suspension.freq => 10
160+
model.excited_suspension.suspension.ks => 30*44000
161+
model.excited_suspension.suspension.cs => 30*4000
162+
model.ms => 1500/4
163+
model.excited_suspension.suspension.springdamper.num_windings => 10
164+
# model.r1.phi => -1.0889
165+
model.excited_suspension.suspension.r2.phi => -0.6031*(mirror ? -1 : 1)
166+
# model.r3.phi => 0.47595
167+
model.body_upright.v => 0.14
168+
]
169+
170+
display(sort(unknowns(ssys), by=string))
171+
172+
prob = ODEProblem(ssys, defs, (0, 2))
173+
# sol.u[1] = [0.17054946565059462, -0.6015077878288565, 0.274732819217001]
174+
sol = solve(prob, FBDF(autodiff=true))
175+
@test SciMLBase.successful_retcode(sol)
176+
rms(x) = sqrt(sum(abs2, x) / length(x))
177+
@test rms(sol(0:0.1:2, idxs=model.body_upright.v)) ≈ 2.740 atol=0.01
178+
```
179+
```@example suspension
180+
import GLMakie
181+
Multibody.render(model, sol, show_axis=false, x=-1, y=0.3, z=0.3, lookat=[0,0.3,0.3], timescale=3, filename="suspension.gif") # Video
182+
```
183+
184+
![suspension](suspension.gif)
185+
186+
187+
# Half-car suspension
188+
189+
In the example below, we extend the previous example to a half-car model with two wheels. We now model the car as a [`BodyShape`](@ref) and attach one suspension system on each side. This time, we let the car rotate around the x-axis to visualize roll effects due to uneven road surfaces. We excite each wheel with similar but slightly different frequencies to produce a beat.
190+
191+
```@example suspension
192+
@mtkmodel DoubleSuspensionWithExcitationAndMass begin
193+
@structural_parameters begin
194+
wheel_base = 1
195+
end
196+
@parameters begin
197+
ms = 1500, [description = "Mass of the car [kg]"]
198+
rod_radius = 0.02
199+
end
200+
@components begin
201+
world = W()
202+
mass = BodyShape(m=ms, r = [0,0,-wheel_base], radius=0.1, color=[0.4, 0.4, 0.4, 0.3])
203+
excited_suspension_r = SuspensionWithExcitation(; suspension.spring=true, mirror=false, rod_radius,
204+
actuation_position.r = [0, 0, (CD+wheel_base/2)],
205+
actuation_rod.r_0 = r_0 = [0, 0.1, 0],
206+
)
207+
excited_suspension_l = SuspensionWithExcitation(; suspension.spring=true, mirror=true, rod_radius,
208+
actuation_position.r = [0, 0, -(CD+wheel_base/2)],
209+
actuation_rod.r_0 = r_0 = [0, 0.1, 0],
210+
)
211+
body_upright = Prismatic(n = [0, 1, 0], render = false, state_priority=1000)
212+
body_upright2 = Revolute(n = [1, 0, 0], render = false, state_priority=1000, phi0=0, w0=0)
213+
# body_upright = Planar(n = [1, 0, 0], n_x = [0,0,1], render = false, state_priority=1000, radius=0.01)
214+
end
215+
@equations begin
216+
connect(world.frame_b, body_upright.frame_a)
217+
218+
connect(body_upright.frame_b, body_upright2.frame_a)
219+
connect(body_upright2.frame_b, mass.frame_cm)
220+
221+
# connect(body_upright.frame_b, mass.frame_cm)
222+
223+
connect(excited_suspension_r.chassis_frame, mass.frame_a)
224+
connect(excited_suspension_l.chassis_frame, mass.frame_b)
225+
end
226+
227+
end
228+
229+
@named model = DoubleSuspensionWithExcitationAndMass()
230+
model = complete(model)
231+
ssys = structural_simplify(IRSystem(model))
232+
233+
defs = [
234+
model.excited_suspension_r.amplitude => 0.05
235+
model.excited_suspension_r.freq => 10
236+
model.excited_suspension_r.suspension.ks => 30*44000
237+
model.excited_suspension_r.suspension.cs => 30*4000
238+
model.excited_suspension_r.suspension.springdamper.num_windings => 10
239+
model.excited_suspension_r.suspension.r2.phi => -0.6031*(1)
240+
241+
model.excited_suspension_l.amplitude => 0.05
242+
model.excited_suspension_l.freq => 9.5
243+
model.excited_suspension_l.suspension.ks => 30*44000
244+
model.excited_suspension_l.suspension.cs => 30*4000
245+
model.excited_suspension_l.suspension.springdamper.num_windings => 10
246+
model.excited_suspension_l.suspension.r2.phi => -0.6031*(+1)
247+
248+
model.ms => 1500/2
249+
250+
model.body_upright.s => 0.17
251+
model.body_upright.v => 0.14
252+
253+
# model.body_upright.prismatic_y.s => 0.17
254+
# model.body_upright.prismatic_y.v => 0.14
255+
]
256+
257+
display(sort(unknowns(ssys), by=string))
258+
259+
prob = ODEProblem(ssys, defs, (0, 4))
260+
sol = solve(prob, FBDF(autodiff=true), initializealg = ShampineCollocationInit())
261+
@test SciMLBase.successful_retcode(sol)
262+
# first(Multibody.render(model, sol, 0, show_axis=true, x=-1.5, y=0.0, z=0.0))
263+
```
264+
265+
```@example suspension
266+
Multibody.render(model, sol, show_axis=false, x=-1.5, y=0.3, z=0.0, lookat=[0,0.1,0.0], timescale=3, filename="suspension_halfcar.gif") # Video
267+
```
268+
269+
## Adding wheels
270+
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.
271+
The connection between the wheels and the ground form two kinematic loops together with the `body_upright` joint, we thus set both wheels to be cut joints using `iscut=true`.
272+
```@example suspension
273+
@mtkmodel ExcitedWheelAssembly begin
274+
@structural_parameters begin
275+
mirror = false
276+
end
277+
@parameters begin
278+
rod_radius = 0.02
279+
amplitude = 0.1, [description = "Amplitude of wheel displacement"]
280+
freq = 2, [description = "Frequency of wheel displacement"]
281+
end
282+
begin
283+
dir = mirror ? -1 : 1
284+
end
285+
@components begin
286+
chassis_frame = Frame()
287+
suspension = QuarterCarSuspension(; spring=true, mirror, rod_radius)
288+
289+
wheel = SlippingWheel(
290+
radius = 0.2,
291+
m = 15,
292+
I_axis = 0.06,
293+
I_long = 0.12,
294+
x0 = 0.0,
295+
z0 = 0.0,
296+
der_angles = [0, 0, 0],
297+
iscut = true, # NOTE: Only used since while we have an "upright joint"
298+
surface = (x,z)->amplitude*(sin(2pi*freq*t)), # Excitation from a time-varying surface profile
299+
)
300+
301+
end
302+
@equations begin
303+
connect(wheel.frame_a, suspension.r123.frame_ib)
304+
connect(chassis_frame, suspension.chassis_frame)
305+
end
306+
307+
end
308+
309+
@mtkmodel HalfCar begin
310+
@structural_parameters begin
311+
wheel_base = 1
312+
end
313+
@parameters begin
314+
ms = 1500, [description = "Mass of the car [kg]"]
315+
rod_radius = 0.02
316+
end
317+
@components begin
318+
world = W()
319+
mass = BodyShape(m=ms, r = [0,0,-wheel_base], radius=0.1, color=[0.4, 0.4, 0.4, 0.3])
320+
excited_suspension_r = ExcitedWheelAssembly(; mirror=false, rod_radius)
321+
excited_suspension_l = ExcitedWheelAssembly(; mirror=true, rod_radius)
322+
body_upright = Prismatic(n = [0, 1, 0], render = false, state_priority=2000, iscut=false)
323+
body_upright2 = Revolute(n = [1, 0, 0], render = false, state_priority=2000, phi0=0, w0=0, iscut=false)
324+
# body_upright = Planar(n = [1, 0, 0], n_x = [0,0,1], render = false, state_priority=100000, radius=0.01)
325+
end
326+
@equations begin
327+
connect(world.frame_b, body_upright.frame_a)
328+
connect(body_upright.frame_b, body_upright2.frame_a)
329+
connect(body_upright2.frame_b, mass.frame_cm)
330+
331+
# connect(body_upright.frame_b, mass.frame_cm)
332+
333+
connect(excited_suspension_r.chassis_frame, mass.frame_a)
334+
connect(excited_suspension_l.chassis_frame, mass.frame_b)
335+
end
336+
337+
end
338+
339+
@named model = HalfCar()
340+
model = complete(model)
341+
ssys = structural_simplify(IRSystem(model))
342+
343+
defs = [
344+
model.excited_suspension_r.amplitude => 0.05
345+
model.excited_suspension_r.freq => 10
346+
model.excited_suspension_r.suspension.ks => 30*44000
347+
model.excited_suspension_r.suspension.cs => 30*4000
348+
model.excited_suspension_r.suspension.springdamper.num_windings => 10
349+
model.excited_suspension_r.suspension.r2.phi => -0.6031*(1)
350+
351+
model.excited_suspension_l.amplitude => 0.05
352+
model.excited_suspension_l.freq => 9.5
353+
model.excited_suspension_l.suspension.ks => 30*44000
354+
model.excited_suspension_l.suspension.cs => 30*4000
355+
model.excited_suspension_l.suspension.springdamper.num_windings => 10
356+
model.excited_suspension_l.suspension.r2.phi => -0.6031*(+1)
357+
358+
model.ms => 1500
359+
360+
model.body_upright.s => 0.17
361+
model.body_upright.v => 0.14
362+
363+
# model.body_upright.prismatic_y.s => 0.17
364+
# model.body_upright.prismatic_y.v => 0.14
365+
366+
367+
# vec(ori(model.mass.frame_a).R .=> I(3))
368+
# vec(ori(model.excited_suspension_r.suspension.r123.jointUSR.frame_a).R .=> I(3))
369+
]
370+
371+
display(sort(unknowns(ssys), by=string))
372+
373+
prob = ODEProblem(ssys, defs, (0, 4))
374+
sol = solve(prob, FBDF(autodiff=false), initializealg = ShampineCollocationInit())#, u0 = prob.u0 .+ 1e-6 .* randn.())
375+
@show SciMLBase.successful_retcode(sol)
376+
first(Multibody.render(model, sol, 0, show_axis=true, x=-1.5, y=0.0, z=0.0))
377+
```

0 commit comments

Comments
 (0)