Symulacja rzutów
(Projekt zaliczeniowy z kursu Mathematica I: kurs wstępny)

Mikołaj Gołuński
Zaawansowane Materiały i Nanotechnologia, II rok
Wydział Fizyki, Astronomii i Informatyki Stosowanej
Uniwersytet Jagielloński
Kraków 2010

Wstęp teoretyczny

Rzut(1)

Rzut jest jednym z podstawowych zagadnień związanych z kinematyką. Rozważamy ruch rzuconego w polu grawitacyjnym ciała posiadającego jakąś prędkość początkową.
	Rozróżniamy rzut poziomy, ukośny i pionowy oraz związany z nimi spadek swobodny. Rzuty te różnią się od siebie kątem nachylenia wektora prędkości początkowej ciała.
	Ogólnym sposobem rozwiązania problemu rzutu jest rozłożenie ruchu ciała na składową x-ową i y-ową. Tego typu równania łatwo można rozwiązać, a następnie znaleźć rozwiązanie dla całego ruchu poprzez ich złożenie.

Tarcie(2)

Siła tarcia jest to siła jaką działa ośrodek na poruszające się ciało, w którym ono się porusza. Siła ta przeciwdziała ruchowi ciała.
	W rozważanym przez nas przypadku mamy do czynienia z tarciem powietrza, które można opisać poprzez równanie na opór aerodynamiczny

(Subscript[F, o])=(((Subscript[ρ, p])·(v^2)·S·(Subscript[C, x]))/2)

gdzie
(Subscript[F, o]) – siła oporu
v–prędkość ciała względem ośrodka
(Subscript[ρ, p])– gęstość ośrodka (w naszym przypadku powietrza)
S – pole rzutu ciała na powierzchnię prostopadłą do kierunku ruchu ciała względem ośrodka (w przypadku kuli pole będzie wynosiło ((πr)^2))
(Subscript[C, x])– współczynnik siły oporu zależny od kształtu ciała

Komentarz do symulacji

Sposób przeprowadzenia symulacji

Ciało (kulka) w symulacji zostaje zrzucone z pewnej wysokości i odbija się doskonale sprężyście od "lusterka", którego wysokość i kąt względem podłoża możemy regulować. Wysokość z jakiej kulka jest zrzucana ustalana jest poprzez wybranie prędkości początkowej z jaką ma się ona poruszać w momencie odbicia od "lusterka". Ponadto w przypadku z uwzględnieniem tarcia można wybrać masę oraz promień kulki.  
  	Przy obliczaniu przypadku, w którym nie bierzemy pod uwagę sił tarcia, wyszedłem od sił działających na ciało (w naszym wypadku tylko siła grawitacji) i poprzez przeprowadzanie operacji całkowania, pamiętając o stałych, doszedłem do równania ruchu.
  	Gdy chciałem dodać siły tarcia nie mogłem jednak użyć tego samego sposobu ponieważ siła tarcia zmienia się w czasie w zależności od prędkości ciała. Z tego powodu tarcie uwzględniłem na etapie dodawania prędkości i, w razie potrzeby, z prędkości wyliczałem siłę jaka w danym momencie działa na ciało.
	Ponadto przy rozważaniu przypadku z tarciem trzeba pamiętać o osiąganiu przez ciało tzw. prędkości granicznej tak na odcinku spadku swobodnego jak podczas samego rzutu.

Stałe przyjęte w symulacji

W poniższej symulacji przyjąłem następujące stałe:

(Subscript[ρ, p]) (gęstość powietrza  w temperaturze 25°C i pod ciśnieniem 1atm)=1.185 (kg/(m^3))(3)

(Subscript[C, x]) (współczynnik siły oporu dla kuli)=0.47(2)

g (przyspieszenie ziemskie)=9.80665 (m/(s^2))

Symulacja

Obliczenia

Clear["Global`*"]
(*czyszczenie wszystkich zmiennych z pamięci*)

(*---STAŁE---*)
<<PhysicalConstants`(*załadowanie pakietu stałych fizycznych*)
g=AccelerationDueToGravity*Second^2/Meter;(*podstawienie pod "g" przyspieszenia grawitacyjnego i usunięcie jego wymiaru, żeby nie przeszkadzał w obliczeniach*)

Cx=0.47;(*współczynnik tarcia dla kuli*)
ρp=1.184;(*gęstość powietrza dla temperatury T=25°C i ciśnienia p=1atm*)

(*---BEZ TARCIA---*)
Fx=0; (*siła działajaca na kulkę wzdłuż osi x*)
Fy=0;(*siła działająca na kulkę wzdłuż osi y (poza siłą grawitacji)*)
Fgy[m_]=-m g;(*siła grawitacji działająca na kulkę*)
ax[m_]=Fx/m;(*przyspieszenie kulki wzdłuż osi x*)
ay[m_]=(Fgy[m]+Fy)/m;(*przyspieszenie kulki wzdłuż osi y*)
vx[m_,t_,v0_,α_]=v0*Cos[α Degree]+Integrate[ax[m],t];(*prędkość kulki wzdłuż osi x*)
vyH[m_,t_]=Integrate[ay[m],t];(*prędkość kulki wzdłuż osi y na odcinku, na którym spada swobodnie*)
vy[m_,t_,v0_,α_]=v0*Sin[α Degree]+Integrate[ay[m],t];(*prędkość kulki wzdłuż osi y po odbiciu się jej od "lusterka"*)
x[m_,t_,v0_,α_]=Integrate[vx[m,t,v0,α],t];(*równanie ruchu kulki wzdłuż osi x*)
y[m_,t_,v0_,h_,α_]=h+Integrate[vy[m,t,v0,α],t];(*równanie ruchu kulki wzdłuż osi y po odbiciu się jej od "lusterka"*)
H[v0_,h_]=v0^2/(2*g)+h;(*wysokość z jakiej musi spadać kulka aby przy odbiciu od "lusterka" osiągnęła prędkość v0*)
yH[m_,t_,v0_,h_]=H[v0,h]+Integrate[vyH[m,t],t];(*równanie ruchu kulki wzdłuż osi y na odcinku, na którym spada swobodnie*)
czasH[v0_,h_]=Sqrt[2*(H[v0,h]-h)/g];(*czas swobodnego spadku kulki*)
czas[m_,v0_,h_,α_]=Solve[y[m,t,v0,h,α]==0,t][[2,1,2]];(*czas od odbicia się kulki od "lusterka" do "wylądowania"*)
zasieg[m_,v0_,h_,α_]=x[m,czas[m,v0,h,α],v0,α];(*zasięg rzutu*)
hmax[m_,t_,v0_,h_,α_]:=FindMaximum[y[m,t,v0,h,α],t][[1]];(*wysokość maksymalna*)
tmax[m_,t_,v0_,h_,α_]:=FindMaximum[y[m,t,v0,h,α],t][[2,1,2]];(*czas po jakim kulka osiągnie wysokość maksymalną*)

xk[m_,t_,v0_,h_,α_]=Piecewise[{{0,t≤czasH[v0,h]},{x[m,t-czasH[v0,h],v0,α],t>czasH[v0,h]}}];(*równanie ruchu kulki wzdłuż osi x dla całego toru ruchu*)
yk[m_,t_,v0_,h_,α_]=Piecewise[{{yH[m,t,v0,h],t≤czasH[v0,h]},{y[m,t-czasH[v0,h],v0,h,α],t>czasH[v0,h]}}];(*równanie ruchu kulki wzdłuż osi y dla całego toru ruchu*)
czask[m_,v0_,h_,α_]=czasH[v0,h]+czas[m,v0,h,α];(*całkowity czas trwania ruchu kulki*)

(*---TARCIE---*)
Fo[r_,v_]=(ρp v^2 Pi r^2 Cx)/2;(*siła tarcia działająca na kulkę*)
ao[r_,v_,m_]=Fo[r,v]/m;(*przyspieszenie kulki wywołane siłą tarcia*)
vo[r_,v_,m_,t_]=Integrate[ao[r,v,m],t];(*prędkość kulki wywołana siłą tarcia*)
voyH[r_,m_,t_]=vyH[m,t]+vo[r,vyH[m,t],m,t];(*prędkość kulki wzdłuż osi y na odcinku, na którym spada swobodnie, z uwzględnienim tarcia*)
yoyH[r_,m_,t_,v0_,h_]=oH[r,m,v0,h]+Integrate[voyH[r,m,t],t];(*równanie ruchu dla kulki spadającej swobodnie z uwzględnieniem tarcia*)
toHmin[r_,m_]=Solve[D[voyH[r,m,t],t]==0,t][[2,1,2]];(*czas po jakim kulka spadając swobodnie osiągnie prędkość graniczną*)
hoHmin[r_,m_,v0_,h_]=yoyH[r,m,toHmin[r,m],v0,h];(*wysokość na jakiej kulka spadając swobodnie osiągnie prędkość graniczną*)
voyHmin[r_,m_]=voyH[r,m,toHmin[r,m]];(*prędkość graniczna jaką osiągnie kulka spadając swobodnie*)
yoyHmin[r_,m_,t_]=Integrate[voyHmin[r,m],t];(*równanie ruchu dla spadającej swobodnie kulki, która uzyskała prędkość graniczną*)
yoHk[r_,m_,t_,v0_,h_]=Piecewise[{{yoyH[r,m,t,v0,h],t<toHmin[r,m]},{hoHmin[r,m,v0,h]+yoyHmin[r,m,t-toHmin[r,m]],t≥toHmin[r,m]}}];(*równanie ruchu dla kulki spadającej swobodnie uwzględniające tarcie i prędkość graniczną*)
toH[r_,m_,h_,H_]=Solve[-Integrate[voyH[r,m,t],t]==H-h,t][[4,1,2]];(*czas jaki jest potrzebny na uzyskanie prędkości, przy odbiciu od "lusterka", równej v0*)
oH[r_,m_,v0_,h_]=Re[Solve[v0==voyH[r,m,toH[r,m,h,H]],H][[2,1,2]]];(*wysokość z jakiej trzeba zrzucić kulkę aby przy odbiciu od "lusterka" uzyskała prędkość v0, z uwzględnieniem tarcia*)
czasoH[r_,m_,v0_,h_]=Piecewise[{{Solve[yoHk[r,m,t,v0,h][[1,1,1]]==h,t][[2,1,2]],v0<-voyHmin[r,m]},{Solve[yoHk[r,m,t,v0,h][[1,2,1]]==h,t][[1,1,2]],v0≥-voyHmin[r,m]}}];(*czas spadku kulki od jej upuszczenia do odbicia się od "lusterka"*)

vox[r_,m_,t_,v0_,α_]=Piecewise[{{vx[m,t,v0,α]-vo[r,vx[m,t,v0,α],m,t],vx[m,t,v0,α]>=vo[r,vx[m,t,v0,α],m,t]},{0,vx[m,t,v0,α]<vo[r,vx[m,t,v0,α],m,t]}}];
(*prędkość kulki wzdłuż osi x po odbiciu się jej od "lusterka" z uwzględnieniem tarcia*)
tox[r_,m_,v0_,α_]=Solve[vox[r,m,t,v0,α][[1,1,1]]==0,t][[1,1,2]];(*czas po jakim kulka przestanie się poruszać wzdłuż osi x*)
xox[r_,m_,t_,v0_,α_]=Integrate[vox[r,m,t,v0,α],t];(*równanie ruchu wzdłuż osi x dla kulki, która odbiła się od "lusterka", z uwzględnienim tarcia*)
ox[r_,m_,t_,v0_,α_]=Piecewise[{{xox[r,m,t,v0,α],vox[r,m,t,v0,α]>0},{xox[r,m,tox[r,m,v0,α],v0,α],vox[r,m,t,v0,α]<=0}}];(*równanie ruchu wzdłuż osi x dla kulki, która odbiła się od "lusterka", z uwzględnieniem tarcia i okresu, w którym prędkość kulki wzdłuż osi x równa się zero*)
voy[r_,m_,t_,v0_,α_]=vy[m,t,v0,α]+vo[r,vy[m,t,v0,α],m,t];(*prędkość kulki wzdłuż osi y, po odbiciu się jej od "lusterka" z uwzględnieniem tarcia*)
oy[r_,m_,t_,v0_,h_,α_]=h+Integrate[voy[r,m,t,v0,α],t];(*równanie ruchu kulki wzdłuż osi y, po odbiciu się jej od "lusterka", z uwzględnieniem tarcia*)
tomin[r_,m_,v0_,α_]=Solve[D[voy[r,m,t,v0,α],t]==0,t][[2,1,2]];(*czas, po którym kulka poruszając się wzdłuż osi y, po odbicu się od "lusterka", osiągnie prędkość graniczną*)
voymin[r_,m_,v0_,α_]=voy[r,m,tomin[r,m,v0,α],v0,α];(*prędkość graniczna kulki poruszającej się wzdłuż osi y, po odbiciu się od "lusterka"*)
oymin[r_,m_,t_,v0_,α_]=Integrate[voymin[r,m,v0,α],t];(*równanie ruchu wzdłuż osi y dla kulki, która odbiła się od "lusterka" i osiągnęła prędkość graniczną*)
ho[r_,m_,v0_,h_,α_]=oy[r,m,tomin[r,m,v0,α],v0,h,α];(*wysokość, na której kulka osiągnie prędkość graniczną*)
koy[r_,m_,t_,v0_,h_,α_]=Piecewise[{{oy[r,m,t,v0,h,α],t<tomin[r,m,v0,α]},{ho[r,m,v0,h,α]+oymin[r,m,t-tomin[r,m,v0,α],v0,α],t≥tomin[r,m,v0,α]}}];(*równanie ruchu wzdłuż osi y dla kulki, która odbiła się od lusterka, z uwzględnieniem tarcia i prędkości granicznej*)
toy[r_,m_,v0_,h_,α_]=Solve[oy[r,m,t,v0,h,α]==0,t][[2,1,2]];(*czas, po którym kulka "wylądowałaby" gdyby uwzględniać tarcie, ale bez prękdości granicznej*)
voyt[r_,m_,v0_,h_,α_]=voy[r,m,toy[r,m,v0,h,α],v0,α];(*prędkość jaką uzyskałaby kulka w momencie "lądowania" gdyby uwzględnić tarcie, ale bez prędkości granicznej*)

oczas[r_,m_,v0_,h_,α_]=Piecewise[{{Solve[oy[r,m,t,v0,h,α]==0,t][[3,1,2]],Abs[voyt[r,m,v0,h,α]]<Abs[voymin[r,m,v0,α]]},{Solve[ho[r,m,v0,h,α]+oymin[r,m,t-tomin[r,m,v0,α],v0,α]==0,t][[1,1,2]],Abs[voyt[r,m,v0,h,α]]≥Abs[voymin[r,m,v0,α]]}}];(*czas po jakim kulka "wyląduje" po odbiciu się od "lusterka" uwzględniając tarcie i prędkość graniczną*)
ozasieg[r_,m_,v0_,h_,α_]=ox[r,m,Abs[oczas[r,m,v0,h,α]],v0,α];(*zasięg rzutu z uwzględnieniem tarcia*)
ohmax[r_,m_,t_,v0_,h_,α_]:=FindMaximum[koy[r,m,t,v0,h,α],{t,0}][[1]];(*wysokość maksymalna rzutu z uwzględnieniem tarcia*)
otmax[r_,m_,t_,v0_,h_,α_]:=FindMaximum[oy[r,m,t,v0,h,α],{t,0}][[2,1,2]];(*czas po jakim kulka osiągnie wysokość maksymalną przy uwzględnieniu tarcia*)

oxk[r_,m_,t_,v0_,h_,α_]=Piecewise[{{0,t<czasoH[r,m,v0,h]},{ox[r,m,t-czasoH[r,m,v0,h],v0,α],t>=czasoH[r,m,v0,h]}}];
(*równanie ruchu dla kulki dla pełnego ruchu wzdłuż osi x, z uwzględnieniem tarcia*)
koyk[r_,m_,t_,v0_,h_,α_]=Piecewise[{{yoHk[r,m,t,v0,h],t<czasoH[r,m,v0,h]},{koy[r,m,t-czasoH[r,m,v0,h],v0,h,α],t>=czasoH[r,m,v0,h]}}];(*równanie ruchu dla kulki dla pełnego ruchu wzdłuż osi y, z uwzględnieniem tarcia*)
oczask[r_,m_,v0_,h_,α_]=czasoH[r,m,v0,h]+Abs[oczas[r,m,v0,h,α]];(*całkowity czas lotu kulki z uwzględnieniem tarcia*)

(*---LUSTERKO---*)
x1[α_,r_]=r (-Cos[α Degree]);(*położenie x-owe pierwszego punktu "lusterka"*)
x2[α_,r_]=r Cos[α Degree];(*położenie x-owe drugiego punktu "lusterka"*)
y1[h_,α_,r_]=h+(Sqrt[(r^2)-(0^2)+2 r 0 (-Cos[α Degree])-(r^2) (((-Cos[α Degree]))^2)]);(*położenie y-owe pierwszego punktu "lusterka"*)
y2[h_,α_,r_]=h-(Sqrt[(r^2)-(0^2)+2 r 0 Cos[α Degree]-(r^2) ((Cos[α Degree])^2)]);(*położenie y-owe drugiego punktu "lusterka"*)
lusterko[h_,α_,r_]=Graphics[Line[{{x1[α,r],y1[h,α,r]},{x2[α,r],y2[h,α,r]}}]];(*funkcja rysująca "lusterko"*)

(*---WYKRESY---*)
beztarcia[v0_,h_,α_,T_]:=
ParametricPlot[
{xk[m,t,v0,h,α],yk[m,t,v0,h,α]},
{t,0,T},
AxesOrigin->{-0.1*zasieg[m,v0,h,α],0},
ImageSize->Medium,
PlotRange->{{-0.1*zasieg[m,v0,h,α],1.1*zasieg[m,v0,h,α]},{0,H[v0,h]}}](*wykres toru ruchu kulki bez uwzględnienia tarcia powietrza*)

ztarciem[r_,m_,v0_,h_,α_,T_]:=
ParametricPlot[
{oxk[r,m,t,v0,h,α],koyk[r,m,t,v0,h,α]},
{t,0,T},
AxesOrigin->{-0.1*ozasieg[r,m,v0,h,α],0},
ImageSize->Medium,
PlotRange->{{-0.1*ozasieg[r,m,v0,h,α],1.1*ozasieg[r,m,v0,h,α]},{0,oH[r,m,v0,h]}}](*wykres toru ruchu kulki z uwzględnieniem tarcia powietrza*)

(*---PEŁNE ANIMACJE---*)
calebeztarcia=
Manipulate[
Column[
{
Show[
beztarcia[v0,h,90-2*α,T],
lusterko[h,α,0.1*zasieg[m,v0,h,90-2*α]]
],
Grid[
{{N[H[v0,h]]"m"," - wysokość, z której spada kulka"},
{N[zasieg[m,v0,h,90-2*α]]"m"," - zasięg rzutu"},
{N[hmax[m,t,v0,h,90-2*α]]"m"," - wysokość maksymalna kulki"},
{N[czas[m,v0,h,90-2*α]]"s"," - czas od odbicia od \"lusterka\" do upadku na ziemię"}},
Alignment->Left,
Frame->{False,All}
]
}
],
{{T,czask[m,v0,h,90-2*α],"czas [s]"},0.0001,czask[m,v0,h,90-2*α]},
{{v0,10,"prędkość przy odbiciu [m/(s^2)]"},0.01,10},
{{h,0,"położenie \"lusterka\" [m]"},0,10},
{{α,22.5,"kąt odbicia [deg]"},0,45},
ControlPlacement->Left,
Paneled->False
];(*połączenie wszystkich części animacji bez tarcia*)

caleztarciem=Manipulate[{
If[v0>-voyHmin[r,m],v0=-voyHmin[r,m],v0=v0],
Column[
{
Show[
ztarciem[r,m,v0,h,90-2*α,T],
lusterko[h,α,0.1*ozasieg[r,m,v0,h,90-2*α]]
],
Grid[
{{N[oH[r,m,v0,h]]"m"," - wysokość, z której spada kulka"},
{N[ozasieg[r,m,v0,h,90-2*α]]"m"," - zasięg rzutu"},
{N[ohmax[r,m,t,v0,h,90-2*α]]"m"," - wysokość maksymalna kulki"},
{N[Abs[oczas[r,m,v0,h,90-2*α]]]"s"," - czas od odbicia od \"lusterka\" do upadku na ziemię"}},
Alignment->Left,
Frame->{False,All}
]
}
]
}[[2]],
{{T,oczask[r,m,v0,h,90-2*α],"czas [s]"},0.0001,oczask[r,m,v0,h,90-2*α]},
{{v0,-voyHmin[r,m],"prędkość przy odbiciu [m/(s^2)]"},0.01,-voyHmin[r,m]},
{{h,0,"położenie \"lusterka\" [m]"},0,10},
{{α,22.5,"kąt odbicia [deg]"},0,45},
{{m,100,"masa kulki [kg]"},0.01,200},
{{r,1,"promień kulki [m]"},0.01,2},
ControlPlacement->Left,
Paneled->False
];(*połączenie wszystkich części animacji z tarciem*)

(*---RYSOWANIE---*)
rysuj=Manipulate[If[tarcie,caleztarciem,calebeztarcia],{{tarcie,False},{True,False}},AppearanceElements->None]; (*funkcja rysująca wizualizację*)

Ilustracja

rysuj

Literatura

(1) http://www.fizykon.org/kinematyka/rzuty_wzory.htm
(2) http://en.wikipedia.org/wiki/Drag_%28 physics %29
(3) W. Mizerski Małe tablice fizyczno-astronomiczne, wyd. II, Adamantan, Warszawa 2004