Vlakec

Obravnavali bomo gibanje vlakca po železnici. Zaradi enostavnosti bomo modelirali vlakec kar kot materialno točko, čeprav smo s tem zanemarili rotacijski del gibanja vlakca.

Obliko železnice definiramo s parametrično podano funkcijo Tir[u], katere tirnica predstavlja sklenjeno krivuljo v prostoru.

Tir[u_]={Cos[u],Sin[u],Cos[2u]};

ParametricPlot3D[Tir[u],{u,0,2 Pi},Axes->False,Boxed->False,MaxRecursion->0,Mesh->All,PlotPoints->100,PlotStyle->Blue]

V vsaki točki na krivulji izračunamo smeri enotske vektorje v smeri tangente, normale in binormale.

Tangenta[u_]=D[Tir[u],u];

EnotskaTangenta[u_]=Tangenta[u]/Sqrt[Tangenta[u].Tangenta[u]];

EnotskaNormala[u_]=D[EnotskaTangenta[u],u]/Sqrt[(D[EnotskaTangenta[u],u].D[EnotskaTangenta[u],u])];

Binormala[u_]=Cross[EnotskaTangenta[u],EnotskaNormala[u]];

Stezo definiramo, tako da tirnico odebelimo v smeri normale na krivuljo. Vlakec predstavimo z dvema kroglama in pa pravokotnikom.

steza=ParametricPlot3D[Tir[u]+t*EnotskaNormala[u],{u,0,2 Pi},{t,-0.1,0.1},MaxRecursion->0,Mesh->All,PlotPoints->{50,4},Axes->False,Boxed->False,PlotStyle->Blue];

vlakec[u_]:=ParametricPlot3D[Evaluate[Tir[s]+v* EnotskaNormala[s]+0.09Binormala[s]],{s,u-.15,u+.15},{v,-0.1,0.1},PlotStyle->Green];

kolesa[u_]=Graphics3D[{Green,PointSize[0.03],Point[Tir[u-0.07]+0.05Binormala[u]],Point[Tir[u+0.07]+0.05Binormala[u]]}];

Tirnico smo parametrizirali s parametrom u. Če želimo simulirati gibanje vlakca po stezi, moramo izračunati, kako se u spreminja v odvisnosti od časa. Za začetek bomo predpostavili, da se vlakec giblje pod vplivom sile teže in brez trenja. Enačbe gibanja izpeljemo s pomočjo 2. Newtonovega zakona v tangencialni smeri. 

                                                                             m (d^2)s / (dt^2) = <(Subscript[F, g]),(Subscript[e, t])>.
                                                                             
S posrednim odvajanjem pridemo do zveze (d^2)s /(dt^2) = s''(u) (((du/dt))^2) + s'(u) (d^2)u / (dt^2), ki nam enačbo prevede v obliko     

						           m (s''(u) (((du/dt))^2) + s'(u) (d^2)u / (dt^2)) = <(Subscript[F, g]),(Subscript[e, t])>.	                                                                        
                                                                                                                                                                                  
Vzeli bomo, da je g = 10m/(s^2) in m = 1kg in enačbo rešili numerično.

g=10;m=1;

ds[u_]=Sqrt[Tangenta[u].Tangenta[u]];

dds[u_]=ds'[u];

Najprej si bomo pogledali primer, ko vlakec postavimo na vrh železnice in ga rahlo porinemo.

NDSolve[{m*(dds[x[t]]*x'[t]^2+ds[x[t]]*x''[t])=={0,0,-m*g}.EnotskaTangenta[x[t]],x[0]==0,x'[0]==0.2},x,{t,0,50},PrecisionGoal->10];

a[t_]=First[Evaluate[x[t]/.%]];

Gibanje vlakca bo v tem primeru periodično, a se bo ves čas premikal v isto smer na krivulji.

Plot[a[t],{t,0,50}]

Animate[Show[{kolesa[a[t]],steza,vlakec[a[t]]},PlotRange->{{-1.1,1.1},{-1.1,1.1},{-1.15,1.15}},Boxed->False,Axes->False],{{t,0,"Čas (s)"},0,60,ControlType->Trigger,AnimationRate->1}]

Sedaj si poglejmo primer, ko vlakec postavimo nekje na sredo železnice in ga spustimo.

NDSolve[{m*(dds[x[t]]*x'[t]^2+ds[x[t]]*x''[t])=={0,0,-m*g}.EnotskaTangenta[x[t]],x[0]==0.95*2Pi,x'[0]==0},x,{t,0,50},PrecisionGoal->10];

a[t_]=First[Evaluate[x[t]/.%]];

Vlakec bo sedaj periodično nihal med dvema legama na železnici.

Plot[a[t],{t,0,50}]

Animate[Show[{kolesa[a[t]],steza,vlakec[a[t]]},PlotRange->{{-1.1,1.1},{-1.1,1.1},{-1.15,1.15}},Boxed->False,Axes->False],{{t,0,"Čas (s)"},0,60,ControlType->Trigger,AnimationRate->1}]

Če želimo modelirati še silo trenja, moramo v enačbo gibanja na desno stran dodati še člen  

                                                                             m (d^2)s / (dt^2) = <(Subscript[F, g]),(Subscript[e, t])> + (Subscript[F, tr]),
                                                                             
kjer je (Subscript[F, tr]) = -k sgn(ds/dt) (Subscript[F, p]), sila podlage (Subscript[F, p]) pa vsebuje samo komponenti sile železnice na vlakec v normalni ravnini na tirnico. Za koeficient trenja bomo vzeli k = 0,02.

Vlakec postavimo na vrh železnice in ga porinemo dovolj močno, da doseže vrh sosednjega hriba, nato pa se izniha v okolici ravnovesne lege pri u = 3π/2.

k=0.02;

Ukrivljenost[u_]=Sqrt[Cross[Tir'[u],Tir''[u]].Cross[Tir'[u],Tir''[u]]]/Sqrt[Tir'[u].Tir'[u]]^3;

NDSolve[{m*(dds[x[t]]*x'[t]^2+ds[x[t]]*x''[t])=={0,0,-m*g}.EnotskaTangenta[x[t]]-k*Sign[x'[t]]*Sqrt[(m*Ukrivljenost[x[t]]*ds[x[t]]^2*x'[t]^2-{0,0,-m*g}.EnotskaNormala[x[t]])^2
+(-{0,0,-m*g}.Binormala[x[t]])^2],x[0]==0,x'[0]==3},x,{t,0,35},PrecisionGoal->10];

a[t_]=First[Evaluate[x[t]/.%]];

Plot[a[t],{t,0,35}]

Animate[Show[{kolesa[a[t]],steza,vlakec[a[t]]},PlotRange->{{-1.1,1.1},{-1.1,1.1},{-1.15,1.15}},Boxed->False,Axes->False],{{t,0,"Čas (s)"},0,35,ControlType->Trigger,AnimationRate->1}]