#! /usr/bin/env python3 # -*- coding: UTF-8 -*- # FY1001/TFY4145 Mekanisk fysikk hosten 2015. # Oving 6 from numpy import loadtxt,array,append,zeros,arctan,sqrt,pi,cos,sin,arange from matplotlib.pyplot import rc,figure,title,xlabel,ylabel,plot,legend,hold,show,subplot,axvline def main(): # Laste inn data fra fil, hvor de to forste linene ikke inneholder data data=loadtxt('bordtennisball.txt',skiprows=2) c=2.0/3 # Bordtennisballens treghetsmoment I=cr^2 g=9.81 # Tyngdens akselerasjon # Separering av data lest inn fra filen texp=data[:,0] xexp=data[:,1] yexp=data[:,2] # N=Antall datapunkt(tidssteg) N=len(texp) # Definere tomme tabeller for eksperimentell hastighet og energi Vxexp=zeros(N) Vyexp=zeros(N) Vexp=zeros(N) Kexp=zeros(N) Uexp=zeros(N) Eexp=zeros(N) # Beregne hastighet i det første tidssteget Vxexp[0]=(xexp[1]-xexp[0])/(texp[1]-texp[0]) Vyexp[0]=(yexp[1]-yexp[0])/(texp[1]-texp[0]) # Beregne hastigheten for hvert tidssteg(untatt forste og siste) for i in range(1,N-1): Vxexp[i]=(xexp[i+1]-xexp[i-1])/(texp[i+1]-texp[i-1]) Vyexp[i]=(yexp[i+1]-yexp[i-1])/(texp[i+1]-texp[i-1]) # Beregne hastighet i det siste tidssteget Vxexp[N-1]=(xexp[N-1]-xexp[N-2])/(texp[N-1]-texp[N-2]) Vyexp[N-1]=(yexp[N-1]-yexp[N-2])/(texp[N-1]-texp[N-2]) # Beregne vinklen phi og radien r+R for hvert tidssteg phiexp=zeros(N) rexp=zeros(N) for i in range(0,N): phiexp[i]=arctan(xexp[i]/yexp[i]) # Vinkel ift. origo rexp[i]=sqrt(xexp[i]**2+yexp[i]**2) Vexp[i]=sqrt(Vxexp[i]**2+Vyexp[i]**2) Kexp[i]=(c+1)*Vexp[i]**2/(2*100*100) # Kinetisk energi ved rulling Uexp[i]=g*yexp[i]/100 # Potensiell energi i tyngdefeltet Eexp[i]=Kexp[i]+Uexp[i] # Total energi=Kin+Pot # Felles parametere g=9.81 # Tyngdens akselerasjon phimax=pi/2 # Maksimal vinkel dt=0.01 # Tidssteg R=0.785 # Sirkelradius # Parametere for bordtennisball mys=0.25 # Statisk friksjonskoeffisient myk=0.25 # Kinematisk friksjonskoeffisient r=0.019 # Bordtennisballens radius c=2.0/3 # Bordtennisballens treghetsmoment I=cr^2 rpR=r+R # Radius til massesenteret tid=array([0.0]) # Starter tiden ved t=0 V=array([0.015]) # Start-hastighet phi=array([0.047]) # Start-vinkel Omega=array(V/rpR) # Vinkelhastighet til CM, om origo n=array(g*cos(phi)) # Normalkraft A=array(g*sin(phi)/(c+1)) # Akselerasjon langs banen Alpha=array(A/rpR) # Vinkelakselerasjon til CM f=array([0.0]) # Friksjonskraft fmax=array([0.0]) # Maksimal statisk friksjonskraft x=array(rpR*sin(phi)) # CM Posisjon, x-koordinat y=array(rpR*cos(phi)) # CM Posisjon, y-koordinat Vx=array(V*cos(phi)) # CM Hastighet, x-komponent Vy=array(V*sin(phi)) # CM Hastighet, y-komponent jmax=20/dt # Maksimalt antall tidssteg jrulling=0 # Siste steg-nummer for rulling jkontakt=0 # Siste steg-nummer for kontakt j=0 # Nå er alle størrelser vi trenger definert og initialisert. # Vi integrerer bevegelsesligningen ved å iterere while-løkken while (n[j] > 0.0) and (phi[j] < phimax) and (j fmax[j+1]): # Sjekke om skliing er aktuelt f[j+1]=myk*n[j+1] # Ny kinematisk friksjonskraft A[j+1]=g*sin(phi[j+1])-f[j+1] # Ny akselerasjon Alpha[j+1]=A[j+1]/rpR # Ny vinkelakselerasjon else: jrulling=j+1 # Tar vare pa siste steg med rulling jkontakt=j+1 # Tar vare pa siste steg med kontakt tid=append(tid,j*dt) # Ny tid j+=1 # Nytt steg # Nå er vi ferdig med å integrere bevegelsesligningen. # Printer ut interessente vinkler og hastigheter print('Ved rulleslutt er:\n vinkelen: %.1f grader\n vinkelhastigheten: %.1f grader/s\n hastigheten: %.1f m/s'%(phi[jrulling]*180/pi,Omega[jrulling]*180/pi,V[jrulling])) print('Ved kontaktslutt er:\n vinkelen: %.1f grader\n vinkelhastigheten: %.1f grader/s\n hastigheten: %.1f m/s'%(phi[jkontakt]*180/pi,Omega[jkontakt]*180/pi,V[jkontakt])) # Her lager vi plot av ulike størrelser for å sammenligne numeriske og # eksperimentelle data. rc('font',**{'family':'sans-serif','serif':'Computer Modern'}) # Figur bestaende av 2 figurer hoyden og 2 i bredden figure('Figur1') # del-figur 1: Vinkel, vinkelhastighet om origo og normalkraft som funksjon av tid subplot(221) title(r'$\phi,\Omega$'+' og '+r'$N$'+' som funksjon av tid') hold(True) xlabel('t(s)') plot(tid,phi*180/pi,label=r'$\phi({\rm grader})$') plot(tid,Omega*180/pi,label=r'$\Omega({\rm grader}/s)$') plot(tid,n,label='$N(N/kg)$') axvline(tid[jrulling]) # Linje som skiller(tid) rulling fra skliing legend(loc='best') # Delfigur 2: Vinkelhastighet, vinkelakselerasjon og normalkraft som funksjon av vinkel subplot(222) title(r'$\Omega, \alpha$'+' og '+r'$N$'+' som funksjon av vinkel') hold(True) xlabel(r'$\phi({\rm grader})$') plot(phi*180/pi,Omega*180/pi,label=r'$\Omega({\rm grader}/s)$') plot(phi*180/pi,Alpha*180/pi,label=r'$\alpha({\rm grader}/s^2)$ ') plot(phi*180/pi,n,label=r'$N(N/kg)$') axvline(phi[jrulling]*180/pi) # Linje som skiller(vinkel) rulling fra skliing legend(loc='best') # Delfigur 3: Hastighet i x- og y-retning som funksjon av tid subplot(223) title(r'$V_x$'+' og '+r'$V_y$'+' som funksjom av tid') hold(True) xlabel('t(s)') plot(tid,Vx*100,label=r'$V_x(cm/s)$') plot(tid,Vy*100,label=r'$V_y(cm/s)$') plot(texp,Vxexp,label=r'$V_x^{exp}$') plot(texp,Vyexp,label=r'$V_y^{exp}$') axvline(tid[jrulling]) legend(loc='best') # Delfigur 4: Hastighet i x- og y-retning som funksjon av vinkel subplot(224) title(r'$V_x$'+' og '+r'$V_y$'+' som funksjom av vinkel') hold(True) xlabel(r'$\phi({\rm grader})$') plot(phi*180/pi,Vx*100,label=r'$V_x(cm/s)$') plot(phi*180/pi,Vy*100,label=r'$V_y(cm/s)$') plot(phiexp*180/pi,Vxexp,label=r'$V_x^{\rm exp}$') plot(phiexp*180/pi,Vyexp,label=r'$V_y^{\rm exp}$') axvline(phi[jrulling]*180/pi) legend(loc='best') show() if __name__=="__main__": main()