DeStabieleFiets_06-04-2014 - TU Delft Institutional Repository

Technische Universiteit Delft
Faculteit Elektrotechniek, Wiskunde en Informatica
Delft Institute of Applied Mathematics
De Stabiele Fiets
Scriptie ten behoeve van het
Delft Institute of Applied Mathematics
als onderdeel ter verkrijging
van de graad van
BACHELOR OF SCIENCE
in
TECHNISCHE WISKUNDE
door
KARIN MIRANDA KETELLAPPER
April 2014, Delft
BSc verslag TECHNISCHE WISKUNDE
“De Stabiele Fiets”
KARIN KETELLAPPER
Technische Universiteit Delft
Begeleider
Dr. B.J. Meulenbroek
Overige commissieleden
Drs. E.M. van Elderen
Prof.dr.ir. C. Vuik
April 2014, Delft
Inhoudsopgave
1 Voorwoord
3
2 Inleiding
5
3 Van co¨
efficienten naar model
3.1 Fysische betekenis van alle parameters . . . . .
3.1.1 De fietsparameters . . . . . . . . . . . .
3.1.2 De bewegingsparameters . . . . . . . . .
3.1.3 Welke parameters zijn constant en welke
3.2 Van co¨efficienten naar energie . . . . . . . . . .
3.3 Van energie naar Lagrange vergelijkingen . . .
3.4 Euler-Lagrange . . . . . . . . . . . . . . . . . .
3.4.1 Uitrekenen van termen met . . . . . .
3.5 Overige Lagrangiaanse vergelijkingen . . . . . .
3.6 Bewegingsvergelijkingen opstellen . . . . . . . .
3.7 Bewegingsvergelijkingen vereenvoudigen . . . .
.
.
.
.
.
.
.
.
.
.
.
6
6
6
8
9
9
11
13
15
17
18
19
4 Structuur onderzoek stabiliteit
4.1 Karakteristieke vergelijking . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
4.1.1 Oplossingen karakteristieke vergelijking . . . . . . . . . . . . . . . . . . . . . . . .
4.1.2 Routh Stability Criteria . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
21
21
22
22
5 Waardes co¨
efficienten invullen
23
6 Resultaten
6.1 Eigenwaarden & Stabiliteit . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
6.2 Toepassen Routh Stability Criteria . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
23
24
25
7 Een eerste conclusie
26
8 Verandering model
8.1 Resultaten na verandering . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
8.2 Conclusie . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
27
28
29
9 Gyroscopisch effect
9.1 Wat is het gyroscopisch effect . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
9.2 Gyroscopisch effect fiets . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
30
30
31
10 Eindconclusie
33
Referenties
35
A Parameters fiets
37
B Betekenis co¨
efficienten
38
C Maple
C.1 Karakteristieke vergelijking berekenen . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
C.2 Stabiliteitsgebied Routh Stability Criteria . . . . . . . . . . . . . . . . . . . . . . . . . . .
C.3 Lagrangiaanse vergelijking voor . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
39
39
39
40
D Matlab code
D.1 Startmodel . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
D.2 M varieren . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
D.3 hf varieren . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
41
41
42
43
E Routh determinant X
45
1
. . . . .
. . . . .
. . . . .
variabel
. . . . .
. . . . .
. . . . .
. . . . .
. . . . .
. . . . .
. . . . .
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
2
1
Voorwoord
Met deze scriptie sluit ik mijn Bachelor opleiding Technische Wiskunde aan de Faculteit EWI van de
Technische Universiteit in Delft af. De afgelopen jaren heb ik vele uren op mijn fiets doorgebracht
waarmee ik van huis naar EWI ging en weer terug. Daarom was het onderwerp voor mijn scriptie ook
snel gekozen. Met deze scriptie kon ik mijn hobby, fietsen, combineren met ´e´en van mijn favoriete vakken,
namelijk Modelleren.
Dit voorwoord wil ik graag gebruiken om een dankwoord uit te spreken naar een aantal mensen. Allereerst
aan Bernard Meulenbroek, die mij begeleid heeft bij het schrijven van mijn scriptie. Daarnaast wil ik
Hendrien Kakebeeke bedanken die heel mijn scriptie gelezen heeft en voorzien heeft van feedback. Als
laatste, zeker niet het minste, wil ik mijn ouders, Henk Ketellapper en Ada Hansum-van den Boomgaard
bedanken. Zij hebben het voor mij mogelijk gemaakt om deze studie te doen en mij altijd gesteund. Daar
zal ik ze altijd dankbaar voor zijn!
Delft, 7 april 2014
3
4
2
Inleiding
Bijna elke Nederlander heeft wel een fiets. Nederland zonder fietsen is bijna ondenkbaar. Het is een
handig vervoersmiddel voor in de stad, je kunt snel overal tussendoor, je kan de fiets makkelijk parkeren,
ideaal dus om van A naar B te komen. De vraag die veel mensen al gesteld hebben - Hoe komt het dat
de fiets overeind blijft? - kwam ook bij mij naar boven, want een stilstaande fiets valt gewoon om, en
als je heel langzaam gaat is het ook moeilijk om je evenwicht te bewaren. De vraag die in deze scriptie
centraal zal staan is: wanneer en waardoor blijft de fiets fietsen? Stel de fiets gaat rechtdoor, dan zijn
het voornamelijk de zijwaartse bewegingen en het sturen die een fiets uit balans kunnen brengen.
Uiteindelijk willen we een model voor de fiets maken die iets zegt over de stabiliteit bij een bepaalde
snelheid. Hiervoor hebben we de bewegingsvergelijkingen van de fiets nodig. Deze vergelijkingen stellen
we op met behulp van de Lagrangiaan, die weer is opgebouwd uit de kinetische en potentiele energie. Voor
het opstellen van die energie¨en hebben we de fietsparameters nodig. Dat zijn alle parameters die specifiek
voor die fiets zijn en ook hoeken die beschrijven wat de positie en stand van de fiets is. De snelheid van
de fiets wordt aangegeven met de parameter v. De hoeken δ en θ, die respectievelijk aangeven hoe schuin
het stuur staat en de afwijking van het achterwiel van het rechte pad, veronderstellen we heel klein en
worden enkel lineair meegenomen.
We komen dan op twee bewegingsvergelijkingen die samen een stelsel van twee 2e orde differentiaalvergelijkingen vormen. De parameters van die differentiaalvergelijkingen zijn en β, waarbij is hoe schuin
het achterwiel staat en β de grootte van de stuurhoek is. De karakteristieke vergelijking van dit stelsel
is een 4e graads vergelijking. Door deze op te lossen berekenen we de eigenwaarden. De fiets is stabiel
als alle re¨ele delen van de eigenwaarden negatief zijn. De eigenwaarden zijn afhankelijk van de snelheid
van de fiets. Door een figuur te maken waar we de snelheid tegenover de eigenwaarden plotten kunnen
we het stabiliteitsgebied uit de figuur aflezen.
Vervolgens kijken we naar het Routh Stability Criteria voor 4e graadsvergelijkingen. Wanneer alle
co¨efficienten uit de karakteristieke vergelijking en de Routhe determinant hetzelfde teken hebben, dan
is de fiets stabiel. Ook hier hangen die co¨efficienten af van de snelheid. We plotten weer de snelheid
tegenover de bijbehorende co¨efficienten en dan kunnen we het stabiliteitsgebied aflezen.
Wanneer we hebben gekeken naar het stabiliteitsgebied, dan kunnen we onderzoeken wat het gyroscopisch
effect voor invloed heeft op de stabiliteit van de fiets. Algemeen wordt aangenomen dat het gyroscopisch
effect noodzakelijk is, maar er is een artikel, zoals [2], waarin beweerd wordt dat het niet noodzakelijk is.
Uiteindelijk zullen we er op uit komen dat voor deze fiets, uit [1] die we in deze scriptie beschouwen, het
gyroscopisch effect er voor zorgt dat er een stabiliteitsgebied is. Zonder het gyroscopisch effect is deze
fiets instabiel voor elke snelheid.
5
3
Van co¨
efficienten naar model
In dit hoofdstuk gaan we de twee bewegingsvergelijkingen opstellen die horen bij de beweging van een
fiets. Om deze op te kunnen stellen beginnen we met het defini¨eren van alle parameters van de fiets.
Dit zijn enerzijds de fietsparameters, die onder een gekozen fiets vast liggen en anderszijds de parameters
die kunnen veranderen, de bewegingsparameters. Met al deze parameters werken we naar de kinetischeen potentiele energie toe, de basis van de Lagrangiaan en de externe krachten die daarbij komen kijken.
Vervolgens kunnen we vanuit de Lagrangiaan de bewegingsvergelijkingen afleiden.
3.1
Fysische betekenis van alle parameters
Voordat we daadwerkelijk beginnen met het benoemen van de parameters is het handig om een korte
situatieschets te maken. Er zijn namelijk een paar aannames die we nemen. De berijder van de fiets zit
vast aan de fiets en kan niet bewegen, dus vormt een starre verbinding met het frame. Maar de berijder
is niet star verbonden aan het stuur, de berijder fiets als het ware zonder handen. Het stuur is vrij om
te bewegen. De kuilen, hobbels of andere obstakels op de weg of de wind van op zij nemen we mee in de
kleine verstoringen van de beweging. Deze verstoringen nemen we niet apart mee in de externe krachten.
Voordat we punten, afstanden en hoeken bij de fiets kunnen defini¨eren moeten we de oorsprong en de
richtingen van de assen vaststellen. Met andere woorden, het assenstelsel moet worden vastgelegd. We
defini¨eren drie verschillende assenstelsels.
Het eerste assenstelsel wat we defini¨eren is het assenstelsel om de locatie van de fiets vast te kunnen
leggen. De oorsprong is vast en de X-as staat in de richting van de fiets in startpositie. De Y -as staat
daar loodrecht op in het grondvlak richting de linkerkant van de fiets. De Z-as staat daar vanzelfsprekend
weer loodrecht op en gaat de lucht in. Kort gezegd, een rechtshandig assenstelsel.
Definieer Sa als het contactpunt van het achterwiel met het grondvlak. x en y zijn de co¨ordinaten van
Sa ten opzichte van de oorsprong, met andere woorden de afstand x meet de afstand evenwijdig aan de
X-as tussen Sa en de Y -as en y is de afstand van Sa tot de X-as.
Omdat we ook de dynamica van de fiets bekijken gebruiken we nog twee extra assenstelsels, deze bewegen
met de fiets mee. Deze twee assenstelsels hebben elk hun oorsprong in een zwaartepunt. We defini¨eren
• Cv is het zwaartepunt van het voorwiel met het stuur en de voorvork, met andere woorden alles
wat met het stuur mee kan bewegen.
• Ca is het zwaartepunt van het achterwiel met de berijder en het achterwiel.
Wanneer de fiets verticaal staat met het voorwiel in het verlengde van het frame (fiets recht en stuur
recht), wijzen de X,Y en Z-as in dezelfde richting van het vaste assenstelsel zoals eerder is gedefinieerd:
X-as naar voren, Y -as naar links en Z-as naar boven.
3.1.1
De fietsparameters
Nu we de assenstelsels vastgelegd hebben kunnen we kijken naar de fietsparameters. Dit zijn de parameters die vastliggen door de keuze van de fiets. Hierbij kan gedacht worden aan de grootte van de wielen,
de hoogte van het stuur enzovoorts.
Deze fietsparameters zijn weergegeven in Figuur 1. De straal van het achterwiel noemen we ra en van
het voorwiel rv . De loodlijn door het middelpunt van het wiel noemen we respectievelijk Ba en Bv .
6
De stuuras is de as in het verlengde van de stang uit het frame waar je stuur aan vast zit. Dan is d de
lengte van de lijn loodrecht op de stuuras tot middelpunt Bv en e is de lengte van de lijn loodrecht op
de stuuras tot middelpunt Ba en tot slot is m de afstand tussen e en d.
Figuur 1: Vereenvoudigde weergave fiets naar Figuur 11 uit [1]
We hebben nu allemaal fietsparameters benoemd die direct zichtbaar zijn op de fiets. Welke fiets we ook
bekijken, we kunnen direct de grootte van de wielen aanwijzen en alle andere maten. Er zijn ook nog
andere fietsparameters die we niet direct met het blote oog kunnen zien. In Figuur 2 zijn die verschillende
fietsparameters te zien. We zullen ze stuk voor stuk aflopen, daarvoor hebben we nog het zwaartepunt
van de gehele fiets, Cf nodig. Het snijpunt van de stuuras met het grondvlak noemen we H.
Figuur 2: Aanvullig op Figuur 11 uit [1]
De volgende lengtes kunnen we nu benoemen. De afstand tussen Sa en Sv noemen we b, de afstand
tussen H en Sv noemen we de voorloop p. Dan defini¨eren we nog de hoogtes van de zwaartepunten, dus
de afstand tussen het zwaartepunt en het grondvlak. De hooogtes ha , hv en hf horen bij resp. Ca , Cv
en Cf .
Vanaf de loodlijn door Ba meten we de horizontale afstand tot Cf en Ca , deze noemen we uf en u. De
horizontale afstand tussen Cv en de loodlijn door Bv noemen we a.
7
De afstand f is de horizontale afstand tussen Cv en de stuuras.
Het snijpunt van de loodlijn door Cv met de stuuras is D. De afstand tussen Ca en D is w, de afstand
van D tot het grondvlak is k en de afstand van D tot de loodlijn door Ba is j.
Als laatste hebben we nog q dit is de afstand tussen H en Sa , enkel als de fiets geen stuurbeweging maakt,
dus rechtuit gaat, dan geldt q = b + p, zoals in Figuur 2.
Figuur 3: Vereenvoudigde weergave fiets naar Figuur 3 uit [1]
3.1.2
De bewegingsparameters
Om de beweging van de fiets te kunnen simuleren hebben we een aantal hoeken nodig, dit zijn de
bewegingsparameters, weergegeven in Figuur 3. We beginnen met de hoeken die aangeven hoe scheef een
wiel staat, dit zijn de hoeken δ en .
• De hoek δ geeft aan hoe schuin het voorwiel staat, met andere woorden δ wordt gemeten vanuit het
punt Sv en is de hoek tussen de verticaal door Sv en de loodlijn langs het voorwielvlak door Bv .
• De hoek geeft aan hoe schuin het achterwiel staat, anders gezegd: wordt gemeten vanuit het
punt Sa en is de hoek tussen de verticaal door Sa en de loodlijn langs het voorwielvlak door Ba .
De hoeken α, β en γ worden gemeten vanuit het punt H.
• α is de hoek tussen de stuuras en de snijlijn van het achterwielvlak met het grondvlak.
• β geeft aan hoever het stuur gedraaid is, met andere woorden β is de hoek tussen de snijlijnen van
het voor- en achterwielvlak met het grondvlak.
• γ is de hoek die aangeeft hoe schuin de fiets is, we meten γ als de hoek tussen de stuuras en de
snijlijn van het voorwielvlak met het grondvlak.
8
Hoek θ is de afwijking van het achterwiel ten opzichte van de X-as, met andere woorden θ is de hoek
tussen de snijlijn van het achterwielvlak met het grondvlak en de X-as. En als laatste, niet zichtbaar in
de figuren, is de hoek ψ die de verdraaiing van het stuur om zijn eigen as aangeeft.
Nu hebben we alle parameters benoemd, in Appendix A staan de complete afbeeldingen met parameters,
overgenomen uit artikel [1].
3.1.3
Welke parameters zijn constant en welke variabel
De twee ’interne’ zwaartepunten Ca en Cv zijn constant, daarmee ook hun afstanden tot het grondvlak,
ha en hv en ook de afstanden tot de bijbehorende wielvlakken: u en a. Hierdoor is f constant. Verder
zijn alle fietsparameters constant.
Variabel zijn alle hoeken, δ, , α, β, γ, θ en ψ en de afstanden p, b en q en de co¨ordinaten x en y van Sa ,
het raakpunt van het achterwiel met de grond.
Hoe groot de variaties zijn en tot welke orde we elke parameter meenemen gaan we nu bespreken.
3.2
Van co¨
efficienten naar energie
We hebben alle parameters weten van de fiets, inclusief de hoeken, afstanden en assenstelsels gedefinieerd.
Nu gaan we naar de energie¨en kijken die werken op de fiets. Hier gaat nog wat natuurkundig rekenwerk
aan vooraf. In deze paragraaf gaan we door dat rekenwerk heen.
Figuur 4: Boldriehoeksmeting, Figuur 4 uit [1]
Door wat algemeen bekende trigonometrische formules1 en Figuur 4, afkomstig uit artikel [1], te bestuderen komen we tot de volgende gelijkheden:
tan(δ)
=
sin(β)
sin(ψ)
=
sin(β)
+ cos(β) tan()
cos()
sin(α)
sin(γ)
=
cos(δ)
cos()
cot(α)
(1)
(2)
We bekijken nu steeds een hele kleine afwijking van de voorwaartse beweging van de fiets. Met andere
woorden y, θ, β, ψ, en δ zijn oneindig klein en we nemen ze tot de 2e orde mee en α, γ, q, p en b zijn ook
oneindig klein en die nemen we mee tot orde 1.
1 Bron:
http://nl.wikipedia.org/wiki/Boldriehoeksmeting
9
We noemen α0 de waarde die α en γ hebben als β = ψ = 0. En we noemen ∆α en ∆γ de kleine
veranderingen van hoeken α en γ. Dan zijn α en γ:
α
= α0 + ∆α
γ
= α0 + ∆γ
(3)
Dit kunnen we substitueren in (2), dan
sin(α0 + ∆α) cos() = sin(α0 + ∆γ) cos(δ)
Beide kanten kunnen we de somregel van sinus toepassen:
sin(α0 ) cos(∆α) + cos(α0 ) sin(∆α) cos() = sin(α0 ) cos(∆γ) + cos(α0 ) sin(∆γ) cos(δ)
We kunnen nu de reeksontwikkeling invullen, maar de termer hoger dan de 2e orde kunnen we verwaarlozen, dan is
2
∆α2
+ . . . + cos(α0 )(∆α − . . .)
1−
+ ...
sin(α0 ) 1 −
2
2
∆γ 2
δ2
= sin(α0 ) 1 −
+ . . . + cos(α0 )(∆γ − . . .)
+ ...
1−
2
2
Dit uitwerken geeft
sin(α0 ) − sin(α0 )
2
2
∆α2
2
∆α2
+ cos(α0 )∆α − sin(α0 ) + sin(α0 )
− cos(α0 )∆α
2
2
2
2
2
= sin(α0 ) − sin(α0 )
∆γ 2
δ2
δ2
∆γ 2
δ2
+ cos(α0 )∆γ −
sin(α0 ) +
sin(α0 )
−
cos(α0 )∆γ
2
2
2
2
2
Nu delen we alles door cos(α0 ), dan
− tan(α0 )
∆α2
2
2 ∆α2
2
+ ∆α − tan(α0 ) +
tan(α0 ) − ∆α
2
2
2 2
2
= − tan(α0 )
∆γ 2
δ2
δ 2 ∆γ 2
δ2
+ ∆γ −
tan(α0 ) +
tan(α0 ) − ∆γ
2
2
2 2
2
We zouden de termen van hogere orde dan twee verwaarlozen en α en γ nemen we maar tot de eerste
orde mee dus
∆α −
2
δ2
tan(α0 ) = ∆γ −
tan(α0 )
2
2
Dit kunnen we herschrijven naar
∆α − ∆γ =
1
tan(α0 )(2 − δ 2 )
2
(4)
Noem nu H0 het punt op de stuuras dat op het grondvlak ligt als β = 0 en het stuur vooruit staat. Dan
zal H0 een klein stukje boven of onder het grondvlak komen te liggen als het stuur een kleine afwijking
maakt.
In [1] wordt nu de volgende observatie gedaan
p∆γ = q∆α
10
(5)
We combineren nu (4) en (5) en dan vinden we
p 1
tan(α0 )(2 − δ 2 )
q−p2
q 1
tan(α0 )(2 − δ 2 )
−
q−p2
∆α
= −
∆γ
=
(6)
Wanneer we nu (1) gelineairiseerd bekijken vinden we ook:
δ = β cos(α) + (7)
Nu substitueren we (7) in (6), bovendien weten we dat q − p = b. We vullen we dit in, dan krijgen we:
∆α
=
∆γ
=
p
1
β + β 2 cot(α)
b
2
1
p+b
β + β 2 cot(α)
b
2
(8)
Merk op: door het lineariseren hebben we wel een kleine fout toegevoegd aan de uitkomst, maar deze
blijft heel klein, hoogstens orde 2.
Hiermee veranderen een aantal parameters die variabel waren in constanten, die we in H3.1.3 hadden
gevonden. Vanaf nu zijn x, y, θ, en β variabel.
De parameters ra , rv , p, b, cot(α), Ixa , Iza , Ixza , Ixv , Izv , Ixzv , Mv en Ma zijn vanaf nu constant.
3.3
Van energie naar Lagrange vergelijkingen
We zijn uiteindelijk op zoek naar de bewegingsvergelijkingen. We gaan deze vinden met behulp van de
Methode van Lagrange. Op de grond, bij de wielen, werken er nog twee reactiekrachten Pa en Pv op de
fiets die onbekend zijn. Door de Methode van Lagrange kunnen we deze krachten elimineren. Verder
maken we eerst nog twee opmerkingen zodat het rekenwerk met Lagrange een stuk eenvoudiger wordt.
1. De co¨
ordinaat x kan weggelaten worden. Want in de X-richting werken alleen krachten die heel erg
klein zijn, namelijk hoogstens van orde 2. Daarmee kunnen we ook de term 12 M V 2 weglaten in de
kinetische energie, want V , de snelheid in de X-richting, is een constante term.
2. Het gyroscopisch effect van de fiets kan ondergebracht worden in de energie, maar ook onder de
momenten. We gaan dit laatste doen.
Voordat we daadwerkelijk aan de Langrange-vergelijkingen gaan beginnen, gaan we eerst de onderdelen
apart berekenen. We berekenen eerst de kinetische energie, vervolgens de potenti¨ele energie en dan nog
de momenten. Dan hebben we genoeg informatie om de Langrangiaan en externe krachten te berekenen.
De kinetische energie, T , kunnen we onderverdelen in vier sub-energi¨en:
T = Tza + Toa + Tzv + Tov
11
(9)
Waarbij de sub-energi¨en zijn:
Tza
=
de energie van het achterste gedeelte van de fiets door translatie van het zwaartepunt
2
= 12 Ma y˙ + ha ˙ + uθ˙
Toa
=
Tzv
=
Tov
=
de energie van het achterste gedeelte van de fiets door rotatie van het zwaartepunt
= 12 Ixa ˙2 + 12 Iza θ˙2 + Ixza ˙θ˙
de energie van het voorste gedeelte van de fiets door translatie van het zwaartepunt
2
= 12 Mv y˙ + f β˙ + (b − a)θ˙ + hv ˙
(10)
de energie van het voorste gedeelte van de fiets door rotatie van het zwaartepunt
˙ 2 + Ixzv (β˙ cot(α) + )(
˙
˙ 2 + 12 Izv (β˙ + θ)
˙ β˙ + θ)
= 12 Ixv (β˙ cot(α) + )
Zoals gegeven in [1]. De totale kinetische energie is dan:
T
=
1
2 Ma
y˙ + ha ˙ + uθ˙
2
2
+ 12 Ixa ˙2 + 12 Iza θ˙2 + Ixza ˙θ˙ + 12 Mv y˙ + f β˙ + (b − a)θ˙ + hv ˙
˙ 2 + Ixzv (β˙ cot(α) + )(
˙
+ 12 Ixv (β˙ cot(α) + )
˙ 2 + 21 Izv (β˙ + θ)
˙ β˙ + θ)
(11)
De potentiele energie, Epot , is zoals in [1]:
Epot
1
1 2
2
= −gMa
β + β cot(α) + ha 2
2
1 2
1
a(b + p)
2
β + β cot(α) + hv ( + β cot(α))
−gMv −
b
2
2
pu
b
Dan zijn er nog de gyroscopische momenten. Deze kunnen we aflezen in Figuur 5
Figuur 5: Zoals Figuur 9 en 10 uit [1]
12
(12)
Dan krijgen we de volgende vergelijkingen:
Mxv
=
=
Mzv
=
=
Mxa
Mza
3.4
=
het moment van het voorste gedeelte van de fiets ten opzichte van de x-as
˙
Iwv ωv (θ˙ + β)
het moment van het voorste gedeelte van de fiets ten opzichte van de z-as
Iwv ωv (β˙ cot(α) + )
˙
=
het moment van het achterste gedeelte van de fiets ten opzichte van de x-as
Iwa ωa θ˙
=
het moment van het achterste gedeelte van de fiets ten opzichte van de z-as
=
Iwa ωa ˙
(13)
Euler-Lagrange
We hebben nu genoeg informatie om de Lagrangiaan op te stellen, want de Lagrangiaan is L = T − Epot .
De Lagrangiaanse vergelijkingen, ten op zichte van de co¨ordinaten qi , zijn:
∂L
d ∂L
−
= Qi
dt ∂ q˙i
∂qi
(14)
Waarbij Qi de externe krachten zijn die niet bevat zijn in de energie¨en. Dit kunnen we herschrijven als:
d
dt
∂T
∂Epot
−
∂ q˙i
∂ q˙i
−
∂T
∂Epot
+
= Qi
∂qi
∂qi
(15)
We scheiden nu de kinetische energie van de potentiele energie, dan
∂T
d ∂Epot
∂Epot
d ∂T
−
=
−
+ Qi
dt ∂ q˙i
∂qi
dt ∂ q˙i
∂qi
(16)
Wanneer we kijken naar de kinetische energie¨en, zoals in (10), zien we dat deze enkel van de veranderingen
∂T
van de co¨
ordinaten afhangt en niet van de co¨ordinaten zelf. Dus
= 0 voor alle qi . Als we naar de
∂qi
potentiele energie kijken, zoals in (12), zien we dat deze enkel van de co¨ordinaten afhangen en niet van
∂Epot
de verandering van de co¨
ordinaten. Dus
= 0 voor alle q˙i . Dan
∂ q˙i
d ∂T
∂Epot
=−
+ Qi
dt ∂ q˙i
∂qi
(17)
Er zijn 4 co¨
ordinaten, namelijk: y, θ, en β. Herinner uit H.3.1.2:
• y is de co¨
ordinaat van de locatie van het contactpunt van het achterste wiel met de grond.
• θ is de afwijking van het achterwiel ten opzichte van de X-as.
• is hoe schuin het achterwiel staat.
• β is hoe ver het stuur gedraaid is.
We hebben dus 4 Lagrangiaanse vergelijkingen. Nu werken er nog twee externe krachten Pa en Pv ,
die voorkomen dat de wielen gaan slippen. En bovendien werken de momenten nog. Deze moeten we
meenemen in de Lagrange vergelijkingen. Dan:
13
d ∂T
dt ∂ y˙
= −
∂Epot
+ Pa + Pv = αy
∂y
(18)
d ∂T
dt ∂ θ˙
= −
∂Epot
+ bPv + Mza + Mzv = αθ
∂θ
(19)
d ∂T
dt ∂ ˙
= −
∂Epot
− Mxv − Mxa = α
∂
(20)
d ∂T
dt ∂ β˙
= −
∂Epot
∂δ
− pPv + Mzv − Mxv
= αβ
∂β
∂β
(21)
Al deze termen kunnen we stuk voor stuk uitrekenen, behalve Pa en Pv . De rechterkant van elke vergelijking geven we een naam: αy , αθ , α , αβ . We refereren aan deze namen bij het uitrekenen van de termen.
Voordat we daar aan toe komen hebben we nog een klein beetje werk te doen, namelijk we moeten de
voorwaarden voor het voorkomen van het slippen van de wielen nog mee nemen. Deze voorwaarden
komen uit [1].
V θ − y˙
V (θ + β) − y˙ − bθ˙ + pβ˙
=
0
=
0
Dit kunnen we schrijven als:
y˙
= Vθ
y˙
= V θ + V β − bθ˙ + pβ˙
Samenvoegen geeft:
V β − bθ˙ + pβ˙ = 0
Herschrijven geeft:
V β + pβ˙
θ˙ =
b
(22)
θ¨ =
V β˙ + pβ¨
b
(23)
y¨ =
V
˙
(V β + pβ)
b
(24)
Dan
˙ θ¨ en
Door deze voorwaarden te gebruiken in de Lagrangiaanse vergelijkingen kunnen we de termen van θ,
y¨ elimineren. Met behulp van artikel [1] hebben we nu al het gereedschap verzameld om de Lagrangiaanse
vergelijkingen te berekenen, de bewegingsvergelijkingen op te beschrijven en daarmee zelf het onderzoek
naar de stabiliteit van de fiets te vervolgen.
We willen nu alle termen van de Lagrangiaanse vergelijking uitrekenen om daaruit de bewegingsvergelijkingen te halen. Omdat de berekeningen voor elke Lagrangiaanse vergelijking hetzelfde zijn, zullen
we het enkel voor de Lagrangiaanse vergelijking van doen, de rest gaat analoog en zullen we enkel
als resultaat weergeven. De verificatie van de berekeningen is gedaan in Maple, zie Appendix C.3, ook
hier zijn enkel de berekeningen voor uitgevoerd, voor de andere Lagrangiaanse vergelijkingen gaat dit
natuurlijk analoog.
14
3.4.1
Uitrekenen van termen met Zoals gezegd gaan we nu beginnen met het oplossen van de Lagrangiaanse vergelijkingen, we beginnen
met het uitrekenen met de linkerkant van (20).
d ∂T
dt ∂ ˙
= Ma y¨ + ha ¨ + uθ¨ ha + Mv y¨ + f β¨ + (b − a)θ¨ + hv ¨ hv + Ixa ¨ + Ixza θ¨
¨
+Ixv (β¨ cot(α) + ¨) + Ixzv (β¨ + θ)
(25)
Dan substitueren we de slipvoorwaarden, (23) en (24) in vergelijking (25), dan
d ∂T
dt ∂ ˙
=
Ma
+Mv
˙
¨
V
˙ + ha ¨ + u V β + pβ
(V β + pβ)
b
b
!
ha
!
V
V β˙ + pβ¨
V β˙ + pβ¨
˙
¨
(V β + pβ) + f β + (b − a)
+ hv ¨ hv + Ixa ¨ + Ixza
b
b
b
+Ixv (β¨ cot(α) + ¨) + Ixzv (β¨ +
V β˙ + pβ¨
)
b
(26)
¨ β˙ en β van elkaar gescheiden zijn.
Dit kunnen we herschrijven zo dat de termen ¨, ,
˙ , β,
d ∂T
dt ∂ ˙
=
Ma h2a + Mv h2v + Ixa + Ixv ¨
p
p
p
b+p ¨
+ Ma ha u + Mv hv f + Mv hv (b − a) + Ixza + Ixv cot(α) + Ixzv
β
b
b
b
b
+ {Ma ha p + Ma ha u + Mv hv p + Mv hv (b − a) + Ixza + Ixzv }
+ {Ma ha + Mv hv }
V2
β
b
V ˙
β
b
(27)
We hebben de linkerkant van de Lagrangiaanse vergelijking (20) berekend, nu is de rechterkant, α aan
de beurt, merk op dat we de potenti¨ele energie al weten uit (12) en ook de momenten zijn al bekend uit
(13).
15
α
∂Epot
− Mxv − Mxa
∂
=
−
=
gMa
=
pu
a(b + p)
˙
{gMa ha + gMv hv } − Iwv ωv β + gMa
+ gMv hv cot(α) − gMv
β
b
b
a(b + p)
β + ha + gMv hv ( + β cot(α)) −
β − Iwv ωv θ˙ + β˙ − Iwa ωa θ˙
b
b
pu
− {Iwv ωv + Iwa ωa } θ˙
(28)
Wanneer we nu het verschil van de linker en de rechterkant bekijken moet daar 0 uit komen, dan scheiden
¨ β˙ en β van elkaar. Dan volgt dat:
we weer direct ¨, ,
˙ , β,
d ∂T
− α
dt ∂ ˙
=
0
(29)
=
0
(30)
Ma h2a + Mv h2v + Ixa + Ixv ¨ − {gMa ha + gMv hv } p
p
b+p ¨
p
β
+ Ma ha u + Mv hv f + Mv hv (b − a) + Ixza + Ixv cot(α) + Ixzv
b
b
b
b
V
+ (Ma ha (p + u) + Mv hv p + Mv hv (b − a) + Ixza + Ixzv ) + Iwv ωv β˙
b
pu
a(b + p)
V2
− gMa
+ gMv hv cot(α) − gMv
β
+ (Ma ha + Mv hv )
b
b
b
+ {Iwv ωv + Iwa ωa } θ˙
Hierin substitueren we nog de slipvoorwaarde (22) in, dan:
d ∂T
− α
dt ∂ ˙
=
0(31)
=
0(32)
Ma h2a + Mv h2v + Ixa + Ixv ¨ − {gMa ha + gMv hv } p
p
p
b+p ¨
+ Ma ha u + Mv hv f + Mv hv (b − a) + Ixza + Ixv cot(α) + Ixzv
β
b
b
b
b
V
b+p
p ˙
+ (Ma ha (p + u) + Mv hv p + Mv hv (b − a) + Ixza + Ixzv ) + Iwv ωv
+ Iwa ωa
β
b
b
b
+ (Ma ha + Mv hv )
V2
V
+ (Iwv ωv + Iwa ωa )
b
b
pu
a(b + p)
β
− gMa
+ gMv hv cot(α) − gMv
b
b
16
3.5
Overige Lagrangiaanse vergelijkingen
We hebben dezelfde berekeningen voor de andere Lagrangiaanse vergelijkingen uitgevoerd en deze ook
weer geverifieerd met Maple. Dit zijn de resultaten:
d ∂T
− αy
dt ∂ y˙
=
0
V ˙
V2
β + {Ma + Mv }
β − P a − Pv
b
b
=
0 (33)
d ∂T
− αθ
dt ∂ θ˙
=
0
V2
β − bP v
b
=
0 (34)
d ∂T
− αβ
dt ∂ β˙
=
0
=
0 (35)
n
p
po ¨
β
{Ma ha + Mv hv } ¨ + Ma u + Mv f + Mv (b − a)
b
b
+ {Ma p + Ma u + Mv p + Mv (b − a)}
{uMa ha + Mv hv (b − a) + Ixza + Ixzv } ¨ − (Iwv ωv + Iwa ωa )˙
n
o
p
p
p
p
+ u2 Ma + (b − a)Mv f + Mv (b − a)2 + Iza + Izv + Izv + Ixzv cot(α) β¨
b
b
b
b
+
V
uMa p + u Ma + (b − a)Mv p + Mv (b − a) + Iza + Izv
− Iwv ωv cot(α) β˙
b
2
2
+ {uMa + (b − a)Mv }
a(b + p)
pu
+ gMv hv cot(α) − gMv
{Mv f hv + Ixzv + Ixv cot(α)} ¨ − Iwv ωv ˙ − gMa
b
b
p
b+p
2b + p
+ Mv f 2 + Mv f (b − a) + Izv
+ Ixzv
cot(α) + Ixv cot(α)2 β¨
b
b
b
p
V
+ (Mv f p + Mv f (b − a) + Izv + Ixzv cot(α)) + Iwv ωv cot(α)
b
b
β˙
V
pu
a(b + p)
V2
+ Mv f
+ Iwv ωv cot(α) − gMa
+ gMv hv cot(α) − gMv
cot(α) β
b
b
b
b
+pPv
17
3.6
Bewegingsvergelijkingen opstellen
We nemen de volgende bewegingsvergelijkingen. De eerste bewegingsvergelijking is de derde Lagrange
vergelijking:
d ∂T
− α = 0
dt ∂ ˙
Als tweede bewegingsvergelijking is een lineaire combinatie van de tweede en vierde Lagrange vergelijking,
namelijk:
d ∂T
p
− αβ +
dt ∂ β˙
b
d ∂T
− αθ
dt ∂ θ˙
=0
We zullen, in de tweede bewegingsvergelijking, zien dat door een lineaire combinatie van die twee Lagrangevergelijkingen te nemen, we de term Pv elimineren. Dit motiveert de keuze van de tweede bewegingsvergelijking.
We zullen nu de lineaire combinatie beter bekijken. We willen de term Pv elimineren, deze term komt
alleen voor in αβ en αθ , dus bekijken we de lineaire combinatie: αβ + pb αθ . Herinner:
αβ
αθ
a(b + p)
pu
+ gMv hv cot(α) − gMv
= Iwv ωv ˙ + gMa
b
b
=
a(b + p)
pu
+ gMv hv cot(α) − gMv
cot(α)β − Iwv ωv cot(α)θ˙ − pPv
+ gMa
b
b
(36)
(Iwv ωv + Iwa ωa )˙ + Iwv ωv cot(α)β˙ + bPv
(37)
We vullen dit in de lineaire combinatie in en dit geeft:
p
pu
a(b + p)
αβ + αθ = Iwv ωv ˙ + gMa
+ gMv hv cot(α) − gMv
b
b
b
pu
a(b + p)
+ gMa
+ gMv hv cot(α) − gMv
b
b
+
cot(α)β − Iwv ωv cot(α)θ˙ − pPv
i
ph
(Iwv ωv + Iwa ωa )˙ + Iwv ωv cot(α)β˙ + bPv
b
(38)
Dan zien we nu al dat er een term bPv en −bPv is. We kunnen Pv dus inderdaad elimineren:
p
αβ + αθ =
b
pu
a(b + p)
gMa
+ gMv hv cot(α) − gMv
+ cot(α)β
b
b
+Iwv ωv
b+p
p
p
˙ + cot(α)β˙ − cot(α)θ˙ + Iwa ωa ˙
b
b
b
(39)
Door de lineaire combinatie te nemen van de Lagrange vergelijkingen, kunnen we nu iets nuttigs zeggen
over de bewegingsvergelijkingen, want daar zitten nu alleen maar termen in waar we de afkomst van
kennen.
18
3.7
Bewegingsvergelijkingen vereenvoudigen
We hebben nu de bewegingsvergelijkingen gevonden. We schrijven ze voor de volledigheid even uit. De
eerste bewegingsvergelijking is:
(Ixa + Ixv + h2v Mv + h2a Ma )¨
− g(ha Ma + hv Mv )
+
np
b
o
Ixza + Ixzv + uha Ma + (b − a)hv Mv + Ixzv + Ixv cot(α) + hv f Mv β¨
b+p
p
V ˙
Iwv ωv + Iwa ωa
β
b
b
b
V2
p
V
V
+ (hv Mv + ha Ma )
− (uMa + (b − a)Mv )g − gf Mv + Iwv ωv + Iwa ωa
β=0
b
b
b
b
+ Ixza + Ixzv + uha Ma + (b − a)hv Mv + p(ha Ma + hv Mv ) +
(40)
De tweede bewegingsvergelijking is:
n
o
p
Ixzv + Ixv cot(α) + hv f Mv + Ixza + Ixzv + uha Ma + (b − a)hv Mv ¨
b
n
o
np
o
p
− Iwv ωv + (Iwv ωv + Iwa ωa ) ˙ −
uMa + (b − a)Mv g + f gMv b
b
p2
+ Ixv cot2 (α) + 2Ixzv cot(α) + Izv + f 2 Mv + 2 Iza + Izv + u2 Ma + (b − a)2 Mv
b
2p
Izv + Ixzv cot(α) + (b − a)f Mv β¨
+
b
n
p
+ Izv + Ixzv cot(α) + (b − a)f Mv + pf Mv + Iza + Izv + u2 Ma + (b − a)2 Mv
b
V ˙
p2
β
+ (uMa + (b − a)Mv )
b
b
V 2
pu
p
v
+
Ma + (b − a)Mv ·
− g cot(α) + Iwv ωv cot(α) β = 0
b
b
b
b
(41)
Dit zijn nog vrij ingewikkelde vergelijkingen, dus we voeren wat notatie in om de vergelijkingen te
vereenvoudigen:
Ixz
= Ixza + Ixzv + uha Ma + (b − a)hv Mv
Ix
= Ixa + Ixv + h2v Mv + h2a Ma
Iz
= Iza + Izv + u2 Ma + (b − a)2 Mv
Sx
= Mv hv + Ma ha
Sz
= Mv (b − a) + Ma u
∆
=
Γ
= Izv + Ixzv cot(α) + (b − a)f Mv = Ias − Ian cot(α) + jf Mv
k
= hv − f cos(α) sin(α)
j
= b − a − f sin2 (α)
ωv
=
ωa
=
Ixzv + Ixv cot(α) + hv f Mv = Ias cot(α) + Ian + kf Mv
b
rv
b
ra
19
Dan zien de bewegingsvergelijkingen er zo uit:
Ix ¨ − gSx +
+
V
b+p
b
p
b ˙
Ixz + ∆ β¨ +
Ixz + pSx +
Iwv + Iwa
β
b
b
b
rv
b
ra
hp
i
Iwa V 2
+
− g Sz + f Mv β = 0
ra
b
b
p
Iwv
Sx +
rv
i
hp
i
b + p Iwv
p Iwa
Ixz + ∆ ¨ − V
+
˙ − g Sz + f Mv b
b rv
b ra
b
h
i
P
Ias
p2
V p
p
¨
I
+
2
+
+
Γ
β
+
I
+
p
S
+
f
M
+
Γ
β˙
z
z
z
v
b
b b
b
sin2 (α) b2
2
hp
i
p
Iwv
V
+
Sz + f Mv +
cot(α)
− g Sz + f Mv cot(α) β = 0
b
rv
b
b
(42)
hp
(43)
In beide vergelijkingen komen β en voor met hun eerste en tweede afgeleide. Dit kunnen we schrijven
als een stelsel:
a1 ¨
+a3 b1 ¨ +b2 ˙ +b3 +a4 β¨ +a5 β˙
+b4 β¨ +b5 β˙
+a6 β
=0
+b6 β
=0
(44)
waarbij de co¨efficienten dan gelijk zijn aan:2
2 dit
a1
=
Ix ;
a2
=
0;
a3
=
−gSx ;
a4
=
a5
=
a6
=
b1
=
b2
=
b3
=
b4
=
b5
=
b6
=
p
Ixz + ∆;
b
V
b+p
b
p
b
Ixz + pSx +
Iwv + Iwa
b
b
rv
b
ra
2
h
i
Iwv
Iwa V
p
Sx +
+
− g Sz + f Mv ;
rv
ra
b
b
p
Ixz + ∆;
b
b + p Iwv
p Iwa
−V
+
;
b rv
b ra
hp
i
−g Sz + f Mv ;
b
Ias
p2
P
+ 2 Iz + 2 Γ;
2
b
sin (α) b
h
i
V p
p
Iz + p
Sz + f Mv + Γ ;
b b
b
2
hp
i
p
Iwv
V
Sz + f Mv +
cot(α)
− g Sz + f Mv cot(α);
b
rv
b
b
zijn dezelfde co¨
efficienten als die in [1] gebruikt zijn
20
(45)
4
Structuur onderzoek stabiliteit
We hebben nu twee differentiaalvergelijkingen, zie (44), die we kunnen schrijven als een stelsel:
"
a1
a4
b1
b4
#"
¨
¨
β
"
#
+
0
a5
b2
b5
#"
˙
β˙
#
"
+
a3
a6
b3
b6
#"
#
β
=0
Dit kunnen we nog iets beter schrijven zodat we de structuur van het stelsel beter kunnen zien. Bekijk
(45), dan zien we dat:
1. in a5 , b2 en b5 overal termen staan met de snelheid.
2. in a3 , a6 , b3 en b6 zien we bij elke term een v 2 en anders een g.
3. in a3 en b3 staan geen termen met v 2 , dus daar komt een 0.
We kunnen het stelsel als volgt schrijven: M¨
q + vCq˙ + (gK + v 2 P)q = 0. Hierbij zijn M, C, K en P 2 × 2
matrices, q een vector, g de gravitatie en v de snelheid van de fiets. De algemene vorm van het probleem
ziet er zo uit:
"
m11
m12
m12
m22
met
#"
"
M=
¨
β¨
#
"
+v
m11
m12
m12
m22
0
c12
c21
c22
#
#"
"
, C=
˙
β˙
#
( "
+
0
c12
c21
c22
g
k11
k12
k12
k22
#
"
, K=
"
#
+v
2
k11
k12
k12
k22
0
p12
0
p22
#
#) "
#
=0
β
"
, P =
0
p12
0
p22
(46)
#
De oplossingen van dit stelsel vergelijkingen zijn van de vorm q = q0 eλt . Deze oplossing vullen we in,
dan krijgen we
Mλ2 + vCλ + gK + v 2 P q0 eλt = 0
We willen uiteindelijk iets over de stabiliteit van de fiets bij een bepaalde snelheid kunnen zeggen. Daarom
kijken we naar de karakteristieke vergelijking. Voor de stabiliteit kijken we eerst naar de oplossingen
van die karakteristieke vergelijking en vervolgens naar het Routh Stabilty Criteria, die kijkt naar de
co¨efficenten van de karakteristieke vergelijking.
In de rest van dit hoofdstuk zullen we eerst de wiskundige aanpak uitleggen. In de volgende hoofdstukken
kijken we naar de toepassing op ons model, de restultaten die we vinden en de conclusies die we daaruit
kunnen trekken.
4.1
Karakteristieke vergelijking
Wanneer we naar de eigenwaarden kijken geldt dat alle eigenwaarden een negatief re¨eel deel moeten
hebben. Dus we moeten de eigenwaarden berekenen. Dit doen we door te kijken naar de determinant
van de matrix A := Mλ2 + vCλ + gK + v 2 P. Wanneer de determinant niet gelijk is aan 0, dan is de
matrix inverteerbaar, maar dan vinden we enkel de triviale oplossing namelijk q0 eλt = 0, want:
Aq0 eλt = 0 ⇒ A−1 Aq0 eλt = A−1 0 ⇒ q0 eλt = 0
Deze vergelijking heeft dus alleen niet-triviale oplossingen als det(Mλ2 + Cλ + K) = 0. De matrix A ziet
er als volgt uit
"
A=
m11 λ2 + gk11
c12 vλ + m12 λ2 + p12 v 2 + gk12
c21 vλ + m12 λ2 + gk12
c22 vλ + m22 λ2 + p22 v 2 + gk22
21
#
De karakteristieke vergelijking, det(A) = 0, is dan
(m11 m22 − m212 )λ4
+v(m11 c22 − (c12 + c21 )m12 )λ3
2
+[(−c12 c21 + m11 p22 − m12 p12 )v + g(k11 m22 − 2k12 m12 + k22 m11 )]λ2
−v[p12 c21 v 2 + ((c12 + c21 )k12 − c22 k11 )g]λ
2
+g[(k11 p22 − k12 p12 )v 2 + g(k11 k22 − k12
)]
4.1.1
=0
Oplossingen karakteristieke vergelijking
De karakteristieke vergelijking is een 4e graads polynoom, er zijn dus 4 eigenwaarden, λ1 , λ2 , λ3 en λ4 .
Voor stabiliteit moeten die allemaal negatief zijn. Merk op dat de eigenwaarden afhankelijk zijn van de
snelheid. De stabiliteit is dus afhankelijk van de snelheid.
4.1.2
Routh Stability Criteria
We bekijken het Routh Stability Criteria zoals uit gelegd wordt in [2] in Chapter 4. Het Routh Stability
Criteria zegt iets over de nulpunten van een 4e graads polynoom. Als we kijken naar een standaard 4e
graads polynoom ziet die er zo uit:
Aλ4 + Bλ3 + Cλ2 + Dλ + E = 0
Het Routh Stability Criteria definieert de Routh-determinant als volgt: X = BCD − ADD − EBB.
Het Routh Stability Criteria zegt dat als voor een gegeven snelheid A, B, C, D, E en X hetzelfde teken
hebben, dus allemaal positief of negatief zijn, dan is de fiets met die snelheid stabiel. Het karakteristieke
vergelijking zoals eerder gevonden is een 4e graads polynoom. Dus kunnen we hier direct het Routh
Stability Criteria op toepassen:
A
= m11 m22 − m212
B
= v(m11 c22 − (c12 + c21 )m12 )
C
=
(−c12 c21 + m11 p22 − m12 p12 )v 2 + g(k11 m22 − 2k12 m12 + k22 m11 )
(47)
3
D
= p12 c21 v + g((c12 + c21 )k12 − c22 k11 )v
E
2
= g(k11 p22 − k12 p12 )v 2 + g 2 (k11 k22 − k12
)
De Routh determinant X staat in Appendix E. Om de Routh determinant iets overzichtelijker te maken
vereenvoudigen we de termen A, B, C, D en E.
A
= A1
B
= B1 v
C
= C1 v 2 + C2
D
= D1 v 3 + D2 v
E
= E1 v 2 + E2
De Routh determinant X is dan
X = B1 C1 D1 − A1 D12 v 6 + B1 C1 D2 + B1 C2 D1 − 2A1 D1 D2 − B12 E1 v 4 + B1 C2 D2 − A1 D22 − B12 E2 v 2
22
Wanneer A, B, C, D, E en X hetzelfde teken hebben is de fiets stabiel voor die snelheid. We zien nu al
dat A onafhankelijk is van de snelheid. Dus A bepaalt welk teken B, C, D, E en X moeten hebben voor
stabiliteit.
5
Waardes co¨
efficienten invullen
Omdat we nu de daadwerkelijke oplossingen van het stelsel gaan berekenen, vullen we eerst de waarden
in, zodat de matrices gevuld worden met getallen. We vullen ze in naar aanleiding van (45) en de waarden
uit B zoals in [1]. We krijgen dan de volgende stelsels:
"
98
1.77
1.77
0.36
#"
# "
#"
# ( "
"
#
¨
0 26.5
˙
−82.5 −2.08
0
2
+v
+ g
+v
¨
˙
β
−0.79 1.37
β
−2.08 −0.97
0
71.71
#) "
2.06
#
=0
β
(48)
Wanneer we dit vergelijken met [3], zien we een soortgelijk stelsel. Die ziet er als volgt uit:
"
81
2.3
2.3
0.3
#"
φ¨
δ¨
#
"
+v
0
34
−0.9
1.7
#"
φ˙
δ˙
#
( "
+
g
−81 −2.6
−2.6
−0.8
#
"
+v
2
0
77
0
2.7
#) "
φ
δ
#
=0
De vector q is een vector bestaande uit 2 elementen, de eerste is φ, de leunhoek van de fiets, hoe schuin
de fiets staat, en de tweede is δ, de stuurhoek van de fiets, hoe ver het stuur draait. Dit komt overeen
met wat wij en β hebben genoemd.
We kunnen dus de stelsels met elkaar vergelijken. En we zien dat alle co¨efficienten wel heel veel op elkaar
lijken. Er zitten maar hele kleine afwijkingen tussen elk paar co¨efficienten. We kunnen dus wel aannemen
dat de resultaten ook vergelijkbaar zullen zijn en ook de beide fietsen, die als uitgangspunt genomen zijn,
erg op elkaar zullen lijken.
Maar omdat we van de co¨efficienten uit [3] niet precies weten waar ze vandaan komen gaan we verder
met de co¨efficienten uit [1], want daarvan kennen we precies de herkomst. We gaan dus met (48) verder.
6
Resultaten
Zoals eerder besproken willen we nu de karakteristieke vergelijking bepalen en daarmee de eigenwaarden
berekenen en het Routh Stability Criteria op toepassen. We hebben al gezien dat de oplossing van de
vorm q = q0 eλt is. Dit vullen we in, dan krijgen we
"
98
1.77
1.77
0.36
#
"
2
λ +v
0
26.5
−0.79
1.37
#
"
λ+g
−82.5
−2.08
−2.08 −0.97
#
"
+v
2
0
71.71
0
2.06
Dit kunnen we samenvoegen tot de matrix A. Dit vullen we in, dan
"
A=
98λ2 − 809.33
1.77λ2 + 26.5vλ − 20.4 + 71.71v 2
1.77λ2 − 0.79vλ − 20.4
0.36λ2 + 1.37vλ − 9.52 + 2.06v 2
23
#
#!
q0 eλt = 0
De karakteristieke vergelijking is dan
det(A) = 32.15λ4 + 88.75vλ3 + (95.89v 2 − 1151.66)λ2 + (56.65v 3 − 584.17v)λ + 7284.94 − 203.98v 2
Deze vergelijking hebben we nodig bij het berekenen van de eigenwaarden en bij het toepassen van het
Routh Stability Criteria. We merken op dat deze vergelijking enkel van de snelheid afhangt, dus ook de
stabiliteit zal afhangen van de snelheid van de fiets.
6.1
Eigenwaarden & Stabiliteit
We hebben nu het karakteristieke polynoom. Door het karakteristiek polynoom gelijk te stellen aan 0,
krijgen we de karakteristieke vergelijking. Hieruit kunnen we de eigenwaarden, λ1 , λ2 , λ3 en λ4 bepalen
door de vergelijking op te lossen. De eigenwaarden die we vinden hangen af van de snelheid. Daarom
berekenen we de eigenwaarden voor een gegeven snelheid, vervolgens plotten we de eigenwaarden tegenover
de snelheid. Wanneer alle re¨ele delen van de eigenwaarden negatief zijn is de beweging van de fiets stabiel.
Wanneer v = 0 verwachten we dat de fiets niet stabiel is, want in de praktijk valt een fiets zonder snelheid
gewoon om. Als de snelheid v = 0, dan is de karakteristieke vergelijking:
32.15λ4 − 1151.66λ2 + 7284.94 = 0
We kunnen nu al snel zien dat niet alle eigenwaarden een negatief re¨eel deel hebben. Wanneer we w = λ2
in de vergelijking substitueren kunnen we de oplossingen berekenen. Dan vinden we 2 oplossingen voor
√
w, beide zijn voor v = 0 positief. Voor λ vinden we dan λ = ± w. Met andere woorden, er is een
positieve en een negatieve eigenwaarde. De eigenwaarden als v = 0 zijn:
λ1 = 2.86, λ2 = −2.86, λ3 = 5.26 en λ4 = −5.26
Als v = 1 is de karakteristieke vergelijking
det(A) = 32.15λ4 + 88.75λ3 − 1055.77λ2 − 527.52λ + 7080.96
De eigenwaarden horende bij v = 1 zijn
λ1 = 3.35 + 0.65i, λ2 = −2.87, λ3 = −6.6 en λ4 = 3.35 − 0.65i
De eigenwaarden hebben we berekend met behulp van Maple, zie Appendix C.1.
Nu willen we voor elke snelheid de bijbehorende eigenwaarden plotten, zodat we een beeld krijgen van
wat de eigenwaarden doen. Dit doen we met behulp van Matlab, zie Appendix D.1. De plot van de
eigenwaarde staat in Figuur 6. Dit komt overeen met de plot van de eigenwaarden in [1].
De fiets is stabiel als de snelheid tussen de 4.5 [m/s] en de 6 [m/s] is.
24
Figuur 6: Plot eigenwaarden tegenover snelheid
6.2
Toepassen Routh Stability Criteria
In het model dat we nu hebben zien A, B, C, D en E er zo uit
A
=
32.15
B
=
88.75v
C
= −1151.66 + 95.89v 2
D
= −584.17v + 56.65v 3
E
=
7284.94 − 203.98v 2
De Routh-determinant is X = BCD − ADD − EBB. Uitrekenen met Maple geeft:
X = (−8.65 · 106 )v 2 − (7.03 · 106 )v 4 + (3.79 · 105 )v 6
Omdat A positief is, onafhankelijk van de snelheid, hoeven we dus enkel het geval te onderzoeken dat
A, B, C, D, E en X positief zijn. Dit reken we uit met Maple, zie Appendix C.2.
We vinden hetzelfde stabileitsgebied als met Matlab, waar we de eigenwaarden gebruiken. A, B, C, D, E
en X zijn functies van de snelheid v. Aan de formules kunnen we zien dat dit vrij ’eenvoudige’ functies zijn.
Door ze in ´e´en plaatje te plotten kunnen we ook zien wat het stabiliteitsgebied is, en welke van de functies
daar invloed op hebben. Wanneer we weten welke functies veel invloed hebben op het stabiliteitsgebied
kunnen we kijken naar het algemene model, welke co¨efficienten zorgen voor de formule.
In figuur 7 plotten we de functies A, B, C, D, E en X en we zien dan het stabiliteitsgebied en we zien ook
door welke functies het stabiliteitsgebied begrensd wordt, namelijk door X en door E.
25
Figuur 7: Maple plot van de functies A, B, C, D, E en X.
In Figuur 8 is te zien dat de Routh determinant X, de roze lijn, inderdaad negatief is buiten het stabiliteitsgebied en dus positief binnen het stabiliteitsgebied.
Figuur 8: Maple plot van de Routh determinant X.
7
Een eerste conclusie
Zowel bij de plot van de eigenwaarden als bij het Routh Stability Criteria hebben we gezien dat de fiets,
met het huidige model, stabiel is als de snelheid tussen de 4,5 [m/s] en 6 [m/s] is. Bovendien hebben we
gezien dat het stabiliteitsgebied begrensd wordt door de vergelijkingen E en X uit het Routh Stability
Criteria.
26
Wanneer we nu kijken naar het algemene model, zie (46), en daar de bijbehorende Routh vergelijkingen
van berekenen, zie (47), dan zien we dat E de volgende vergelijking is:
2
E = g((k11 p22 − k12 p12 )v 2 + g(k11 k22 − k12
))
De Routh determinant X wordt een grote term die op dit moment nog geen extra(/nuttige) informatie
toevoegt, deze staat in Appendix E. Maar E is een veel eenvoudigere uitdrukking, en we weten dat deze
het stabiliteitsgebied begrenst. De uitdrukking voor E is een tweede graads vergelijking met nulpunten:
s
v0 =
2 )g
(k11 k22 − k12
(k11 p22 − k12 p12 )
(49)
Merk op dat we alleen naar de positieve wortel kijken omdat er geen negatieve snelheid is, maar we
spreken van nulpunten om dat de negatieve oplossing er ook is. Wanneer we het stabiliteitsgebied willen
2
vergroten willen we dat |v0 | groter wordt. Dit kan door k11 k22 − k12
groter te maken of k11 p22 − k12 p12
kleiner te maken.
Nu is het interessant om te kijken wat die co¨efficienten betekenen. Hiervoor hebben we de co¨efficienten
a3 , a6 , b3 en b6 nodig zoals in (45), het verband met de matrices K en P zoals beschreven in H.4. We
vinden dan
"
K=
k11
k12
k12
k22
"
P=
#
0
p12
0
p22
"
=
#
−Sx
− pb Sz + f Mv
− pb Sz + f Mv − pb Sz + f Mv cot(α)

#
=
h
0
0
h
p
b Sz
Sx +
+ f Mv +
i
Iwa
Iwv
rv + ra i
Iwv
rv cot(α)
1
b

1
b

Wanneer we de betekenis invullen in (49) krijgen we:
v
u
u
u
v0 = t
p
2 S
+
f
M
cot(α)
−
−
S
+
f
M
g
v
v
b z
b z
h
i h
(−Sx ) pb Sz + f Mv + Irwv
cot(α) 1b − − pb Sz + f Mv
+
Sx + Irwv
v
v
(−Sx ) −
p
Iwa
ra
i 1
b
Nu weten we de natuurkundige uitdrukking van de bovengrens van de stabiliteit van de fiets. Het is
alleen nog niet direct duidelijk welke parameter nu direct invloed uitoefent op de stabiliteit. Daarom
gaan we eerst kijken naar welke matrixco¨efficient invloed heeft op de stabiliteit. We doen dit door steeds
een matrixco¨efficient te vari¨eren en vervolgens naar het stabiliteitsgebied te kijken. Hierbij moeten we
wel meteen opmerken dat we hiermee enkel wiskundig kijken naar wat er gebeurt en laten we los wat de
betekenissen zijn van de co¨efficienten. In werkelijkheid hangen al deze co¨efficienten af van de parameters
en dus soms ook van elkaar. Maar om een idee te krijgen wat er gebeurt kunnen we dat, met het zojuist
opgemerkte in het achterhoofd, wel doen.
8
Verandering model
We laten nu de betekenis van de co¨efficienten uit de matrix los, we bekijken puur wiskundig wat een
verandering in het model betekent voor het stabiliteitsgebied. Zo verkrijgen we inzicht wat de co¨efficient
voor invloed heeft op het model. Daarna betrekken we weer de betekenis van de co¨efficienten erbij om te
kijken hoe we op een juiste manier wijzigingen in het model aan kunnen brengen.
27
Als eerste gaan we co¨efficienten k11 ,k12 en k22 om beurten vari¨eren. In Maple berekenen we de determinant
die steeds afhangt van v en een variabele in de matrix. Dan passen we het Routh Stability Criteria toe
en plotten we het stabiliteitsgebied. Zie Figuur 9, 10 en 11.
Figuur 9: Hier is k11 als variabele neergezet, de rest is als het
startmodel. In het startmodel is
k11 = −82.5. Wanneer we richting 0 gaan zien we dat het stabiliteitsgebied verder reikt, maar
ook dat het steeds later zal starten.
Figuur 10: Hier is k12 als variabele neergezet, de rest is als
het startmodel. In het startmodel is k12 = −2.08. Wanneer we
k12 steeds kleiner kiezen zal het
startmodel steeds later starten.
Figuur 11: Hier is k22 als variabele neergezet, de rest is als
het startmodel. In het startmodel is k11 = −0.97. Wanneer we
k2 2 kleiner maken zal het stabiliteitsgebied later beginnen en later eindigen.
Door deze plots proberen we eerst k11 te vari¨eren, deze geeft met een kleine aanpassing al een grote
verandering. Met het huidige model is de fiets stabiel als de snelheid tussen de 4.5 en de 6 [m/s] is.
Wanneer we nu co¨efficient k11 van matrix K veranderen van -82.5 naar -75 dan is de fiets stabiel als de
snelheid groter is dan 4.7 [m/s]. Daarmee vergroten we het stabiliteitsgebied.
We bekijken nu wat de betekenis is van k11 en we gaan kijken wat er gebeurt als we die parameters
veranderen. We hebben gezien dat k11 = −Sx . In Figuur 20 lezen we af dat:
Sx = M hf = massa van de fiets incl. berijder × afstand van Ca tot het grondvlak (tot de x-as)
waarbij Ca het zwaartepunt van het achterste gedeelte van de fiets is, inclusief de berijder. De stabiliteit
van de fiets is dus te be¨ınvloeden door de massa van zowel de fiets als de berijder en de verdeling van de
massa over de fiets.
8.1
Resultaten na verandering
De nieuwe situatie veranderen we k11 van 82, 5 naar 75. Het is natuurlijk interessant om te kijken naar
de veranderingen die dat met zich mee brengt voor de stabiliteit van de fiets. Daarom betrekken we nu
weer de parameters van de fiets er bij. Wat betekent het dat we k11 willen veranderen? Dit betekent
dat we Sx willen veranderen, maar Sx zit niet alleen in k11 maar ook in andere co¨efficienten. Bovendien
wordt Sx berekend door M en hf . We pakken dit aan door naar de twee uitersten te kijken, bij de
eerste veranderen we M van 82.5 naar 75 en als tweede veranderen we hf van 1 naar 0.90909, want
82.5 · 0.90909 = 75. We gaan met Matlab steeds kijken naar de lengte van het stabiliteitsgebied. In ons
oorspronkelijke model is de fiets stabiel als de snelheid tussen 4.5 en de 6 [m/s] is, dus de lengte van het
stabiliteitsgebied is dan 1.5.
28
Met behulp van Matlab kijken we wat er gebeurt als we de massa, M , van de fiets incl. berijder veranderen.
In Figuur 12 is te zien wat de lengte van het stabiliteitsgebied is bij elke waarde van M . Zie de Matlabcode
in Appendix D.2.
Figuur 12: Lengte van stabiliteitsgebied per waarde van M
En zo kijken we ook wat er gebeurt als we de hoogte, hf , veranderen. In Figuur 13 is te zien wat de
lengte van het stabiliteitsgebied is bij elke waarde van hf . Zie de Matlabcode in Appendix D.3.
8.2
Conclusie
Het blijkt dat hoe lager de massa van de berijder is, hoe groter het stabiliteitsgebied wordt. Hier kunnen
we niet direct iets mee, want het gewicht van de berijder is niet specifiek voor het ontwerp van de fiets.
Het is wel een test naar de werkelijkheid, want de fiets met zware berijder zal al sneller vallen als de fiets
schuin komt te staan dan bij een fiets met lichte berijder.
Daarnaast zien we dat hoe kleiner hf is, dus hoe lager het zwaartepunt bij de grond is, hoe stabieler de
fiets is. Dit hangt wel samen met het ontwerp van de fiets. Door het zwaartepunt lager proberen te krijgen
kan het stabiliteitsgebied misschien vergroot worden. Let op dat dit misschien kan, andere parameters
zullen veranderen als het zwaartepunt verlaagd wordt en dit zou de vergroting van het stabiliteitsgebied
kunnen tegengaan.
Hier hebben we een referentiepunt naar de werkelijkheid, want bij kinderfietsen zit het zwaartepunt lager
dan bij een fiets voor volwassene. Kinderfietsen zijn dus stabieler dan ’grote’ fietsen, want kinderen, die
nog niet altijd heel goed kunnen fietsen, kunnen het zo makkelijker leren omdat ze minder snel uit balans
zijn.
29
Figuur 13: Lengte van stabiliteitsgebied per waarde van hf
9
Gyroscopisch effect
Over het algemeen wordt aangenomen dat een fiets stabiel is door het gyroscopisch effect. Er zijn berichten
dat zonder het gyroscopisch effect de fiets stabiel kan zijn, zo ook in [2]. Voordat we gaan kijken of het
gyroscopisch effect wel of niet zorgt voor stabiliteit, gaan we eerst bekijken wat het gyroscopisch effect
nu eigenlijk is.
9.1
Wat is het gyroscopisch effect
Om een idee te krijgen wat het gyroscopisch effect is, beschrijven we een experiment om te laten zien wat
het gyroscopisch effect doet3 . Wanneer we een wiel laten draaien oefenen we een kracht, F , uit op het
wiel. Het draaimoment, de torsie, is dan t = r × F , waarbij r de straal van het wiel is. Het draaimoment
staat loodrecht op r en F , door de rechterhandregel weten we welke richting. Ter illustatie zie Figuur
14A.
Hoe harder we het wiel laten draaien hoe hoe groter het draaimoment is. Nu laten we het wiel niet zomaar
draaien, we doen een staaf door het middelpunt van het wiel en aan het uiteinde van de staaf maken we
een touw vast dat aan het plafond hangt. Als we het wiel rechtop houden dan werkt de zwaartekracht
naar beneden en de arm is de lengte van de staaf tussen het wiel en het touw. Wanneer we het wiel
loslaten zal het wiel opzij vallen en gaan slingeren, zie Figuur 14B. Nu geven we eerst het wiel een flinke
draai, zodat het draaimoment heel groot wordt en dan laten we het wiel los, dan valt het wiel niet naar
beneden maar gaat het rondjes draaien om het touw heen. Dit effect, dat een draaiend wiel om zijn as
gaat draaien noemen we het gyroscopisch effect, zie Figuur 14C.
3 https://www.youtube.com/watch?v=ty9QSiVC2g0
30
Figuur 14: uitleg gyroscopisch effect
9.2
Gyroscopisch effect fiets
Wanneer we nu kijken naar het model van de fiets zien we dat het gyroscopische effect in de momenten
zit, zie Figuur 5. De traagheidsmomenten Iwv en Iwa zorgen voor het gyroscopische effect. We hebben
een model gemaakt op basis van alle parameters, dat betekent dat we nu kunnen kijken wat er gebeurt
als er geen gyroscopisch effect zou zijn, dat wil zeggen Iwv = 0 en Iwa = 0.
Allereerst kijken we naar wat er gebeurt met de plot van de eigenwaarden. Dan zien we in Figuur 15 dat er
geen snelheid v is zodat alle eigenwaarden negatief zijn. Maar we zien ook dat de positieve eigenwaarden
wel heel dicht bij 0 komen. Met andere co¨efficienten worden de eigenwaarden misschien wel negatief.
Figuur 15: Eigenwaarde plot zonder gyroscopisch effect
En wanneer we naar het Routh Stability Criteria kijken krijgen we de volgende vergelijkingen:
A(v)
=
32.695
B(v)
=
89.023v
C(v)
=
49.49v 2 − 1157.5
D(v)
=
−586.06v
E(v)
=
7295.1
X(v)
= −2.58204177 · 106 v 4 − 8.6527525 · 106 v 2
(50)
Hieruit kunnen we afleiden dat voor stabiliteit moet gelden dat A, B, C, D, E en X positief zijn, omdat
31
A en E beide constant zijn, dus onafhankelijk van v, en positief zijn. In Figuur 16 zien we dat de blauwe
lijn, D nooit positief is bij een positieve snelheid, dit kunnen we ook zien aan de vergelijking want D is
lineair en gaat door de oorsprong met negatieve richting.
Figuur 16: Routh Stability Criteria zonder gyroscopisch effect
Ook de lijn X komt nooit positief uit, er is ´e´en re¨eel nulpunt, namelijk v = 0. Verder is de functie
negatief, ook te zien in Figuur 17.
Figuur 17: Routh Stability Criteria zonder gyroscopisch effect
32
10
Eindconclusie
Wat hebben we nu met dit alles bereikt? We hebben bereikt dat we nu een aantal conclusies kunnen
trekken op basis van het onderzoek en de resultaten.
Als eerste merken we op dat we het model voor de voorwaartse beweging van de fiets afgeleid hebben, dus
van elke co¨efficient in het model weten wat de herkomst is en dus wat de samenstelling van fietsparameters
is. Dat wil zeggen dat we nu voor elke fiets, waarvan we de fietsparameters hebben, een model kunnen
opstellen om te kijken naar de stabiliteit van de fiets.
Daarnaast kunnen we aan de hand van het model de rand van het stabiliteitsgebeid vastleggen, namelijk
deze fiets is stabiel tussen de 4,5 en de 6 [m/s]. En aan de hand daarvan kunen we bepalen welke
fietsparameters de stabiliteit beperken. Bij de fiets die wij bekeken hebben bleek het stabiliteitsgebied
onder andere beperkt te zijn door de hoogte van het zwaartepunt van de fiets en de massa van de berijder.
Daarvan hebben we gezegd dat deze gelijk een referentie naar de werkelijkheid vormen, want wanneer de
berijder zwaarder is zal het zwaartepunt sneller verplaatsen en de potentiele energie groter worden dan
wanneer de berijder lichter is. En kinderfietsjes zijn sneller stabiel dan fietsen voor volwassenen.
En we kunnen concluderen, in tegenstelling tot [2], het gyroscopisch effect voor de stadsfiets wel degelijk
essentieel is om stabiliteit te krijgen. Zonder het gyroscopisch effect kunnen we geen stabiliteitsgebied
aanwijzen voor deze fiets, maar het is niet onwaarschijnlijk dat met andere fietsparameters er wel een
stabiliteitsgebied aan te wijzen is.
33
34
Referenties
[1] Ir. B.D. Herfkens and vertaling door J.D.G. Kooijman. De stabiliteit van het rijwiel. -, 1949, vertaling
2006.
[2] J.D.G. Kooijman, J.P. Meijaard, J.M. Papadopoulos, A. Ruina, and A.L. Schwab. A bicycle can be
self-stable without gyroscopic or caster effects. Science Magazine, 332:339–442, 2011.
[3] A.L. Schwab and J.P. Meijaard. De stabiele fiets. Nederlands tijdschrift voor Natuurkunde, 74:24–28,
2008.
35
36
A
Parameters fiets
Figuur 18: Figuur 3 uit [1]
Figuur 19: Figuur 11 uit [1]
37
B
Betekenis co¨
efficienten
We nemen de coefficienten aan zoals beschreven in artikel [1]:
Figuur 20: Betekenis en waarden coefficienten, uit [1]
38
C
C.1
Maple
Karakteristieke vergelijking berekenen
Figuur 21: Karakteristieke vergelijking en eigenwaarden
C.2
Stabiliteitsgebied Routh Stability Criteria
Figuur 22: Stabiliteitsgebeid Routh Stability Criteria m.b.v. Maple
39
C.3
Lagrangiaanse vergelijking voor Figuur 23: Epsilon termen Lagrangiaanse vergelijking berekenen
40
D
D.1
Matlab code
Startmodel
Het eerste deel van de code spreek voor zich, de parameters worden gedefinieerd, de matrices ingevoerd,
de determinant berekend en vervolgens het reele en imaginaire deel van de eigenwaarden. Daarna wordt
gezocht wanneer de eigenwaarden allemaal negatief zijn. De eerste en laatste keer dat dat gebeurt wordt
uitgevoerd op het scherm en de eigenwaarden worden geplot tegenover de snelheid, zie Figuur 6.
clear all;
h = 0.05;
v = [0:h:10];
n = length(v);
ReS = zeros(4,n);
ImS = zeros(4,n);
syms w;
stab = zeros(n,1);
g = 9.81;
%stapgrootte
%vector met snelheden van 0 tot 10 en stapgrootte h
%lengte van de sneheidsvector
%oplossingsmatrix reele deel
%oplossingsmatrix imaginaire deel
%variabele lambda
%stabiliteitsvector
%gravitatie
for i=1:n
M = [98, 1.77; 1.77, 0.36];
C = [0, 26.5; -0.79, 1.37];
K = [-82.5, -2.08; -2.08, -0.97];
P = [0, 71.71; 0, 2.06];
A = M*w^2 + C*v(i)*w + K*g + P*v(i)^2;
Determinant = det(A);
%determinant berekenen
EwVec = solve(Determinant == 0);
%oplossingen determinant opslaan in vector
AantEw = length(EwVec);
%aantal oplossingen tellen
for j = 1:AantEw
ReS(j,i) = real(EwVec(j));
%voor elke opl het reele en
ImS(j,i) = imag(EwVec(j));
%imaginaire deel opslaan
end
if ReS(1,i) < 0 && ReS(2,i) < 0 && ReS(3,i) < 0 && ReS(4,i) < 0
stab(i) = 1;
%vector geeft 1 als alle ew < 0
end
end;
c = find(stab, 1, ’first’);
%c is het eerste niet-nul element
d = find(stab, 1, ’last’);
%d is het laatste niet-nul element
StartStab = roundn(v(c),-1);
%begin stabiliteitsgebied
EindStab = roundn(v(d),-1);
%eind stbailiteitsgebeid
disp([’De fiets is stabiel als de snelheid tussen:’])
disp([’v = ’,num2str(StartStab),’[m/s] en v = ’,num2str(EindStab), ’[m/s].’]);
Y = sort(ImS,1);
%imaginaire deel plotten
plot(v,Y(4,:),’k--’)
hold on;
Z = zeros(1,n);
%nullijn plotten
plot(v,Z,’k’);
hold on
X = sort(ReS,1);
%reele deel plotten
plot(v,X(1,:),’k’,v,X(2,:),’k’,v,X(3,:),’k’,v,X(4,:),’k’,’LineWidth’,1.5)
axis([0 10 -10 10])
xlabel(’snelheid van de fiets in [m/s]’);
ylabel(’Rele en Imaginaire deel eigenwaarden’);
hold off;
grid on;
41
D.2
M varieren
clear all;
eindt = 40;
StartStab = zeros(1,eindt);
EindStab = zeros(1,eindt);
LengteStabGebied = zeros(1,eindt);
%begin stabiliteitsgebied
rv = 0.35; ra = 0.35; b
= 1.17; p
= 0.08;
a
=0.12; uf = 0.33; j
= 1.00;
k
= 0.52;
d
= 0.08; ww
= 0.03; f
= 0.04;
alfa = 65*pi/180;
g
= 9.81; hf = 1.00; Mv = 5.55;
Mwa = 2.65;
Mwv = 2.65; Iwa = 0.245; Iwv = 0.245;
Ix = 98;
Iz = 14.7; Ixz = 23.5; Ias = 0.195;
Ian = -0.04;
Is = 0.05; Sz = 27.2;
Ss = 0.005;
Delta = Ias*cot(alfa) + Ian + k*f*Mv;
Gamma = Ias - Ian*cot(alfa) +j*f*Mv;
for t = 1:eindt
h = 0.1;
v = [0:h:10];
n = length(v);
syms w;
stab = zeros(1,n);
ReS = zeros(4,n);
ImS = zeros(4,n);
c = 1;
d = 1;
M
Sx
%stapgrootte
%vector met snelheden van 0 tot 10 en stapgrootte h
%lengte van de sneheidsvector
%variabele w (eigenwaarde)
%oplossingsmatrix reele deel
%oplossingsmatrix imaginaire deel
= 70 + 0.5*t;
= M*hf;
%oorspronkelijk 82.5
for i=1:n
a1 = Ix;
a2 = 0;
a3 = -g*Sx;
a4 = (p/b)*Ixz + Delta;
a5 = (v(i)/b)*(Ixz + p*Sx + (b+p)/rv*Iwv + (p/b)*Iwa);
a6 = (Sx + Iwv/rv + Iwa/ra)*((v(i)^2)/b) - g*((p/b)*Sz + f*Mv);
b1 = (p/b)*Ixz + Delta;
b2 = -v(i)*(((b+p)/b)*(Iwv/rv) +(p/b)*(Iwa/ra));
b3 = -g*((p/b)*Sz + f*Mv);
b4 = Ias/(sin(alfa)^2) + ((p^2)/(b^2))*Iz +2*p*Gamma/b;
b5 = v(i)/b*((p/b)*Iz + p*((p/b)*Sz + f*Mv) + Gamma);
b6 = ((p/b)*Sz + f*Mv + (Iwv/rv)*cot(alfa))*((v(i)^2)/b)-g*((p/b)*Sz + f*Mv)*cot(alfa);
%gevolg overgenomen
A = a1*b4 - b1*a4;
B = a1*b5 - b1*a5 C = a1*b6 + a3*b4 D = a3*b5 - b2*a6 E = a3*b6 - b3*a6;
van blz. 16 van artikel uit 1949
b2*a4;
b1*a6 - b2*a5 - b3*a4;
b3*a5;
%4e orde differentiaalvgl oplossen
vgl = A*w^4 + B*w^3 + C*w^2 + D*w + E;
Z = solve(vgl == 0);
AantEw = length(Z);
EwVec = Z; %sort(Z);
for j = 1:AantEw
42
%aantal oplossingen tellen
ReS(j,i) = real(EwVec(j));
ImS(j,i) = abs(imag(EwVec(j)));
%voor elke opl het reele en
%imaginaire deel opslaan
end
if ReS(1,i) < 0 && ReS(2,i) < 0 && ReS(3,i) < 0 && ReS(4,i) < 0
stab(i) = 1;
%vector geeft 1 als alle ew < 0
end
end;
if stab == zeros(1,n)
LengteStabGebied(t) = 0;
else
c = find(stab, 1, ’first’);
%c is het eerste niet-nul element
d = find(stab, 1, ’last’);
%d is het laatste niet-nul element
StartStab(t) = roundn(v(c),-1);
%begin stabiliteitsgebied
EindStab(t) = roundn(v(d),-1);
%eind stbailiteitsgebeid
LengteStabGebied(t) = EindStab(t)-StartStab(t);
end
end
figure
s1 = length(LengteStabGebied);
s = 1:s1;
plot(s,LengteStabGebied)
xlabel(’Waarde van M in [kg] (oorspronkelijk 82.5 kg is bij x=25)’);
ylabel(’Lengte stabiliteitsgebied’);
grid on;
D.3
hf varieren
clear all;
eindt = 20;
StartStab = zeros(1,eindt);
EindStab = zeros(1,eindt);
LengteStabGebied = zeros(1,eindt);
%begin stabiliteitsgebied
rv = 0.35; ra = 0.35; b
= 1.17; p
= 0.08;
a
=0.12; uf = 0.33; j
= 1.00;
k
= 0.52;
d
= 0.08; ww
= 0.03; f
= 0.04;
alfa = 65*pi/180;
g
= 9.81; M
= 82.5; Mv = 5.55;
Mwa = 2.65;
Mwv = 2.65; Iwa = 0.245; Iwv = 0.245;
Ix = 98;
Iz = 14.7; Ixz = 23.5; Ias = 0.195;
Ian = -0.04;
Is = 0.05; Sz = 27.2;
Ss = 0.005;
Delta = Ias*cot(alfa) + Ian + k*f*Mv;
Gamma = Ias - Ian*cot(alfa) +j*f*Mv;
for t = 1:eindt
h = 0.1;
%stapgrootte
v = [0:h:10];
%vector met snelheden van 0 tot 10 en stapgrootte h
n = length(v);
%lengte van de sneheidsvector
syms w;
%variabele w (eigenwaarde)
stab = zeros(1,n);
ReS = zeros(4,n);
%oplossingsmatrix reele deel
ImS = zeros(4,n);
%oplossingsmatrix imaginaire deel
c = 1;
d = 1;
hf
Sx
= 0.1*t; %(oorspronkelijk: 1.00)
= M*hf;
43
for i=1:n
a1 = Ix;
a2 = 0;
a3 = -g*Sx;
a4 = (p/b)*Ixz + Delta;
a5 = (v(i)/b)*(Ixz + p*Sx + (b+p)/rv*Iwv + (p/b)*Iwa);
a6 = (Sx + Iwv/rv + Iwa/ra)*((v(i)^2)/b) - g*((p/b)*Sz + f*Mv);
b1 = (p/b)*Ixz + Delta;
b2 = -v(i)*(((b+p)/b)*(Iwv/rv) +(p/b)*(Iwa/ra));
b3 = -g*((p/b)*Sz + f*Mv);
b4 = Ias/(sin(alfa)^2) + ((p^2)/(b^2))*Iz +2*p*Gamma/b;
b5 = v(i)/b*((p/b)*Iz + p*((p/b)*Sz + f*Mv) + Gamma);
b6 = ((p/b)*Sz + f*Mv + (Iwv/rv)*cot(alfa))*((v(i)^2)/b)-g*((p/b)*Sz + f*Mv)*cot(alfa);
%gevolg overgenomen
A = a1*b4 - b1*a4;
B = a1*b5 - b1*a5 C = a1*b6 + a3*b4 D = a3*b5 - b2*a6 E = a3*b6 - b3*a6;
van blz. 16 van artikel uit 1949
b2*a4;
b1*a6 - b2*a5 - b3*a4;
b3*a5;
%4e orde differentiaalvgl oplossen
vgl = A*w^4 + B*w^3 + C*w^2 + D*w + E;
Z = solve(vgl == 0);
AantEw = length(Z);
%aantal oplossingen tellen
EwVec = Z; %sort(Z);
for j = 1:AantEw
ReS(j,i) = real(EwVec(j));
%voor elke opl het reele en
ImS(j,i) = abs(imag(EwVec(j)));
%imaginaire deel opslaan
end
if ReS(1,i) < 0 && ReS(2,i) < 0 && ReS(3,i) < 0 && ReS(4,i) < 0
stab(i) = 1;
%vector geeft 1 als alle ew < 0
end
end;
if stab == zeros(1,n)
LengteStabGebied(t) = 0;
else
c = find(stab, 1, ’first’);
%c is het eerste niet-nul element
d = find(stab, 1, ’last’);
%d is het laatste niet-nul element
StartStab(t) = roundn(v(c),-1);
%begin stabiliteitsgebied
EindStab(t) = roundn(v(d),-1);
%eind stbailiteitsgebeid
LengteStabGebied(t) = EindStab(t)-StartStab(t);
end
end
figure
s1 = length(LengteStabGebied);
s = 1:s1;
plot(s,LengteStabGebied)
xlabel(’Waarde van hf in [m] (oorspronkelijk 1.00 m)’);
ylabel(’Lengte stabiliteitsgebied’);
grid on;
44
E
X
Routh determinant X
=
c212 c221 m12 p12 + c12 c321 m12 p12 − c12 c221 c22 m11 p12 − c12 c21 m11 m12 p12 p22
+ c12 c21 m212 p212 − c221 m11 m12 p12 p22 − c221 m11 m22 p212 + 2c221 m212 p212
+ c21 c22 m211 p12 p22 − c21 c22 m11 m12 p212 v 6
+
c312 c21 gk12 m12 + 2c212 c221 gk12 m12 − c212 c21 c22 gk11 m12 − c212 c21 c22 gk12 m11
− c212 gk11 m212 p22 − c212 gk12 m11 m12 p22 + 2c212 gk12 m212 p12 + c12 c321 gk12 m12
− c12 c221 c22 gk11 m12 − c12 c221 c22 gk12 m11 + c12 c21 c222 gk11 m11 − 2c12 c21 gk11 m212 p22
− c12 c21 gk11 m12 m22 p12 − 2c12 c21 gk12 m11 m12 p22 − 2c12 c21 gk12 m11 m22 p12 + 8c12 c21 gk12 m212 p12
− c12 c21 gk22 m11 m12 p12 + 3c12 c22 gk11 m11 m12 p22 − c12 c22 gk11 m212 p12 + c12 c22 gk12 m211 p22
− 3c12 c22 gk12 m11 m12 p12 − c221 gk11 m212 p22 − c221 gk11 m12 m22 p12 − c221 gk12 m11 m12 p22
− 2c221 gk12 m11 m22 p12 + 6c221 gk12 m212 p12 − c221 gk22 m11 m12 p12 + 3c21 c22 gk11 m11 m12 p22
+ 3c21 c22 gk11 m11 m22 p12 − 3c21 c22 gk11 m212 p12 + c21 c22 gk12 m211 p22 − 5c21 c22 gk
12 m11 m12 p12
+ c21 c22 gk22 m211 p12 − 2c222 gk11 m211 p22 + c222 gk11 m11 m12 p12 + c222 gk12 m211 p12 v 4
2
2
−
c212 g 2 k11 k12 m12 m22 − c212 g 2 k11 k22 m212 − c212 g 2 k12
m11 m22 + 4c212 g 2 k12
m212 − c212 g 2 k12 k22 m11 m12
2
2
− 2c12 c21 g 2 k11 k12 m12 m22 − 2c12 c21 g 2 k11 k22 m212 − 2c12 c21 g 2 k12
m11 m22 + 8c12 c21 g 2 k12
m212
2
2 2
2
2
− 2c12 c21 g k12 k22 m11 m12 + c12 c22 g k11 m12 m22 + 3c12 c22 g k11 k12 m11 m22 − 4c12 c22 g k11 k12 m212
2
+ 3c12 c22 g 2 k11 k22 m11 m12 − 4c12 c22 g 2 k12
m11 m12 + c12 c22 g 2 k12 k22 m211 − c221 g 2 k11 k12 m12 m22
2 2
2
2 2 2
2
2
m12 m22
− c21 g k11 k22 m12 − c21 g k12 m11 m22 + 4c221 g 2 k12
m212 − c221 g 2 k12 k22 m11 m12 + c21 c22 g 2 k11
2
2
2
2
2 2
+ 3c21 c22 g k11 k12 m11 m22 − 4c21 c22 g k11 k12 m12 + 3c21 c22 g k11 k22 m11 m12 − 4c21 c22 g k12 m11 m12
2
2 2 2
2 2 2
2
2 2
2 2
2
+ c21 c22 g 2 k12 k
22 m11 − 2c22 g k11 m11 m22 + c22 g k11 m12 + 2c22 g k11 k12 m11 m12 − 2c22 g k11 k22 m11
2
+ c222 g 2 k12
m211 v 2
45

Download Report