Study of self excited periodic oscillations of a woodpecker toy

IN

DEGREE PROJECT TECHNOLOGY, FIRST CYCLE, 15 CREDITS

,

STOCKHOLM SWEDEN 2019

Study of self excited periodic

oscillations of a woodpecker toy

INOM EXAMENSARBETE TEKNIK, GRUNDNIVÅ, 15 HP , STOCKHOLM SVERIGE 2019

Studie av självexciterade

periodiska svängningar hos en

hackspettsleksak

Abstract

Sammanfattning

Innehåll

3.3.1 Lagen om den kinetiska energin . . . 4

3.3.2 Potentiella energin och mekaniska energilagen . . . 6

3.4 Tvång på position och rörelse . . . 7

3.4.1 Holonoma tvång . . . 8 3.4.2 Tvång på hastighet . . . 8 3.5 Lagrangesk mekanik . . . 9 4 Ordinära differentialekvationer 10 4.1 Explicit ODE . . . 10 4.2 Implicit ODE . . . 11 4.3 DAE . . . 12 5 Numeriska metoder 16 5.1 Fixa punkter och fixpunktsiteration . . . 17

5.2 Kontraktiva sekvenser och Picards algoritm . . . 17

5.3 Numerisk lösning av DAEs . . . 19

6 Diskontinuerliga system och växlingsalgoritmer 19 Modeller och algoritmer 21 7 Simpel modell 21 7.1 Systembeskrivning och härledning . . . 22

7.1.1 Fritt fall . . . 23

7.1.2 Fast i staven . . . 25

7.1.3 Händelser och fasövergångar . . . 26

7.2 Algoritm . . . 27

7.3 Implementation i Matlab . . . 28

8 Komplex modell 33

8.1 Systembeskrivning . . . 34

8.1.1 Stadie I - fritt fall . . . 35

8.1.2 Stadie II - ring fast i nedre läge . . . 36

8.1.3 Stadie III - ring fast i övre läge . . . 36

8.1.4 Händelser och fasövergångar . . . 37

8.2 Algoritm . . . 38

8.3 Implementation i Matlab . . . 39

8.4 Beräkning av perioden . . . 45

Resultat och diskussion 47 9 Resultat simpel modell 47 9.1 Simulationsresultat och periodberäkning . . . 47

10 Resultat komplex modell 50 10.1 Simulationsresultat och periodberäkning . . . 50

11 Diskussion 53 11.1 Diskussion av resultaten . . . 53

11.2 Slutsatser . . . 58

1

Inledning

Varför är det intressant att modellera och simulera en hackspettleksak? Ett svar på den frågan är att systemet är komplext och uppvisar en dynamik vi även hittar hos växlingslådor i bilar och borrmaskiner i vår vardag, näm-ligen oscillationer. Hackspettsleksaken är ett komplext system, men inte tillräckligt komplicerat för att vi inte skall kunna modellera dess rörelse. Analys av simulationsresultaten kan på så sätt ge oss större insikt generellt i hur oscillativa dynamiska med uppför sig. Speciellt för hackspettsleksaken är dess självsvängning kring ett jämviktsläge, där rörelsen inkluderar stö-tar och låsningar vilket leder till impulser i hackspettens rotationshastighet. Denna diskontinuerliga dynamik måste delas upp i kontinuerliga intervall som kan modelleras med system av differentialekvationer, och speciellt lås-ningarna, sätter krav på de numeriska metoderna vi använder oss av. I den här studien ges en inkörsport till detta intressanta område av beräknings-mekaniken: Hur komplexa mekaniska system modelleras, implementeras och löses numeriskt.

2

Problemformulering

Teoretisk bakgrund

I de kommande avsnitten behandlas de teoretiska grunderna som studien vilar på, från klassisk mekanik och system av ordinära differentialekvationer till numeriska metoder för att lösa de problem som uppstår. Allt avslutas med ett avsnitt om diskontinuerliga system och växlingsalgoritmer.

3

Stelkroppsmekanik

Följande avsnitt inleds med definitionen av en stel kropp och hur ett dyna-miskt system kan modelleras såsom flera stela kroppar. Så kallade holonoma

tvång kommer behandlas vilka sätter villkor på positionen för systemets

rö-relse och kallas även geometriska restriktioner. De ekvationer som beskriver systemets dynamik ges entydigt av Lagranges ekvationer, vilken har två ek-vivalenta framställningar: Lagranges ekvationer av första respektive andra typen. När systemet är konservativt är ofta första framställningen mer an-vändbar.

3.1 Stel kropp

En stel kropp definieras som en massbelagd domän i vilket avståndet mellan två godtyckliga punkter A och B är konstant (se [6] s. 61). Antagandet om det konstanta avståndet innebär att en stel kropp är oderformerbar, vilket i verkligeheten är en bra approximation av kroppar med försumbart små deformationer.

En stel kropps läge i rummet bestäms vanligtvis av den linjära positionen tillsammans med rotationsläget. Den linjära positionen bestäms av läget av en fix punkt på kroppen relativt ett rumsfixt koordinatsystem där den kroppsfixa punkten oftast väljs i masscentrum. Om ett koordinatsystem fixe-ras i kroppen så kan rotationsläget beskrivas av kroppens axlar i förhållande till det rumsfixa koordinatsystemet. Ett elegant sätt att beskriva rotationslä-get på är med hjälp av Eulers vinklar som visas nedan i Figur 1. Tillsammans med de kartesiska rumskoordinaterna ges en kropps läge i tre dimensioner av (x, y, z, θ, φ, ψ). För en kropp som utför plan rörelse i xy-planet, så ges kroppens läge av (x, y, θ), där θ är vinkeln mellan den rumsfixa respektive den kroppsfixa x-axeln.

Figur 1: Eulers vinklar 3.2 Flerkroppssystem

Ett flerkroppsystem består av flera sammanlänkade och/eller fria kroppar. För att illustrera skillanden mellan begreppen sammanlänkande och fria kroppar ges i Figur 2 nedan ett exempel på hur två fria partiklar samman-länkas av en masslös led (notera att partiklar antas inte ha rotationsläge, enbart linjär position).

Läget för ett flerkroppssystem beskrivs som alla ingeånde kroppars lägen. För ett flerkroppssytem av n fria kroppar så beskrivs alltså läget av 6n

fri-hetsgrader och detta gäller för alla flerkroppssystem i tre dimensioner. När

en eller flera kroppar är sammanlänkade så kan koordinatframställningen reduceras till färre koordinater i och med att kropparnas rörelse relativt varandra begränsas. I exemplet ovan såg vi att en led mellan partiklarna reducerade antalet frihetsgrader från fyra till tre, där det nya systemet be-traktas som en kropp som utför planrörelse. De kartesiska koordinaterna är placerade i masscentrum för kroppen och rotationsläget framgår av Fi-gur 2. Det gemensamma masscentrumet för två eller fler masscentrum kan beräknas med formeln

rG = 1 m n X k=1 mkrk

(a) Fria partiklar _{(b) Sammanlänkade partiklar} Figur 2: I det sammanlänkade tillståndet är partiklarnas rörelse begränsad relativt varandra till skillnad från det fria tillståndet.

masscentrum. Det är med andra ord ett viktat medelvärde och gäller för alla dimensioner upp till tre.

Ett viktigt specialfall fås när antalet masscentrum är två, det vill säga anta-let kroppar är två, och origo placeras i den ena eller andra kroppens mass-centrum. Resultatet blir avståndet till det gemensamma masscentrumet från den kropp vi sätter origo i. Om kropparna betecknas S1 och S2, det

gemen-samma masscentrumet betecknas G och origo sätts i masscentrum för S1 så

ges avståndet av

rG= r m2

m1+ m2

där r är avståndet mellan S1och S2 och m1, m2 är deras massor, respektive.

Detta kallas för Barycentrum. Genom att cykliskt permutera indexet fås avståndet från S2 till G.

3.3 Konservativa system

Vad som kännetecknar konservativa system är att de bevarar energin. Man talar om konservativa krafter vars uträttade arbete lagras i systemet som

potentiell energi och som sedan kan omvandlas till kinetisk energi, alltså

rörelse. För ett dynamiskt system som är konservativt gäller att summan av den potentiella respektive kinetiska energin, den så kallade mekaniska

energin, bevaras. Härunder sammanfattas energilagarna och ett par viktiga

exempel på konservativa krafter redogörs för.

3.3.1 Lagen om den kinetiska energin

är kroppens linjära förflyttning per tidsenhet och spinnet är dess rotations-hastighet runt en given axel. För planrörelse gäller att spinnet alltid sker runt axeln vinkelrät mot planet. Vidare kommer enbart fallet för planrörel-se behandlas. Notera att kinetiska energin är ett realtivt mått som beror av koordinatsystemet och är alltså inte invariant.

Translations- och rotationshastigheten för en kropp ger två separata energi-bidrag som tillsammans bildar den totala kinetiska energin enligt formeln

T = Ttrans+ Trot= 1

2mv2+ 1 2IG˙θ

där T representerar den totala kinetiska energin och m, v, ˙θ är kroppens massa, translationshastighet och rotationshastighet, respektive. IG är krop-pens tröghetsmoment mätt relativt masscentrum för kroppen, vilket är ett krav för att framställningen skall vara korrekt.

För ett flerkroppssystem som består av flera fria och/eller sammanlänkande kroppar ges den totala kinetiska energin för systemet som summan av alla enskilda kroppars energier, alltså

Ttot = X

i

Ti

där indexet i itereras över alla ingående kroppar i systemet och där Ttot betecknar systemets totala kinetiska energi.

Innan vi kan formulera lagen om den kinetiska energin måste även begreppet arbete definieras. Om en kraft uträttar ett arbete på en kropp så ändras dess kinetiska energi, alltså är arbete en storhet som mäts i energi. Kraften uträttar enbart ett arbete om den verkar tangentiellt med rörelseriktningen; krafter vinkelräta mot rörelseriktningen utträttar inget arbete. Matematiskt kan vi skriva detta på integralform som

W =

Z

C

F · dr

där C representerar banan vi väljer att integrera utefter. Lagen om den kinetiska energin skrives matematiskt som

W = T2− T1

där T1, T2 representerar den kinetiska energin före respektive efter arbetet

utförts. Med ord säger lagen att det uträttade arbetet på en kropp är lika

med skillnanden av kroppens rörelseenergi före och efter. Notera att arbetet

3.3.2 Potentiella energin och mekaniska energilagen

Eftersom den totala energin i ett konservativt system bevaras, så betyder det att energi måste lagras när den kinetiska energin ändras, annars skulle ju inte den totala energin bevaras! Denna lagrade energi kallas för den

potenti-ella energin. Beroende på vilka konservativa krafter som uträttar arbete på

en kropp, så ser uttrycken för den potentiella energin olika ut. Till skillnad från den kinetiska energin så kan denna inte konstrueras på ett entydigt sätt. Nedan följer potentialerna för två vanligt förekommande konservativa krafter: Tyngdkraften och fjäderkraften.

Innan vi redogör för dessa potentialer så måste två viktiga egenskaper hos konservativa krafter belysas som uttrycks av arbetsintegralen nedan

W = Z C F · dr= Z (2) (1) −∇V · dr= V₁− V₂

där integrationsgränserna (1) och (2) representerar ändarna på banan vi väljer att integrera utefter. Vi tolkar detta resultat som att det fordrade

arbetet för att förflytta en kropp mellan två punkter är oberoende av vägen den tar. Med andra ord kan vi beräkna arbetet en konservativ kraft utför

bara genom att veta potentialens värde i ändpunkterna. Detta är den förs-ta egenskapen. Den andra egenskapen är att den konservativa kraften kan skrivas på formen

F = −∇V

där V betecknar den korresponderande potentialen. Genom att känna till potentialen kan alltså kraften härledas, vilket i sin tur kan ge oss dynamiken för ett system.

Tyngdkraften

För en kropp med massan m som påverkas av jordens gravitation vid ytan, alltså en jordradie (6371 km) från masscentrum, så fås följande potential: Vi har dr = dxex+ dyey+ dzez och F = −mgez vilket ger

V(r) = − Z r r1 F · dr= − Z r r1 −mgez· dr= Z z z1 mgdz= mgz + C

där g ≈ 9.82m/s2 _{är jordens tyngdacceleration vid ytan och z är höjden}

över ytan.

Fjäderkraften

För en fjäder som antingen töjs ut longitudinellt eller böjs i något plan så uppkommer en fjäderkraft proportionell mot förskjutningen som ges av

Hoo-kes lag, nämligen F = −kxex. Med dr = dxex+dyey+dzezfås potentialen

för fjäderkraften till V = − Z r r1 F · dr= − Z x x1 −kxdx= 1 2kx2+ C

Konstanten kan väljas fritt och om C sätts till 0 kan potentialen skrivas

V = 1

2kx2

där k betecknar en proportionalitetskonstant, fjäderkonstanten, och x är för-skjutningen från det ospända tillståndet.

Härledningarna för potentialerna ovan är nedkortade och förklaras mer ut-förligt i [6] ss. 253-256

3.4 Tvång på position och rörelse

Ett system, e.g. en kropp, som är fri har inga tvång och kan utföra otvungen

rörelse som beskrivs matematiskt medelst Newtons andra lag m¨r = fa(t, r)

där m är kroppens massa, r ortsvektorn och fa kraftresultanten på krop-pen från yttre verkande krafter. Denna ekvation kallas rörelseekvationen och när tvång införs på position och/eller rörelsen så ändras denna beskrivning vilket vi kommer bekanta oss med nedan. Begreppet tvång innebär här be-gränsningar på kroppens position eller rörelse och beskrivs matematiskt av algebraiska ekvationer kallade bivillkor. I det allmänna fallet för ett system bestående av flera kroppar tar rörelseekvationerna formen

M¨p = fa(t, p, ˙p)

där M är massmatrisen, p en kolumnvektor med de generaliserade

koordi-naterna och fa en vektorvärd funktion med kraftresultanter. Notera att fa

3.4.1 Holonoma tvång

Holonoma tvång är inget annat än tvång på position. Låt oss betrakta ett

godtyckligt dynamiskt system med sådana tvång. Ett enkelt exempel är en partikel som endast tillåts röra sig inom en kub. I detta fall så utgör de sex planen som avgränsar kuben holonoma tvång på partikeln. Om partikeln rör sig innanför så gäller Newtons andra lag som vi redan beskrivit den, men om partikeln istället rör sig på sidorna eller kanterna så aktiveras de holonoma tvången varigenom antalet frihetsgrader reduceras. I själva verket ändras därmed rörelseekvationen och måste kompletteras med så kallade

kontaktkrafterför att tillfredsställa dessa tvång (så att partikeln inte rymmer

sin väg ut genom kuben). Den nya rörelseekvationen skrives

m¨p = fa(p) − fc(p, λ)

0 = g(p)

där fc representerar kontaktkraften vilken alltid är parallell med normalen till ytan partikeln färdas på, det vill säga kontaktkrafter uträttar inget

ar-bete. λ betecknar Lagrangemultiplikatorn och 0 = g(p) är planets ekvation

som utgör partikelns holonoma tvång som alltid måste uppfyllas sålänge kontaktkraften är skild från noll. Notera att här beror fa, de yttre verkande krafterna, varken av tiden eller (den generaliserade) hastigheten. (Tidsobe-roende system kallas vanligtvis för autonoma).

I allmänhet så gäller för system bestående av flera kroppar med holonoma tvång att rörelseekvationerna skrives

M¨p = fa(t, p, ˙p) − G(p)Tλ 0 = g(p)

där g är en vektorvärdfunktion som representerar samtliga holonoma tvång och G är funktionalmatrisen av g som definieras som

G ≡ ∂g

∂p

3.4.2 Tvång på hastighet

Tvång på hastighet kan enkelt åstadkommas genom att derivera befintliga tvång på position

Genom att integrera båda leden i ekvationen återfås såklart de holonoma tvången. Tvång på hastighet som kan integreras till tvång på position kallas även de för holonoma tvång. I de fall då detta inte är möjligt kallas tvången för icke-holonoma. Ett klassiskt exempel på ett icke-holonomt tvång är ett mynt som rullar på en yta med tvångsvillkoret rullning utan glidning. Det här tvånget ger oss ingen information om det fullständiga läget ty myntet kan rulla längs många vägar, och eftersom rullningsvinkeln, vilket represen-terar läget, beror av vägen finns inget entydigt svar.

3.5 Lagrangesk mekanik

Lagrangesk mekanik är ett alternativ till den newtonska formuleringen av

dynamiken där rörelseekvationerna ges av Newtons andra lag. Denna lag postulerar att kraftresultanten på ett objekt är lika med objektets massa

mul-tiplicerat med dess acceleration och skrives matetmatiskt

F = md

2_r

dt2

I den Lagrangeska mekaniken har krafternas roll i systemet ersatts av de ki-netiska och potentiella energierna. Utifrån dessa kan den centrala storheten i Lagrangesk mekanik definieras enligt

L(p, ˙p, t) ≡ T (p, ˙p, t) − V (p, t)

där L betecknar Lagrangefunktionen. Lagrangesk mekanik kan endast tilläm-pas på system vars tvångsvillkor, om sådana existerar, är holonoma. Den Lagrenska framställningen är baserad på variationskalkyl och Hamiltons va-riationsprincip, där den senare är att betrakta som en naturlag i likhet med Newtons andra lag. Genom variationsprincipen erhålls Lagranges ekvationer

av första typen och motsvarar rörelseekvationerna för konservativa system.

Dessa är bara en omformulering av rörelseekvationerna som fås ur Newtons andra lag och måste ge samma dynamik, men de visar sig mer praktiska att arbeta med i komplexa flerkroppssystem.

Lagranges ekvationer av första typen

Genom att tillämpa variationskalkylen och Hamiltons variationsprincip som postulerar att systemets rörelse mellan tidpunkterna t1 och t2 är sådan att

verkningsintegralen

S=

Z t2

t1

L(p, ˙p, t)dt

∂L ∂pi − d dt _∂L ∂˙pi  = 0

där i itereras över de generaliserade koordinaterna. Dessa ekvationer till-fredsställer alltså kravet att S får ett extremum (verkningsintegralen blir

minimal), alltså beskriver de systemets dynamik.

När systemet även har holonoma tvångsvillkor ges rörelseekvationerna av

∂L ∂pi − d dt _∂L ∂˙pi  + nc X k=1 λk ∂gk ∂pi = 0

där ncbetecknar antalet tvångsvillkor och gk är den k:te komponten i g som representerar samtliga holonoma tvång.

4

Ordinära differentialekvationer

En ordinär differentialekvation (ODE) är en differentialekvation som inne-håller en eller funktioner av en och samma oberoende variabel, e.g. tids-variablen t. Detta är fallet för mekaniska system vars dynamik beskrivs av rörelseekvationer, vilka i själva verket är ODEs, där varje frihetsgrad är tids-beroende. Vi skall skilja på två fall av ODEs i det här avsnittet, explicita och implicita ODEs. Sedan redogörs för vad en differential algebraisk

ek-vation (DAE) är och hur den löses. Eftersom en DAE är en blandning av

explicita och implicita ODEs (och algebraiska ekvationer), så ges samtidigt en genomgång i hur sådana löses.

4.1 Explicit ODE

En explicit ODE är ett samband mellan en beroende variabel, e.g. x(t), och dess derivator och skrives på normalform som

x(n)= F (x(n−1), ..., x, t)

(x(n−1)_{(0), ..., x(0)) = (x}(n−1)

0 , ..., x0)

där F är en funktion av derivatorna till x och tiden t, x(n) _{betecknar den}

n:te derivatan med avseende på tiden och (x(n−1)₀ , ..., x0) är begynnelsevärdet

˙x(t) = x(t), t ∈ R

x(0) = 0

där x(0) är begynnelsevärdet och lösningen ges av

x(t) = et, t ∈ R

Om tillståndet (det vill säga x i exemplet ovan) är en vektorvärdfunktion fås den mer generella formen av ett explicit system av ODEs som

x(n)₁ = F1(x(n−1)1 , ..., x1, t) ... x(n)_k = Fk(x(n−1)_k , ..., xk, t) (x(n−1) i (0), ..., xi(0)) = (x (n−1) 0,i , ..., x0,i)

där x = (x1(t), ..., xk(t)) och i itererar över komponenterna i x.

En explicit ODE där ordningen överstiger 1 kan med fördel transformeras till ett system av explicita ODEs varigenom ordningen för varje ODE reduceras till 1 x(n)= F (x(n−1), ..., x, t) =⇒              ˙x = x1 ˙x1 = x2 ... ˙xn−1 = xn= F (x(n−1)_{, ..., x, t)} 4.2 Implicit ODE

En implicit ODE är ett samband mellan flera beroende variabler och deras derivator, e.g. Newtons andra lag

F = md

2_r

dt2

där r = (x(t), y(t), z(t)) är ortsvektorn och där F kan vara en funktion av ortsvektorn och dess derivator. En n:te ordningens implict ODE skrives i allmänhet som

x(n)(t) = F (x(n−1), ..., x, t)

(x(n−1)_{(0), ..., x(0)) = (x}(n−1)

men kan mer generellt uttryckas på vektorform som F( ˙x, x, t) = 0

x(0) = x0

genom att utnyttja transformationen från en n:te ordningens implicit ODE till ett implicit system av n ODEs. Här är x = (x1(t), ..., xk(t))

tillståndsvek-torn och x0 begynnelsevärdet som entydigt bestämmer lösningen. På samma

sätt som en n:te ordningens explicit ODE kan skrivas om till ett system av

nODEs, så Låt oss nu avgränsa oss till system av implicita ODEs som kan

skrivas på formen

E(x) ˙x = f(x, t)

där E är en kvadratisk matris som antingen kan vara reguljär eller singulär och f är en vektorvärdfunktion. Om E är reguljär så existerar en invers till E och det implicita systemet kan då omvandlas till ett explicit system. Om E är singulär är detta inte möjligt, men ett specialfall uppstår om den är singulär på formen

E(x) = E1(¯x1) 0

0 0

!

, x= (¯xT₁,¯xT₂)T

där de nollskilda raderna bildar ett reguljärt block E1 och där vektorn x

skrivits om som två delar. Här ger de första nollskilda raderna ett implicit system av ODEs medan de sista nollraderna ger oss ickelinjära ekvationer, med andra ord vanliga algebraiska ekvationer. Följaktligen kallas ett sådanat blandat system för en differential algebraisk ekvation (DAE).

4.3 DAE

Ovan nämndes DAEen som ett specialfall av implicita ODEs. Vi ska nu dis-kutera vad för egenskaper som karaktärisrar en DAE och skiljer den från en ODE. Som ett löpande exempel skrivs DAEen på formen av ett mekaniskt system med holonoma restriktioner, där p, v och λ betecknar kolonnvektorer av generaliserade positioner och hastigheter samt Lagrange multiplikatorer, respektive (se [3] ss. 139-143 för en liknande genomgång). M är den regul-jära massmatrisen, g kolonnvektor med holonoma restriktioner och G är funktionalmatrisen av g

G ≡ dg

dp

˙p = v

M(p) ˙v = fa(p, v) − G(p)Tλ 0= g(p)

med begynnelsevärdena p(t0) = p0, v(t0) = v0 och λ(t0) = λ0. Konsistenta begynnelsevärden

Kravet på lösningen x till DAEen är att p ∈ C2_[t

0, te], v ∈ C1[t0, te], λ ∈

C0[t0, te], där Cn betecknar n kontinuerliga derivator. Detta kallas för ett

glatthetsvillkor. Begynnelsevärdet x(t₀) = x₀ måste därmed även uppfylla

bivillkoren 0 = g(p) liksom resten av DAEen. Sådana begynnelsevärden kallas då för konsistenta. Vid derivering av bivillkoret fås

0= G(p) ˙p = G(p)v 0= G(p) ˙v + _d dtG(p)  v ≡ G(p) ˙v + ξ(p, v)

En lösning exsisterar enbart om begynnelsevärdet x0 tillfredsställer dessa

två ekvationer. Med hjälp av resultaten fås att ˙v och λ kan bestämmas som lösningen till M(p) G(p)T G(p) 0 ! ˙v λ ! = _−ξfa(p, v)_{(p, v)} ! Vi får ˙v = M(p)−1 fa(p, v) − G(p)Tλ λ=G(p)M(p)−1G(p)T−1G(p)M(p)−1fa(p, v) + ξ(p, v) där

G(p)M(p)−1G(p)T är reguljär ty G antogs ha maximal rang. Kon-sistenta begynnelsevärden till DAEen kräver att x0 = (p0, v0, λ0), alltså

Index och DAEs med konstanta koefficienter

En linjär DAE med konstanta koefficienter kan skrivs på formen E˙x(t) = Ax(t) + f, x(t0) = x0

där E och A är kvadratiska matriser som kan vara singulära. I fallet för mekaniska system med holonoma tvång tar matriserna formen

E ≡    I 0 0 0 M 0 0 0 0   , A ≡    0 I 0 −K −D −GT G 0 0   

Lösningen för ekvationerna delas upp i tre fall som beaktas nedan i punkt-listan (se gärna [Eich-Soellner och Führer, 2002, ss. 56-60] för en liknande genomgång).

1. E är reguljär;

2. E är singulär och A är reguljär; 3. Både E och A är singulära.

Fall 1

I det här fallet är lösningen trivial. Eftersom E är reguljär kan vi vänster-multiplicera båda sidor med inversen E−1 _{och därmed omvandla systemet}

på explicit form

˙x(t) = E−1_{Ax(t) + E}−1_{f(t) = ˜Ax(t) + ˜f(t)} Fall 2

För detta fall gör vi följande omskrivningar av systemet A−1E˙x(t) = x(t) + A−1f(t)

där R korresponderar till de Jordan block med nollskilda egenvärden och N korresponderar till de med nollegenvärden. Här gäller att R är reguljär och N nilpotent R 0 0 N ! ˙¯x1 ˙¯x2 ! = ¯x1 ¯x2 ! + T−1_A−1_f | {z } ¯ f

För N gäller att det minsta talet k som uppfyller att Nk_{= 0 kallas} matri-sens index av nilpotens: k = ind(N). Indexet indikerar hur välställd DAEen är och hur många tidsderiveringar som fordras för att transformera DAEen till en ODE, vilket visas härnäst. Nu, systemet blir alltså

R ˙¯x1 = ¯x1+ ¯f1

N ˙¯x2 = ¯x2+ ¯f2

Eftersom R är reguljär, så är första ekvationen en explicit ODE, den så kal-lade tillståndsformen för systemet. Från den andra ekvationen får vi genom att tidsderivera k − 1 gånger

¯x2= − ¯f2+ N ˙¯x2 = − ¯f2+ N(− ˙¯f2+ N ¨¯x2) = − ¯f2+ N(− ˙¯f2+ N(−¨¯f2+ N ... ¯x2)) ... = − k−1 X i=0 Nif¯₂(i)

Detta är ett explicit uttryck för dellösningen ¯x2. Denna tillsammans med

första ekvationen leder till en ODE som kallas den underliggande ODEen (UODE) ˙¯x1= R−1( ¯f1+ ¯x1) ˙¯x2= − k−1 X i=0 Nif¯₂(i+1)

Fall 3

Om både E och A är singulära så är lösbarheten (att det existerar en lös-ning) beroende av om matrisknippet µE − A är reguljär eller inte. Det går att visa att om determinanten av matrisknippet är skilt från noll, alltså då

µE − Aär reguljär, så existerar en lösning. För ett utförligt bevis av detta

påstående se [3] ss. 59-64.

Indexreduktion för DAEs

Indexet för en DAE talar om hur välställt systemet är och hur många gånger den vektorvärda funktionen ¯f2 behöver tidsderiveras för att systemet skall

transformeras till ett system av ODEs. Vad ¯f2 representerar är inget annat

än de holnoma tvången (begränsningar på position) som vi sett innan ¯

f2 = g(p)

och är alltså algebraiska ekvationer. För att göra om en DAE till en ODE så tidsderiveras alltså begränsningarna på positionen. Som vi såg i delsektionen om konsistenta begynnelsevärden utfördes två tidsderiveringar varigenom en ODE för v erhölls. Genom att tidsderivera de positionella begränsningarna ytterligare en gång erhålls en ODE även för λ, varvid DAEen kan skrivas om som ett system av ODEs för alla de bereonde variablerna på formen

˙p = v ˙v = M(p)−1 fa(p, v) − G(p)T_λ ˙λ = d dt  G(p)M(p)−1G(p)T−1G(p)M(p)−1fa(p, v) + ξ(p, v) Detta visar att rörelseekvationerna för mekaniska system med holonoma tvång beskrivs av index-3 DAEs och metoden att gå från DAE till ODE kallas indexreduktion. I fallet ovan har indexet reducerats från 3 till 0. Det går även på analogt sätt visa att om enbart begränsningar på hastighet (eller acceleration) förekommer, så beskrivs dynamiken av en index-2 (index-1) DAE. För system med blandade begränsningar så gäller att indexet ges av det högst förekommande indexet, e.g. ett system med både positions- och hastighetsbegräsningar beskrivs av en index-3 DAE se [3] ss. 145-146.

5

Numeriska metoder

beräknas efter varje diskontinuitet, vilka kan betraktas som

fixpunktsitera-tioner i de fall systemet är attraherande. I det här avsnittet ges en kort

bakgrund om vad fixpunktsiteration är och hur konvergensen för en sådan iteration kan analyseras. Sist presenteras Picards algoritm, en numerisk me-tod för att finna fixpunkter.

5.1 Fixa punkter och fixpunktsiteration

Som ett första trivialt exempel av fixpunktsiteration och vad en fixpunkt är så betraktar vi cosinusfunktionen g(x) = cos x, se [4] ss. 30-40. Låt oss använda cosinusfunktionen för att bygga upp en sekvens av nummer genom följande regel, nämligen att utdatat blir det nya indatat. Om första data-punkten är given, säg x0, så är detta första talet i sekvensen och nästa ges

av g(x0) = cos x0 = x1. Talet efter det blir cos x1 = x2 och så vidare. Det

visar sig att den här sekvensen faktiskt konvergerar mot ett tal, nämligen

x∞ = 0.7390.... Den här punkten kallas för en attraherande fixpunkt till

g(x) = cos x och regeln vi tillämpade för att blottlägga den kallas

fixpunkt-siteration. Fixpunkten är attraherande i den mening att för vilket x0 vi än

väljer som första tal så konvergerar sekvensen till 0.7390... . Generellt gäller att ett reellt tal (eller vektor) r kallas fix om

g(r) = r

där g är en godtycklig funktion. En fixpunktsiteration beskrivs i allmänhet av den rekursiva formeln

x0= initial gissning

xi+1= g(xi), i = 0, 1, 2, ...

5.2 Kontraktiva sekvenser och Picards algoritm

Konvergensen för en iterationssekvens kan uppskattas med hjälp av det

re-lativaoch absoluta felet. Här förutsätts att fixpunkten från början känns till,

men så är sällan fallet. Det vi däremot kan beräkna är ifall sekvensen av erhållna värden från en fixpunktsiteration är kontraktiv eller inte. Begreppet kontraktiv definerias matematiskt i följande stycke:

g: [a, b] → R kallas kontraktiv på [a, b] iff för alla x och y i [a, b] och r < 1

|g(x) − g(y)| ≤ r|x − y|

sekvenser, e.g. g(x) = x2_{, x >0 med fixpunkten x = 1, således fungerar}

in-te vanlig fixpunktsiin-teration. (Däremot kan fixpunkin-ten beräknas analytiskt). Algoritmen är relativt simpel och lyder

Data: g= g(x)

Result: Fixpunkt r (g(r) = r)

1 initialisering;

2 gör initialgissning x₀ i [a, b] och sätt  > 0; 3 skapa x1 = g(x0); 4 while |x₁− x₀| >  do 5 if |g(x₁) − x₁| >  then 6 avsluta loop; 7 else 8 x₀← x₁; 9 x1← g(x0); 10 end 11 end

Algorithm 1:Picards algoritm

Picards algoritm tillämpad på ett diskontinuerligt system

För att beräkna perioden av ett diskontinuerligt system är det inte självklart hur den ifrågavarande algoritmen skall lösa problemet, vilken är avsedd för kontinuerliga intervall. För ett periodiskt stabilt och attraherande system bestående av flera var för sig kontinuerliga rörelsefaser vet vi att om begyn-nelsevärdet för någon rörelsefas är detsamma som efter en period/cykel, då har systemet konvergerat fullständigt till den stabila omploppsbanan. För ett attraherande system sker fullständig konergens i limes av oändligheten för den beroende variabeln, ofta tiden t. Därav är det praktiskt att an-vända ett intervall runt den stabila banan inom vilket lösningen sägs ha konvergerat, ett så kallat toleransintervall. Efter systemet konvergerat är det en enkel uppgift att beräkna perioden (ta tiden på nästa period/cykel bara!). Antagandent om att omloppsbanan är attraherande implicerar att fixpunktsiterationen är konvergent och måste bevisas matematiskt. (Fastän omloppsbanan tycks vara attraherande för några godtyckligt valda begyn-nelsevärden, vilket visas i Figurerna 7 och 11 i Sektion 9, så visar inte detta att det gäller i allmänhet!) Ingen tid kommer att ägnas åt att bevisa det här, utan läsaren hänvisas till [1] för en redogörelse av det.

nästkommande cykel tar att slutföras, se Sektion 7.4 och 8.4 för implemen-tationen i Matlab.

5.3 Numerisk lösning av DAEs

Matlab erbjuder inbyggda numeriska lösare som kan lösa DAEs av in-dex 1 eller lägre, bland annat ode15s som skall användas i simulationen av hackspettsleksaken i denna studie. Eftersom mekaniska system vanligtvis beskrivs av DAEs med index högre än 1 så är det nödvändigt att omformu-lera DAEen först genom indexreduktion vilket diskuterades i Sektion 4.3. Det är möjligt att lösa DAEs med index högre än 1, men dessa ger upphov till ytterligare problem såsom hur mängden konsistenta begynnelsevärden skall bestämmas och metoder för feluppskattningen, se [12], vilket gör om-formuleringsmetoden mer effektiv. Varför konsistenta begynnelsevärden är ett problem beror av gömda algebraiska begränsningar i högre index system. Låt oss betrakta en DAE som ges av

x1= q(t)

x2= ˙x1

För att göra om systemet till explicita ODEs krävs att första och andra ekvationen deriveras en gång: ˙x1 = ˙q(t) och ˙x2 = ¨x1 = ¨q(t). DAEen har

alltså indexet 2. Vad finns för begränsningar på begynnelsevärdena? Uppen-barligen måste x1(0) = q(0) enligt första ekvationen, men det finns även en

gömd algebraisk begränsning nämligen x2(0) = ˙q(0). Detta är en viktig

skill-nad mellan index 1 och högre index DAEs vilket ger upphov till uppenbara svårigheter för numeriska lösare.

6

Diskontinuerliga system och växlingsalgoritmer

avstånden de emellan blir noll! Låt säga att funktionen som ger avståndet mellan kropparn är

gN = gN(p)

där N betecknar normalavståndet, då kallas gN för en växlingsfunktion eller

händelsefunktion. När gN byter tecken (från + till − i detta fall), så vet vi att

Modeller och algoritmer

I denna del formuleras två modeller för hackspettsleksakens rörelse. I första fallet modelleras systemet med två frihetsgrader och i andra fallet av tre fri-hetsgrader. Holonoma begränsningar förekommer i båda fallen och i andra fallet tas även hänsyn till kontaktkrafter. Den första modellen kan härledas med Lagranges ekvationer av typ två vilket visas. Respektive modell löses numeriskt i Matlab med hjälp av en växlingsalgoritm, sedan beräknas re-spektive modells period med hjälp av Picards algoritm som presenterades i Sektion 5.

Figur 3: Hackspettsleksaken. Den här bilden får användas med tillåtelse av Remco I. Leine, Nonlinear Dynamics Of The Woodpecker Toy, 2001.

7

Simpel modell

7.1 Systembeskrivning och härledning

Figur 4: Hackspettsleksaken. Den här bilden får användas med tillåtel-se av Edda Eich-Soellner, Projizierende Mehrschrittverfahren zur

numeri-schen Lösung von Bewegungsgleichungen technisher Mehrkörpersysteme mit Zwangsbedingungen und Unstetigkeiten Dissertation, 1991.

Hackspettens rörelse kan sammanfattas av fem faser: 1. θ > θK1:

Hackspetten svänger medurs tills den når sin maximala amplitud och därefter tillbaka i motursriktning. Den faller inte nedåt.

2. −θK1< θ < θK1:

Fas 1 avslutas då θ = θK1och hackspetten faller nedåt som konsekvens av tyngdkraften. Hackspetten roterar moturs ( ˙θ < 0).

3. −θK2< θ < −θK1:

För θ = −θK1så sker en krock mellan ring och stav varvid hackspetten slutar falla, men dess rotation är fortfarande moturs.

4. −θK2< θ < −θK1:

För θ = −θK2 så krockar hackspetten med sin näbb i staven och rotationen byter riktning från moturs till medurs. Hackspetten rör sig fortfarande inte nedåt, utan är fast.

5. −θK1< θ < θK1:

Det förekommer alltså två typer av rörelse, antingen befinner sig hackspetten och ringen i fritt fall samtidigt som båda kroppar roterar, eller så är ringen fast i staven medan enbart hackspetten roterar. I det första fallet beskrivs rörelsen av två frihetsgrader nämligen hackspettens vinkelutslag θ och läge vinkelrät staven z, och i det andra fallet ges rörelsen bara av frihetsgraden

θeftersom z antas vara konstant när ringen är låst till staven.

Koordinatsy-stemet visas i Figur 4, där positiva riktningen för z-axeln är uppåt (notera att pilen bredvid z i Figur 4 som pekar nedåt visar hastighetsriktningen!). Nu, utifrån dessa antaganden kan rörelseekvationerna härledas vilket görs enklast med Lagrange ekvationer av andra typen. Notera att denna metod bara kan användas då begränsningarna är holonoma, vilket är uppfyllt för det här systemet. I övergångsfasen mellan dessa två stadier sker en

hän-delse som sätter villkor på systemets tillstånd efteråt. Händelserna sker för

vinklarna ±θK1 och θK2

7.1.1 Fritt fall

När hackspetten faller fritt utför masscentrum en vertikal nedåtgående rö-relse till följd av tyngdkraften. Då systemet antas vara friktionslöst påver-kar inte ringens och hackspettens rotation masscentrum i vertikalled, ty det finns ingen kraft som kan ge upphov till sådan rörelse. Däremot utför masscentrum en viss rörelse i horisontalled, men denna kan med god nog-grannhet försummas för små vinklar (vid stora vinklar fås däremot ett för lågt masströghetsmoment i modellen). Vi utgår från systemets masscent-rum i G och sätter höjdkoordinaten där som vi betecknar med ξ istället för z vilket är hackspettens höjdkoordinat. Detta underlättar härledningen eftersom translationshastigheten nu ges av en term som vi skall se. Målet är nu att formulera energiuttrycken för systemet, alltså uttrycken för kinetisk respektive potentiell energi. Om ringens vridmoment försummas erhålls

T = 1 2m ˙ξ2 | {z } (1) +1₂I2˙θ2 | {z } (2) +1₂m2(r2˙θ)2 | {z } (3) +1₂m1(r1˙θ)2 | {z } (4)

där T representerar den kinetiska energin. Termerna är numrerade i den ord-ning de kommer. Den första termen beskriver energibidraget från masscent-rumets translationshastighet nedåt, där m = m1+m2 och m1 respektive m2

är ringens respektive hackspettens massa. Andra termen ges av hackspettens rotationsrörelse kring sig själv som antas vara detsamma som rotationsrö-relsen runt G (denna rörelse kallas för låst rotation) där I2 är hackspettens

masströghetsmoment. Tredje respektive fjärde termen uppkommer från att hacskpetten respektive ringen utför en centralrörelse runt G där r1 och r2

r1 = b m2 m1+ m2 r2 = b m1 m1+ m2

där b är avståndet mellan ringens och hackspettens masscentrum. Notera att ringen inte ger något energibidrag i form av rotationsrörelse runt sig själv då vi antog att den inte utför något vridande moment. För den potentiella energin erhålls (när ξ-axeln är riktad nedåt!)

V = −mgξ | {z } (1) +1₂cθ2 | {z } (2)

där V representerar systemets potenteilla energi. Första termen är naturligt-vis den energi som är lagrad i systemet på grund av tyngdkraften; här sätts nollnivån där systemet börjar falla ifrån. Eftersom den positiva riktningen är nedåt, så fordras ett minustecken framför termen, ty annars skulle syste-mets totala energi öka. Den andra termen är den lagrade energin i fjädern där c är fjäderkonstanten.

Vi har nu motiverat och formulerat energiutrycken, nästa steg är att ställa upp Lagrangefunktionen och använda Lagranges ekvationer för att slutligen beräkna systemets rörelseekvationer. Från Sektion 2.6.4 fås

L= T − V och d dt _∂L ∂˙pk  − ∂L ∂pk = 0

där L är Lagrangefunktionen och p de generaliserade koordinaterna, det vill säga ξ och θ i vårt fall. Med uttrycken för T och V insatta erhålls

L= 1 2m˙z2+ 1 2I2˙θ2+1 2m2(r2˙θ)2+1 2m1(r1˙θ)2+ mgξ −1 2cθ2 Lagranges ekvationer ger då

Sista steget är att höjdkoordinaten ska utgå från hackspetten istället för

G, alltså måste ξ uttryckas med z. Genom att ställa upp geometrin i

fi-guren och använda att för små vinklar gäller sin θ ≈ tan θ ≈ θ så erhålls transformationen

ξ= r2θ+ z ⇒ ¨ξ = r2¨θ+ ¨z

vilket insatt i ekvationerna ovan slutligen ger, efter lite omskrivning  I2+ m2b2  1 − m2 m1+ m2  = −cθ (m1+ m2)¨z + m1b¨θ= (m1+ m2)g 7.1.2 Fast i staven

För vissa vinklar fastnar ringen i staven och translationsrörelsen nedåt orsa-kad av tyngdkraften upphör momentant. Dessa vinklar betecknas ±θK1och utanför intervallet mellan dem är systemet fast. Då ändras systembeskrv-ningen eftersom en frihetsgrad förloras, eller är överflödig, nämligen höjd-koordinaten z som förblir konstant under rörelsen. Den observante läsaren skulle invända att hackspetten faktikst rör sig i höjdled när den roterar och det är sant, det är inget modellen bortser ifrån, men höjden kan uttryckas med θ. På samma sätt som innan så härleder vi först energiuttrycken för att sedan kunna formulera Lagrangefunktionen och med Lagranges ekvationer erhålla rörelseekvationerna för systemet. Den kinetiska energin skrives

T = 1 2I2˙θ2 | {z } (1) +1₂m2(b ˙θ)2 | {z } (2)

där termerna (1) och (2) motiveras på samma sätt som termerna (2) re-spektive (3) för uttrycket av den förra kinetiska energin. Notera att r2 → b

i (2) med anledning av att hackspetten utför en centralrörelse runt kan-ten av ringen nu istället för runt det gemensamma masscentrumet G. Den potentiella energin blir

V = 1 2cθ2 | {z } (1) − m₂bθg | {z } (2)

där första termen (1) spelar samma roll (2) i det förra uttrycket för den potentiella energin. Andra termen i uttrycket ovan kommer utav att hack-spetten återigen lagrar energi på grund av den tyngdkraft som påverkar den. Skillnaden från förut är att systemets totala massa enbart ges av hack-spettens massa m2 vars masscentrum är beläget avståndet b från fjäderns

Figur 5: Geometri Lagrangefunktionen blir således

L= 1

2I2˙θ2+1

2m2(b ˙θ)2−1

2cθ2+ m2bθg

och med hjälp av Lagranges ekvationer fås slutligen

d dt _∂L ∂ ˙θ  = ∂L ∂θ ⇒ (I2+ m2b 2_{)¨θ = −cθ + m} 2bg

Notera att ingen koordinattransformation krävdes i det här fallet eftersom systemets masscentrum redan sammanföll med höjdkoordinaten för hack-spetten.

7.1.3 Händelser och fasövergångar

För att gå från det ena stadiet till det andra måste en händelse inträffa. Det finns fem händelser varav tre sätter speciella villkor på tillståndet efteråt eftersom en stöt inträffar vilket ger en impuls i vinkelhastigheten. Dessa impulser har redan beräknats i tidigare arbeten [Pfeiffer, 1984] och vinkel-hastigheten före och efter för dessa tre händelser ges av

• ringens överkant stöter i staven: ˙θ+_{= (1 − d} 3)  ˙θ−₊ m2b I2+m2b2 ˙z −_{≡ ˙θ}+ 2

• ringens underkant stöter i staven: ˙θ+_{= (1 − d} 1)  ˙θ−₊ m2b I2+m2b2 ˙z −_{≡ ˙θ}+ 4 m1 = 0.003 kg I2 = 7 · 10−7 kgm2 d1= 0.18335 θK2= 12◦ m2 = 0.0045 kg c= 0.0056 Nm θK1= 10◦ b= 0.015 m g= 9.81 m/s2 d3 = 0.04766

• hackspettens näbb stöter i staven: ˙θ+_{= − ˙θ}−_{≡ ˙θ}+

3

där ˙θ−_{, ˙z}− _{är hastigheter innan stöt och ˙θ}+_{, ˙z}+ _{hastigheterna efter stöt.}

Konstanterna d1, d3 är studstalen för stötarna vilka beräknats

experimen-tellt. För de resterande händelserna sker en jämn fasövergång, det vill säga inget hopp i vinkelhsatigheten. För sakens skull listas alla händelser i den ordning de följer varandra, vilka fasövergångar det leder till och hur vinkel-hastigheten ändras.

• ( ˙θ > 0)

1. θK1 : Tillstånd II → Tillstånd I (jämn övergång: ˙θ+_{= ˙θ}−₎

2. −θK1: Tillstånd I → Tillstånd II ( ˙θ+_{= ˙θ}+ 2)

3. −θK2: Tillstånd II → Tillstånd II ( ˙θ+_{= ˙θ}+ 3)

• ( ˙θ < 0)

4. −θK1: Tillstånd II → Tillstånd I (jämn övergång) 5. θK1 : Tillstånd I → Tillstånd II ( ˙θ+_{= ˙θ}+

4)

Värdena för samtliga parametrar i modellen visas i Tabell 1. 7.2 Algoritm

Data: ˙x = f(x), x(0) = x0 Result: x= x(t)

1 while antalet inträffade händelser n < N do 2 beteckna den nyss inträffade händelsen med e0;

3 if e₀ = 1 then

4 lös Tillstånd I med B.V. x₀ tills nästa händelse e₁ inträffar; 5 beräkna tillståndet x1 efter händelsen;

6 end

7 if e0 = 2 then

8 lös Tillstånd II med B.V. x₀ tills nästa händelse e₁ inträffar; 9 beräkna tillståndet x₁ efter händelsen;

10 end

11 if e0 = 3 then

12 lös Tillstånd II med B.V. x0 tills nästa händelse e1 inträffar;

13 beräkna tillståndet x₁ efter händelsen; 14 end

15 if e₀ = 4 then

16 lös Tillstånd I med B.V. x0 tills nästa händelse e1 inträffar;

17 beräkna tillståndet x₁ efter händelsen; 18 end

19 if e₀ = 5 then

20 lös Tillstånd II med B.V. x₀ tills nästa händelse e₁ inträffar; 21 beräkna tillståndet x1 efter händelsen;

22 end 23 e0 ← e1;

24 n ← n+ 1; 25 x0← x1;

26 end

Algorithm 2:Algoritm för hackspettsleksaken

7.3 Implementation i Matlab

Varje steg i implementeringen visas och motiveras. För att lösa systemet så skriver vi om alla ekvationer på explicit form, således kan den inbyggda Matlab lösaren ode45 tillämpas. Första steget är att definiera konstanter, initialvärden, integrationsintervall och avsätta listor för utdatat för varje variabel. Notera att integrationsintervallet skall göras tillräckligt långt för att en händelse skall inträffa under integrationen.

% General constants

I2 = 7*10^( -7) ; % [ kgm ^2]

c = 0.0056; % [Nm] d1 = 0.18335; d3 = 0.04766; g = 9.81; % [m/s ^2] K1 = 0.174532925; % [ rads ] K2 = 0.20943951; % [ rads ] p = [I2 m1 m2 b c d1 d3 g K1 K2 ]; % parametervector % Constants to the function firstODE

A1 = I2 + m2*b ^2; B1 = m2*b*g;

% Constants to the function secondODE

A2 = I2 + m2*b ^2*(1 - m2 /( m1+m2)); B2 = m1 + m2; C2 = m2*b; D2 = (m1 + m2)*g; IC = [0.1745 10.2547 0 0]; % Initial Condition tstart = 0; tfinal = 10; y0 = IC; tout = tstart ;

% Lists for different statevariables output

thout = y0 (1:1) ; thvout = y0 (2:2) ; zout = y0 (3:3) ; zvout = y0 (4:4) ; ieout = 1; yeout = [];

Vi vill utnyttja att ode45 kan hantera händelser och använder odeset för att ställa in lösaren på evaluera händelsefunktionen för att upptäcka händelser under integrationen.

options = odeset ('Events ', @myEventsFcn );

Nästa steg är att koda huvudprogrammet. Det är denna del algortimen beskriver (fast i pseudokod).

% Main program

N = 10

for i = 1:N % Change this to 5 or 4 for one cycle

event = ieout (end); switch event

case 1

[t,y,te ,ye ,ie] = ode45(@(t,y) secondODE (t,y,A2 ,B2 ,C2 ,D2 ,←-c) ,[0 tfinal ],y0 , options );

[thout , thvout ,zout ,zvout , tstart , tout ]= Update (y,t,thout ,←-thvout ,zout ,zvout , tstart ,tout ,2) ;

y0 = ICgenerate (ye ,ie ,p); case 2

y00 = [y0 (1) y0 (2) ];

[t,y,te ,ye ,ie] = ode45(@(t,y) firstODE (t,y,A1 ,B1 ,c) ,[0 ←-tfinal ],y00 , options );

[thout , thvout ,zout ,zvout , tstart , tout ]= Update (y,t,thout ,←-thvout ,zout ,zvout , tstart ,tout ,1) ;

ie = ie(end); ieout = [ ieout ie ];

ye = [ye(end,:) zout (end) zvout (end)]; y0 = ICgenerate (ye ,ie ,p);

case 3

y00 = [y0 (1) y0 (2) ];

[t,y,te ,ye ,ie] = ode45(@(t,y) firstODE (t,y,A1 ,B1 ,c) ,[0 ←-tfinal ],y00 , options );

[thout , thvout ,zout ,zvout , tstart , tout ]= Update (y,t,thout ,←-thvout ,zout ,zvout , tstart ,tout ,1) ;

ie = ie(end); ieout = [ ieout ie ];

ye = [ye(end,:) zout (end) zvout (end)]; y0 = ICgenerate (ye ,ie ,p);

case 4

[t,y,te ,ye ,ie] = ode45(@(t,y) secondODE (t,y,A2 ,B2 ,C2 ,D2 ,←-c) ,[0 tfinal ],y0 , options );

[thout , thvout ,zout ,zvout , tstart , tout ]= Update (y,t,thout ,←-thvout ,zout ,zvout , tstart ,tout ,2) ;

ie = ie(end); ieout = [ ieout ie ]; ye = ye(end,:) ;

y0 = ICgenerate (ye ,ie ,p); case 5

y00 = [y0 (1) y0 (2) ];

[t,y,te ,ye ,ie] = ode45(@(t,y) firstODE (t,y,A1 ,B1 ,c) ,[0 ←-tfinal ],y00 , options );

[thout , thvout ,zout ,zvout , tstart , tout ]= Update (y,t,thout ,←-thvout ,zout ,zvout , tstart ,tout ,1) ;

ie = ie(end); ieout = [ ieout ie ];

ye = [ye(end,:) zout (end) zvout (end)]; y0 = ICgenerate (ye ,ie ,p);

end end

Härefter definieras alla funktioner: Händelsefunktion, ODE-funktioner, ini-tialvärdesfunktion och listornas uppdateringsfunktion.

function [thout , thvout ,zout ,zvout , tstart , tout ] = Update (y,t,thout ,←-thvout ,zout ,zvout , tstart ,tout ,i)

% This function updates the variables in a correct way since the % degrees of freedom change , thus some variables are constant

or←-zero .

nt = length(t);

if i == 1 % first ODE

zstat = repmat ( zout (end) ,1,nt -1) ; zout = [ zout zstat ];

zvstat = zeros(nt -1 ,1); zvout = [ zvout zvstat '];

elseif i == 2 % second ODE

zout = [ zout y (2: nt ,3) ']; zvout = [ zvout y (2: nt ,4) '];

end end

Funktionen ovan uppdaterar listorna. Beroende av vilken fas systemet är i, så uppdateras variablernas listor annorlunda. Detta kommer utav att antalet frihetsgrader varierar.

function dy = firstODE (t,y,A1 ,B1 ,c) dy = zeros(2 ,1);

dy (1) = y (2) ;

dy (2) = 1/( A1)*(B1 -c*y (1) );

end

function dy = secondODE (t,y,A2 ,B2 ,C2 ,D2 ,c) dy = zeros(4 ,1); dy (1) = y (2) ; dy (2) = -c/A2*y (1) ; dy (3) = y (4) ; dy (4) = D2/B2 + c*C2 /( A2*B2)*y (1) ; end

Ovan definieras systemen av ODEs som vi är intresserade att lösa. Nästa funktion i ledet är händelsefunktionen som definieras nedan.

function [value , isterminal , direction ] = myEventsFcn (t,y) K1 = 0.174532925; % 0.174532925 rads K2 = 0.20943951; % 0.20943951 rads value = [y (1) -K1 , y (1) +K1 , y (1) +K2 , y (1) +K1 , y (1) -K1 ]; isterminal = [1, 1, 1, 1, 1]; direction = [-1, -1, -1, 1, 1]; end

Slutligen har vi initialvärdesfunktionen som generar nya begynnelsevärden inför varje integrationssteg.

function y0 = ICgenerate (ye ,ie ,p)

% Generates new Initial Conditions based on which event happend

case 1

tpos = thetapos ; tvel = thetavel ; zpos = zetapos ; zvel = 0;

y0 = [ tpos tvel zpos zvel ]; case 2 tpos = thetapos ; tvel = (1-p (7) )*( thetavel +p (3) *p (4) /(p (1) +p (3) *p (4) ^2) *←-zetavel ); zpos = zetapos ; zvel = 0;

y0 = [ tpos tvel zpos zvel ]; case 3

tpos = thetapos ; tvel = -thetavel ; zpos = zetapos ; zvel = 0;

y0 = [ tpos tvel zpos zvel ]; case 4

tpos = thetapos ; tvel = thetavel ; zpos = zetapos ; zvel = 0;

y0 = [ tpos tvel zpos zvel ]; case 5 tpos = thetapos ; tvel = (1-p (6) )*( thetavel +p (3) *p (4) /(p (1) +p (3) *p (4) ^2) *←-zetavel ); zpos = zetapos ; zvel = 0;

y0 = [ tpos tvel zpos zvel ];

end end

7.4 Beräkning av perioden

periodlista = [];

Nästa steg är att skapa en if-sats i huvudprogrammet under någon av de olika case-satserna, låt säga case = 5 (då hackspetten roterar moturs i fritt fall). if-satsen kräver att simulationen körts minst en fas, därav i > 0; på så sätt registreras starttiden för cykeln exakt.

if i > 1

periodlista = [ periodlista tout (end)];

end

Perioderna för varje cykel fås då med hjälp av en inbyggd Matlab funktion som tar skillnanden mellan grannkomponenterna i periodlistan vari start-och sluttiderna för cyklerna finns.

periods = diff( periodlista );

Perioderna kan sedan plottas mot cyklerna. På det här sättet ser vi även hur snabbt sekvensen konvergerar.

8

Komplex modell

På samma sätt som asvnittet innan presenteras här den komplexa model-len av hackspettsleksaken, en mer sofistikerad modell som tar hänsyn till

kontaktkrafter för att avgöra när en fasövergång sker. (Kontaktkrafter är

8.1 Systembeskrivning

Figur 6: Hackspettsleksaken. Den här bilden får användas med tillåtelse av Fatemeh Mohammadi, Starters for Multistep Methods in the Solution of

Discontinous ODEs, 2016.

Modellen som skall presenteras nu är hämtad från [Mohammadi, 2016]. För-utom att den här modellen även tar hänsyn till kontaktkrafter till skillnad från den simpla modellen, så har systemet även fler frihetsgrader. Ring-en har två frihetsgrader, dess rotation φS och dess vertikala translation

z. Hackspetten har bara en frihetsgrad nämligen dess rotation φB. Dessa

tre frihetsgrader utgör systemets generaliserade koordinater (den minimala koordinatframställningen). Hackspettens rörelse kan sammanfattas på ett liknande sätt som för den simpla modellen med hälp av ett cykelschema. Skillnaden från innan är att fasen systemet befinner sig i även är beroende utav ringens rotation φS och kontaktkrafterna λ1 och λ2. Detta beroende

betonas i cykelschemat nedan.

Kontaktkrafterna är endast aktiva när systemet inte befinner sig i fritt fall (λ1= λ2 = 0 då) och ges då av λ1= cp(φB− φS) − rSλ2− mBlSlGφ¨B− mBlSg hS λ2= −(mS+ mB)g − mBlGφ¨B 1. Start vid φS = −rS−r0 hS , ˙φB<0:

Syste-2. Start då λ1= 0, ˙φB>0:

Ringen lösgör sig från staven och som konsekvens faller systemet nedåt. Hackspetten rör sig moturs och inducerar samma rotation hos ringen. Fas 2 avslutas då φS = rS−r0

hS .

3. Start vid φS = rS−r0

hS , ˙φB>0:

Ringens övre högra och vänstra nedre hörn låser sig mot staven vilket ger en impuls i hackspettens rotationshastighet. Systemet slutar falla nedåt (kontaktkraften aktiveras: λ1 6= 0). Fas 3 avslutas då φB =

lS+lG−l−B−r0

hB .

4. Start φB = lS+lG−l−B−r0

hB , ˙φB>0:

Hackspetten stöter i staven med näbben i en motursrörelse för att därefter momentant vända riktning. Systemet är fortfarande fast. Fas 4 avslutas då λ1 = 0.

5. Start då λ1= 0, ˙φB<0:

Ringen lösgör sig och systemet börjar återigen falla. Hackspetten ro-terar medurs och inducerar samma rotation hos ringen såsom innan. Fas 5 avslutas då φS = −rS−r0

hS varvid fas 1 börjar om.

I följande delavsnitt ges rörelseekvationerna för varje fas som matematiskt ges av system av DAEs. Vidare måste dessa system av DAEs indexreduceras eftersom de numeriska lösarna som tillämpas här fordrar index ≤ 1. Vi kommer upptäcka att vissa faser beskrivs av samma rörelseekvationer, vilket påbjuder oss att beskriva systemets dynamik i stadier snarare än faser. Nedan följer en kort förklaring till varje stadie.

8.1.1 Stadie I - fritt fall

Systemet (ringen och hackspetten) faller fritt nedför staven, ingen blocke-ring sker. Stadie I beskrivs av

(mS+ mB)¨z + mBlSφ¨S+ mBlGφ¨B = −(mS+ mB)g

(mBlS)¨z + (JS+ mBl2S) ¨φS+ mBlSlG = cp(φB− φS) − mBlSg

mBlG¨z + mBlSlGφ¨S+ (JB+ mBl2G) = cp( ¨φS− ¨φB) − mBlGg 0 = λ1

0 = λ2

8.1.2 Stadie II - ring fast i nedre läge

Systemet är fast, där ringens undre högra och övre vänstra hörn är inspända i staven; kontaktkraftern får nollskilda värden. Stadie II beskrivs av

(mS+ mB)¨z + mBlSφ¨S+ mBlGφ¨B= −(mS+ mB)g − λ2

(mBlS)¨z + (JS+ mBl2_S) ¨φS+ mBlSlGφ¨B= cp(φB− φS) − mBlSg − hSλ1− rSλ2

mBlG¨z + mBlSlGφ¨S+ (JB+ mBl2G) ¨φB= cp(φS− φB) − mBlGg 0 = (rS− r₀) + hSφS 0 = ˙z + rSφ˙S

Den sista ekvationen är en icke-holnom begränsning och innehåller tidsde-rivator av de generaliserade koordinaterna. Detta indikerar att DAEn har index = 2. Genom att derivera den nästsista respektive sista ekvationen två gånger vardera erhåller vi

¨z = 0 ¨

φS = 0 Detta transfomerar DAEn till formen

¨z = 0 ¨

φS = 0

(JB+ mBl_G2) ¨φB = cp(φS− φB) − mBlGg

8.1.3 Stadie III - ring fast i övre läge

På samma sätt som innan utförs en indexreduktion i syfte att göra systemet numeriskt lösbart. DANn transfomeras till formen

¨z = 0 ¨

φS = 0

(JB+ mBlG) ¨2 φB = cp(φS− φB) − mBlGg

vilket är exakt samma form som för andra stadiet. Vi drar slutsatsen att andra och tredje stadiet beskrivs av samma rörelseekvationer, men priset vi betalat är att information om ringens (konstanta) höjd och rotation förlo-rats. Detta innebär däremot inte ett problem vid simuleringen som vi skall se.

8.1.4 Händelser och fasövergångar

Händelserna och fasövergångarna har redogjorts för i cykelschemat ovan. Målet nu är att bestämma impulsen för hackspettens vinkelhastighet vid samma händelser som för den simpla modellen. Skillnaden mellan model-lerna är att impulserna beräknas genom antagandet om att rörelsemängden bevaras

I−= I+

där I− _{och I}+ _{representerar rörelsemängden efter respektive innan.}

Rörel-semängden innan ges av

I−= mBlG˙z−+ (mBlSlG) ˙φS

−

+ (JB+ mBl2G) ˙φB

−

och rörelsemängden efter ges av

I+= mBlG˙z++ (mBlSlG) ˙φS

+

+ (JB+ mBl2G) ˙φB

+

Eftersom ringens rotations- och translationshastighet är noll efter stöt, alltså ˙

φS = ˙z = 0, så fås hackspettens vinkelhastighet efter stöt som

I+= I−⇒ ˙φB + = mBlG˙z−+ (mBlSlG) ˙φS − + (JB+ mBl2_G) ˙φB − (JB+ mBl2 G) Här är ˙z−_{, ˙φ} S − och ˙φB −

hastigheterna innan stöt. Stöt förekommer i tre fall: Ringen stöter i staven i uppåt och nedåt riktning samt hackspettens näbb stöter i staven. När ringen stöter i staven gäller alltså ovannämnda formel i båda fall och när hackspetten söter i staven gäller ˙φB

+

= − ˙φB

−

mS = 3 · 10−4 kg JB = 7 · 10−7 kgm2 lG= 1.5 · 10−2 m rS = 3.1 · 10−3 m

mB = 0.0045 kg cp= 5.6 · 10−3 N/rad lS = 10−2 m r0 = 2.5 · 10−3 m

g= 9.81 m/s2 hS = 2 · 10−2 m JS = 5 · 10−9 kgm2 hB= 2 · 10−2

lB= 2.01 · 10−2 m

Tabell 2: Parametrar som används för den komplexa modellen.

8.2 Algoritm

Algortimen har principiellt samma struktur som den för den simpla model-len. Skillnaden blir här att x(t) = (z, φS, φB,˙z, ˙φS, ˙φB).

Data: ˙x = f(x), x(0) = x0 Result: x= x(t)

1 while antalet inträffade händelser n < N do 2 beteckna den nyss inträffade händelsen med e₀; 3 if e₀ = 1 then

4 lös TillståndI med B.V. x0 tills nästa händelse e1 inträffar;

5 beräkna tillståndet x₁ efter händelsen; 6 end

7 if e₀ = 2 then

8 lös TillståndII med B.V. x₀ tills nästa händelse e₁ inträffar; 9 beräkna tillståndet x1 efter händelsen;

10 end

11 if e0 = 3 then

12 lös TillståndII med B.V. x₀ tills nästa händelse e₁ inträffar; 13 beräkna tillståndet x₁ efter händelsen;

14 end

15 if e0 = 4 then

16 lös TillståndI med B.V. x₀ tills nästa händelse e₁ inträffar; 17 beräkna tillståndet x₁ efter händelsen;

18 end

19 if e₀ = 5 then

20 lös TillståndII med B.V. x0 tills nästa händelse e1 inträffar;

21 beräkna tillståndet x₁ efter händelsen; 22 end

23 e₀ ← e₁; 24 n ← n+ 1; 25 x0← x1;

26 end

8.3 Implementation i Matlab

% A l l m n n a konstanter mS = 3*10^( -4) ; % [kg] JS = 5*10^( -9) ; % [ kgm ] mB = 4.5*10^( -3) ; % [kg] g = 9.81; % [m/s ^2] JB = 7*10^( -7) ; % [ kgm ] rS = 3.1*10^( -3) ; % [m] lS = 10^( -2) ; % [m] r0 = 2.5*10^( -3) ; % [m] hS = 2*10^( -2) ; % [m] lG = 1.5*10^( -2) ; % [m] lB = 2.01*10^( -2) ; % [m] hB = 2*10^( -2) ; % [m] cp = 5.6*10^( -3) ; % [N/ rad ] % Massmatriser

% y = (z ,\ phiS ,\ phiB ,\xi ,\ny ,\ Lambda ) % \ dot {z }=\ xi

% \ dot {\ phiS }=\ ny % \ dot {\ phiB }=\ Lambda

MI = [1 ,0 ,0 ,0 ,0 ,0; 0 ,1 ,0 ,0 ,0 ,0 0 ,0 ,1 ,0 ,0 ,0; 0 ,0 ,0 ,( mS+mB),mB*lS ,mB*lG; 0 ,0 ,0 ,( mB*lS) ,(JS+mB*lS ^2) ,mB*lS*lG; 0,0,0, mB*lG ,mB*lS*lG ,( JB+mB*lG ^2) ]; MII = [1 ,0 ,0 ,0 ,0 ,0; 0 ,1 ,0 ,0 ,0 ,0 0 ,0 ,1 ,0 ,0 ,0; 0 ,0 ,0 ,1 ,0 ,0; 0 ,0 ,0 ,0 ,1 ,0; 0 ,0 ,0 ,0 ,0 ,( JB+mB*lG ^2) ];

% ODEset for the two different states % Initial conditions y0 = [0 0 -0.35 0 0 0]; % I n i t i a l v r d e tstart = 0; tfinal = 2; tspan = [0 tfinal ]; ieout = 2;

% Preallocated lists for each solution to each variable

tout = tstart ; % Lista avsatt f r tiden

zout = y0 (1:1) ; % Lista avsatt f r z- koordinater f r hackspetten B

sout = y0 (2:2) ; % Lista avsatt f r vinkeln f r ringen S

bout = y0 (3:3) ; % Lista avsatt f r vinkeln f r hackspetten B

zvout = y0 (4:4) ; % Lista avsatt f r z- hastigheter f r hackspetten ←-B

svout = y0 (5:5) ; % Lista avsatt f r vinkelhastigheten f r ringen S

bvout = y0 (6:6) ; % Lista avsatt f r vinkelhastigheten f r ←-hackspetten B

garante-respektive massmatris och händelsfunktion in.

optionsI = odeset ('MaxStep ',10^( -5) ,'Mass ',MI ,'Events ', @myEventsFcnI←-);

optionsII = odeset ('MaxStep ',10^( -5) ,'Mass ',MII ,'Events ' ,←-@myEventsFcnII );

optionsIII = odeset ('MaxStep ',10^( -5) ,'Mass ',MII ,'Events ' ,←-@myEventsFcnIII );

% S t t Max Stepsize !

Nästa steg är att skapa huvudprogrammet, det vill säga det som beskrivs i algortimen.

% Main program

N = 100

for i = 1:N

event = ieout (end); switch event

case 1

[t,y,te ,ye ,ie] = ode15s ( @DAEII ,[0 1],y0 , optionsII ); [zout ,sout ,bout ,zvout ,svout ,bvout , tstart , tout

]=←-Updatelists (y,t,zout ,sout ,bout ,zvout ,svout ,bvout ,←-tstart , tout );

ie = ie(end); ieout = [ ieout ie ]; ye=ye(end,:) ; y0 = NyaIV (ye ,ie); case 2

[t,y,te ,ye ,ie] = ode15s ( @DAEII ,[0 1],y0 , optionsIII ); [zout ,sout ,bout ,zvout ,svout ,bvout , tstart , tout

]=←-Updatelists (y,t,zout ,sout ,bout ,zvout ,svout ,bvout ,←-tstart , tout );

ie = ie(end); ieout = [ ieout ie ]; ye=ye(end,:) ; y0 = NyaIV (ye ,ie); case 3

[t,y,te ,ye ,ie] = ode15s (@DAEI ,[0 1],y0 , optionsI ); [zout ,sout ,bout ,zvout ,svout ,bvout , tstart , tout

]=←-Updatelists (y,t,zout ,sout ,bout ,zvout ,svout ,bvout ,←-tstart , tout );

ie = ie(end); ieout = [ ieout ie ]; ye=ye(end,:) ; y0 = NyaIV (ye ,ie); case 4

[t,y,te ,ye ,ie] = ode15s (@DAEI ,[0 1],y0 , optionsI ); [zout ,sout ,bout ,zvout ,svout ,bvout , tstart , tout

]=←-Updatelists (y,t,zout ,sout ,bout ,zvout ,svout ,bvout ,←-tstart , tout );

ie = ie(end); ieout = [ ieout ie ]; ye=ye(end,:) ; y0 = NyaIV (ye ,ie); case 5

[zout ,sout ,bout ,zvout ,svout ,bvout , tstart , tout ]=←-Updatelists (y,t,zout ,sout ,bout ,zvout ,svout ,bvout ,←-tstart , tout );

ie = ie(end); ieout = [ ieout ie ]; ye=ye(end,:) ; y0 = NyaIV (ye ,ie);

end end

Det här följer logiken i Matlab då huvudprogrammet alltid måste komma före funktionerna. Varje case (fall) representerar det stadie som systemet hamnar i efter en given händelse. Om händelsen är 1, så iterar programmet igenom case = 1 och så vidare. Härefter definieras funktionerna: Händelse-funktionerna, ODE-funktionerna och initialvärdesfunktionen.

function [value , isterminal , direction ] = myEventsFcnI (t,y)

% A l l m n n a konstanter rS = 3.1*10^( -3) ; % [m] r0 = 2.5*10^( -3) ; % [m] hS = 2*10^( -2) ; % [m] if y (6) <0 value = [y (2) +(rS -r0)/hS , 1]; isterminal = [1, 0]; direction = [-1, 1]; else value = [1, y (2) -(rS -r0)/hS ]; isterminal = [0, 1]; direction = [1, 1]; end % value = [y (2) +(rS -r0)/hS , y (2) -(rS -r0)/hS ]; % isterminal = [1, 1]; % direction = [-1, 1]; end

function [value , isterminal , direction ] = myEventsFcnII (t,y)

% A l l m n n a konstanter mS = 3*10^( -4) ; % [kg] mB = 4.5*10^( -3) ; % [kg] g = 9.81; % [m/s ^2] JB = 7*10^( -7) ; % [ kgm ] rS = 3.1*10^( -3) ; % [m] lS = 10^( -2) ; % [m] hS = 2*10^( -2) ; % [m] lG = 1.5*10^( -2) ; % [m] cp = 5.6*10^( -3) ; % [N/ rad ] yB = (cp *(y (2) -y (3) )-mB*lG*g)/( JB+mB*lG ^2) ; Z2 = -(mS+mB)*g-mB*lG*yB; Z1 = (cp *(y (3) -y (2) )-rS*Z2 -mB*lS*lG*yB -mB*lS*g)/hS;