TFY4106 Fysikk. Institutt for fysikk, NTNU. Våren 2016. Numerisk øving med obligatorisk innlevering på itslearning. Veiledning: 1. mars - 7. april. Innleveringsfrist: 10. april kl 23.59.

Matematisk/fysisk pendel med friksjon

I denne oppgaven skal vi prøve å reprodusere bevegelsen til en pendel ved å løse bevegelsesligningen (dvs Newtons 2. lov) numerisk.

Figuren til venstre viser en matematisk pendel, bestående av ei kule (punktmasse) med masse m i enden av ei masseløs stang med lengde L. Stanga kan svinge omkring festepunktet (A), slik at kula følger en sirkelbane med radius L. Dersom kula svinger fritt, uten å påvirkes av andre ytre krefter enn tyngdekraften og ”stangdraget”, er bevegelsen bestemt av ligningen

\frac{d^{2}θ}{dt^{2}} + \frac{g}{L} sin θ = 0.

(N2 tangentielt til sirkelbanen, der vi har brukt sammenhengen v = Lθ̇ mellom hastighet og vinkelhastighet ved sirkelbevegelse.)

For små utsving fra likevekt, |θ| ≪ 1, kan sin θ erstattes med θ, og vi har en enkel harmonisk oscillator, med kjent løsning, og med vinkelfrekvens (”egenfrekvens”) ω_{0} = \sqrt{g/L}. Men for større utsving fra likevekt må vi beholde sin θ, og da kan ligningen ovenfor ikke løses analytisk. Numerisk er det imidlertid ingen problemer! En enkel og intuitiv numerisk måte å bestemme θ(t) på er den såkalte Euler-metoden. N2 kan skrives på formen

dv = a dt = \frac{F}{m} dt.

Det betyr at hastigheten ved tidspunktet t+dt kan bestemmes dersom vi kjenner hastigheten ved tidspunktet t:

v(t + dt) = v(t) + a dt = v(t) + dv = v(t) + \frac{F}{m} dt.

Med et endelig tidssteg Δt blir ligningen

v(t + Δt) = v(t) + Δv = v(t) + \frac{F (t)}{m} Δt

bare tilnærmet riktig, men formodentlig en riktig god tilnærmelse dersom Δt velges tilstrekkelig liten. Samme oppskrift kan vi i neste omgang bruke for å bestemme posisjonen, eller som her, vinkelen θ(t). Vi har v = ds/dt = L dθ/dt, dvs

dθ = \frac{v}{L} dt,

og med endelig tidssteg gir dette

θ(t + Δt) = θ(t) + Δθ = θ(t) + \frac{v(t)}{L} Δt.

Hvis vi nå, som her, kjenner θ(t = 0) og v(t = 0), dvs θ(0) = θ_{0} og v(0) = v_{0}, kan ligningene over brukes til å finne θ(Δt) og v(Δt), deretter θ(2Δt) og v(2Δt) osv. I vårt konkrete tilfelle er det i første rekke tyngdens tangentialkomponent som bestemmer akselerasjonen, og dermed hastigheten tangentielt,

a = \frac{F}{m} = −g sin θ.

Bevegelsen til en virkelig pendel bremses av ulike former for friksjon: Luftmotstanden vil med brukbar tilnærmelse være proporsjonal med kvadratet av massens hastighet, f_{D} = Dv^{2}, der koeffisienten D har enheten kg/m. I tillegg er det friksjon i festeanordningen ved rotasjonsaksen (A). Vår pendelstang er festet til en ring med indre radius R litt større enn radien r til akslingen, som rett og slett er ei metallstang. Vi antar at ringen og metallstanga har kontakt i ett punkt, og at det der virker en friksjonskraft f_{μ} = μN. Her er μ kinetisk friksjonskoeffisient mellom ring og metallstang, og N er normalkraften som virker mellom disse to. Normalkraften N fastlegges ved å benytte N2 normalt på sirkelbanen:

N − mg cos θ = m \frac{v^{2}}{L} ,
N = m (g cos θ + \frac{v^{2}}{L} ) ,

så vi ser at N, og dermed f_{μ}, avhenger av pendelens utsving og hastighet. Pendelens bevegelse er filmet med hurtigkamera (100 bilder pr sekund) og analysert i Tracker. Tre filmklipp er lagret i video1.mp4, video2.mp4 og video3.mp4. Tilsvarende tider og posisjoner er lagret i track1.txt, track2.txt og track3.txt, med fire kolonner for tid t, horisontal posisjon x (positiv mot høyre), vertikal posisjon y (positiv oppover) og vinkelen θ (relativt positiv x-akse, og med positiv vinkelendring for rotasjon mot klokka). Enheter i disse filene er s for t, mm for x og y, og grader for θ. Samsvar med vinkelen θ i figuren ovenfor oppnås dermed ved å legge 90 grader til verdiene i de tre txt-filene. Og omregning fra grader til radianer oppnås ved å multiplisere med faktoren π/180. I video1.mp4 roterer pendelen flere ganger rundt i samme retning før den etter hvert mister så mye mekanisk energi at den begynner å svinge fram og tilbake. Merk at vinkelen θ øker kontinuerlig så lenge pendelen fortsetter å rotere mot klokka. Maksimalverdien 1883.1 grader inntreffer etter ca 6.9 s.

Oppgaver:

• Lag et (eller tre) program i Matlab (eller Python) som beregner θ(t). Tilpass verdier for startvinkel og -hastighet, samt parametre μ og D for friksjonskreftene slik at de beregnede banene med god tilnærmelse reproduserer pendelens bevegelse i de tre filmklippene. Det fysiske systemet er det samme i de tre filmklippene, så det bør være mulig å få godt samsvar mellom eksperiment og numerikk for alle tre baner, med samme verdi for D og samme verdi for μ.

• Kall programfilene dine brukernavn_pendel1.m og tilsvarende med pendel2 og pendel3μ og D hvis du lager et program for hvert filmklipp. (Hvis du programmerer i Python, vil programfilene ha endelse .py.) Når du har et program som fungerer bra for ett filmklipp, er det ikke så store endringer som skal til for å få det til å fungere for de to andre filmklippene.

• Framstill resultatene i tre figurer, en figur for hvert filmklipp, der hver figur sammenligner eksperimentell og beregnet vinkel θ, i enheten grader, som funksjon av tiden t. Figurene skal ha tekst på aksene, ’t (s)’ på horisontal akse og ’θ (grader)’ på vertikal akse. Figurene lagres som PDF-filer, med navn brukernavn_FIG1.pdf og tilsvarende for FIG2 og FIG3.

• Lever programfil(er) og figurfiler på itslearning. Der leverer du også ei fil (PDF produsert med LaTeX, Word e.l.) der du i korte trekk beskriver det du har gjort og diskuterer resultatene: Er det tørr friksjon eller luftmotstand som er den dominerende friksjonskraften i dette eksperimentet? Er perioden (svingetiden) for de svakt dempede svingningene (i video2 og video3) som forventet? (Dvs, som for en matematisk pendel med lengde L.) Dersom en eksponentielt avtagende amplitude θ_{0} exp(−γt) skulle tilpasses de svakt dempede svingningene, hva blir da omtrent dempingsfaktoren γ? Gir video2 og video3 samme verdi for γ? Hvis ikke, hva kan årsaken være?

• Gjør øvingen alene eller i samarbeid med en medstudent. Hvis dere jobber to og to sammen, kan dere gjerne levere en felles besvarelse på itslearning. Skriv da tydelig navn på begge to på det som leveres inn.

Tips:

• Det er flere metoder som kan brukes i Matlab for å lese inn data fra fil. En kommando som fungerer er

[tideksp,xeksp,yeksp,thetaeksp] = textread('track1.txt','%f %f %f %f','headerlines',2);

Da blir tideksp en tabell med verdiene i 1. kolonne i fila track1.txt osv. Det er 2 linjer med tekst i starten av fila, derfor 'headerlines',2.

• Jeg anbefaler å snu y-aksen på de målte verdiene, samt å gjøre om fra mm til m og fra grader til radianer. (Selv om vinkler skal angis i grader i figurene til slutt.)

• Programmet kan finne antall elementer i en tabell med kommandoen length:

N = length(t)

• I forkant av en for-løkke lønner det seg å preallokere (initialisere) tabeller for de størrelsene som beregnes (oppdateres) i for-løkka. Eksempel:

dt = 0.00001;
t = linspace(0:dt:9.65);
theta = 0*t;

Her blir theta en tabell med like mange elementer som t, og med verdi null for alle elementer i utgangspunktet.

• Hastigheter regnes enklest ut fra posisjoner med kommandoen diff, som ganske enkelt gir differansen mellom to påfølgende elementer i tabellen. Eksempel:

Vxeksp = diff(xeksp)/dteksp;
Vxeksp(Neksp) = Vxeksp(Neksp-1);

Her er dteksp det eksperimentelle tidssteget. Den 2. linja tas med for at Vxeksp skal få like mange elementer som xeksp.

• Det er selvsagt viktig at alle bidrag til det totale dreiemomentet på pendelen, relativt rotasjonsaksen A, får riktig fortegn i programmet ditt. Det er klart at både f_{μ} (tørr friksjon mellom aksling og ”pendelring)”) og f_{D} (luftmotstand) må ha en retning som bremser bevegelsen. Dette kan en forsikre seg om ved å inkludere faktoren -sign(v(n)) i uttrykkene for f_{μ} og f_{D} i programmet, og deretter for eksempel oppdatere hastigheten slik:

v(n+1) = v(n) + (dt/(m*L))*(fG(n)*L + fD(n)*L + fmu(n)*r);

Her er L armen til både luftmotstanden fD(n) og tyngdekraften fG(n) (i tidssteg n), mens r er armen til fmu(n). Merk at tyngdekraften f_{G} ikke alltid virker bremsende, så faktoren -sign(v(n)) kan ikke inkluderes i f_{G}.

• Eksempel på plotting i Matlab:

figure(1);
plot(t,theta,tideksp,thetaeksp);
xlabel('t (s)');
ylabel('\theta (grader)');

• Aller siste del av video1 kan være vanskelig å beskrive godt. Det skyldes ”slingring” normalt pendelens plan, noe som ikke er så godt synlig i filmen. Det betyr at du kan være fornøyd hvis du klarer å få godt samsvar mellom eksperiment og numerikk fram til ca 7.0 eller 7.5 sekunder.

Opplysninger:

• Pendelen består av ei tynn metallstang med 4 lodd, hver med masse 50 g, festet nær stangas ende. Akslingen er ei metallstang med radius r = 8 mm. Avstanden fra akslingens senterakse til de 4 loddenes massesenter er målt til L = 636 mm. Dette er strengt tatt en fysisk pendel, noe vi tar hensyn til når friksjonskreftene tas i betraktning. (f_{μ} har arm r, mens de andre kreftene har arm L; derfor må Newtons 2. lov for rotasjon brukes.) Vi har imidlertid ikke tatt hensyn til at pendelstanga har en viss masse, og at de fire loddene har en viss utstrekning. Det kan derfor tenkes at du i beregningene dine må variere pendellengden L litt omkring den målte verdien for å få best mulig samsvar mellom eksperiment og beregninger.

• De tre filmklippene starter alle med t = 0. De to første har 0.01 s mellom bildene, det tredje har 0.05 s mellom bildene.