Analysis and Solution of the Mechanical Systems with a Varying Mass in Teaching Physics

František Látal, Jan Říha, Lukáš Richterek, David Smrčka
Department of Experimental Physics, Faculty of Science, Palacký University in Olomouc, Czech Republic

Join the Conversation #WolframTechConf

Introduction

Palacký University

The Palacký University in Olomouc is the second oldest university in the Czech lands (after Charles University in Prague).
The history of higher education in Olomouc began in 1566.
With almost 24,000 undergraduate students on eight faculties, it provides quality education and excellent background for research in the wide range of academic disciplines.

Faculty of Science

The Faculty of Science is a research-based faculty that provides professional education in Mathematics, Physics, Chemistry, Biology, and Earth Sciences.
The Faculty of Science trains also future secondary school teachers of sciences, mainly in full-time courses.
Currently, the Faculty of Science has more than 4000 students.

Mechanical Systems with a Varying Mass

Variable Mass Systems

With using the software Mathematica, we solved the problem in which the mass of the body changes during the motion, that is, m is a function of t, i.e. m(t).


We choose three examples of this  variable mass system:

	- 	Condensation of a Water Droplet
	
	- 	The Rocket Equation
	
	- 	Buquoy´s Problem

Condensation of a Water Droplet

A dust particle of negligible mass being under the influence of gravitation begins to fall at time t = 0 by saturated water vapor. The steam condensates with constant rate and a water droplet of steadily increasing mass.

We solved this nonlinear differential equation:

x(t)x''(t) + (((x'(t)))^2) - gx(t) = 0

Initial conditions:

(* m0 = 0.000001; (* just some small value to avoid the singularity *) *)
(Subscript[x, 0])= 0.000001; (* initial position – just some small value to avoid the singularity *)
(Subscript[v, 0])= 0; (* initial velocity *)
tM = 5; (* maximum time for numerical integration *)
g = 9.81; (* acceleration of gravity *)

Solution of the nonlinear differential equation with using NDSolve:

sd=NDSolve[{x''[t]*x[t] + x'[t]^2-g*x[t]== 0,x[0]==(Subscript[x, 0]),x'[0]==(Subscript[v, 0])},x,{t,0,tM}];
xm=NMaximize[{Evaluate[x[t]/.sd[[1]]],0<t≤ tM},t][[1]]; (* Find maximum for PlotRange *)
Plot[Evaluate[x[t]/.sd],{t,0,tM},PlotRange->{Automatic,{-5,1.1*xm}},Filling-> Axis,BaseStyle->{FontFamily->"Times",FontSize-> 14},FrameLabel->{"t","x"},GridLines->Automatic,PlotStyle->{Thickness[0.005]},Frame->True]

Comparison with analytic solution x = gt2/6:

Plot[{Evaluate[x[t]/.sd],g*t^2/6},{t,0,tM},PlotRange->{Automatic,Automatic},Filling-> Axis,BaseStyle->{FontFamily->"Times",FontSize-> 14},AxesLabel->{"t","x"},GridLines->Automatic,PlotStyle->{Thickness[0.005]}]

We can also use a simple Euler-Cromer method with the step h = 0.05 s (and compare with the analytic solution and play with the step size)

h=0.01;
xe[0]=(Subscript[x, 0]);
ve[0]=(Subscript[v, 0]);
For[i=0,i≤ tM/h,i++,
ve[i+1] = ve[i]+(g-ve[i]^2/xe[i])*h;
xe[i+1]=xe[i]+ve[i+1]*h];
graf1e=ListPlot[Table[{i*h,xe[i]},{i,0,tM/h}],AxesStyle->{Arrowheads@.03,Arrowheads@.03},AxesLabel->{"t","x(t)"},BaseStyle->{FontFamily->"Times",FontSize-> 12},GridLines->Automatic,PlotStyle->{Thickness[0.005],Darker[Red]}];
(* graf1=ListLinePlot[Table[{i*h,xe[i]},{i,0,(Subscript[t, m])/h}],AxesStyle->{Arrowheads@.03,Arrowheads@.03},AxesLabel->{"t","x(t)"},BaseStyle->{FontFamily->"Times",FontSize-> 12},GridLines->Automatic,PlotStyle->{Thickness[0.005]}] *)
graf2e=Plot[g*t^2/6,{t,0,tM},AxesStyle->{Arrowheads@.03,Arrowheads@.03},AxesLabel->{"t","x(t)"},BaseStyle->{FontFamily->"Times",FontSize-> 12},GridLines->Automatic,PlotStyle->{Thickness[0.005],Darker[Green]}];
Show[graf1e,graf2e]

Finally, we can show the “droplet” on the chart (increasing its size during its flight):

Manipulate[Show[Plot[g*t^2/6,{t,0,time},PlotStyle->Darker[Green],PlotRange-> {{0,1.1*tM},{0,1.1*Max[Table[xe[i],{i,0,tM/h}]]}}],Graphics[{PointSize[0.01+0.001*xe[Round[time/h]]],Darker[Red],Point[{time,g*time^2/6}]}]],{time,0.0001,tM,Appearance->"Open"}]

The Rocket Equation

A rocket of initial mass m(0) per unit time expels the quantity of gas Δm/Δt > 0 with the constant velocity v(0). We look for equation of motion.

Initial conditions :

r0=6378*^3;(* initial distance from the center of the earth r=Rz in meters *)
tr[0]=0;
xr[0]=r0;
vr[0]=0;(* initial velocity *)
g=9.81; 
q=q0=120;(* mass of burned fuel per one second *)
m[0]=13*^3;(* initial mass of rocket *)
mk=4.25*^3;(* finally mass of rocket *)
c1=2.2*^3;(* velocity of burning fuels *)
ar[0]=(c/m[0])*q-g*r0*r0/x[0]/x[0];

Euler method

dt=0.05; (* integration step *)
For[i=0,xr[i]-r0≥0,i++,
(* after burning the fuel weight change equal to 0 *)
If[ m[i]<mk,q=0,q=q0];
m[i+1]=m[i]-q*dt;
ar[i]=(c1/m[i])*q-g*r0^2/xr[i]^2;
vr[i+1]=vr[i]+ar[i]*dt;
xr[i+1]=xr[i]+vr[i]*dt;
tr[i+1]=tr[i]+dt;
]
mez=i;

graf1r=ListPlot[Table[{j*dt,xr[j]},{j,0,mez}],AxesStyle->{Arrowheads@.03,Arrowheads@.03},BaseStyle->{FontFamily->"Times",FontSize-> 14},AxesLabel->{"t/s","h"},GridLines->{Automatic,{r0,1.01*r0,1.02*r0,1.03*r0}},PlotStyle->{Darker[Blue,0.5],Thickness[0.008]},Ticks->{Automatic,{{r0,"(Subscript[R, z])",.05},{1.01*r0,"1.01(Subscript[R, z])",.05},{1.02*r0,"1.02(Subscript[R, z])",.05},{1.03*r0,"1.03(Subscript[R, z])",.05}}}];
graf2r=ListPlot[Table[{tr[j],vr[j]},{j,0,mez}],BaseStyle->{FontFamily->"Times",FontSize-> 14},Frame->True,FrameLabel->{"t/s",Row[{"v","/m","·","s",""}]},GridLines->Automatic,PlotStyle->{Darker[Red,0.4],Thickness[0.008]},FrameTicks-> {Automatic,Table[{1000*i,i},{i,-2,1}]}];
graf3r=ListPlot[{Table[{tr[j],g*r0^2/xr[j]^2},{j,0,mez}],Table[{tr[j],g},{j,0,mez}]},PlotRange->{Automatic,{0.92*g,1.02*g}},BaseStyle->{FontFamily->"Times",FontSize-> 14},Frame ->True,FrameLabel->{"t/s",Row[{"|(Subscript[a, g])|","/m","·","s",""}]},GridLines->Automatic,PlotStyle->{Thickness[0.008]},Epilog->{Text[Style["g",Italic],{tr[mez]/2,1.005*g}]}];
graf4r= ListPlot[Table[{tr[j],m[j]},{j,0,mez}],PlotRange->{Automatic,{0,13800}},BaseStyle->{FontFamily->"Times",FontSize-> 14},AxesLabel->{"t/s","m/t"},AxesStyle->{Arrowheads@.03,Arrowheads@.03},GridLines->Automatic,PlotStyle->{Darker[Green,0.4],Thickness[0.008]},Ticks-> {Automatic,Table[{1000*i,i},{i,0,14,2}]}];
GraphicsGrid[{{graf1r,graf4r},{graf2r,graf3r}},Frame->All,ImageSize->1000]

tmax=mez*dt;
Manipulate[ListPlot[Table[{j*dt,m[j]},{j,0,time/dt}],PlotStyle->Darker[Green],PlotRange-> {{0,1.1*tmax},{0,1.01*Max[Table[m[i],{i,0,tmax/dt}]]}},PlotLabel->m[Round[time/dt]],BaseStyle->{FontFamily->"Times",FontSize-> 14},AxesLabel->{"t/s","m/t"},AxesStyle->{Arrowheads@.03,Arrowheads@.03},GridLines->Automatic,Ticks-> {Automatic,Table[{1000*i,i},{i,0,14,2}]},Epilog->{PointSize[Large],Darker[Red,0.3],Point[{time,m[Round[time/dt]]}]}],{time,0,tmax,Appearance->"Open"}]

Buquoy’s Problem

Historical Background

Graf Georg von Buquoy, Czech aristocrat, mathematician and gifted inventor (1781-1851)

Bucquoy was the first who investigated mechanical systems with a varying mass

In 1812 he explicitly formulated the correct dynamical equations of motion for the case when the mass of the moving object is changing

In 1815 he presented his results at the Paris Academy of Sciences

Formulation of the Problem

The first explicit example of a system with a varying mass suggested by Buquoy in 1814 is the following:

Consider an ideally flexible fibre lying reeled on a horizontal plane. Determine its motion when a constant vertical force (directed upward) is exerted on the end of the fibre.

Denoting the position of the end of the fibre above the horizontal plane by x>0 we make the following natural simplifying assumptions:

the vertical gravitational field is homogenous

the fibre is thin and its linear density η is constant

the fibre reels off without friction at the origin of the x-axis during the upward motion

the fibre “smoothly disappears” at the origin of the x-axis during the motion downward, i.e., the part that has already landed on the plane does not move

experimental model: the vertical motion of a balloon with a heavy rope hanging down

The Equation of Motion

The motion of the fibre is determined by Newton’s law:

(OverDot[p])=F-m g

p ... the momentum of the moving part of the fibre

F>0 ... vertically upward oriented constant force

m g ... the oppositely exerting weight of the reeling fibre

Supposing the mass is given by m =η x where η is a constant linear density of the fibre

(((x(OverDot[x])))^·)=(F/η)-g x

(Overscript[x, ¨])=(F/(η x))-g-(((OverDot[x])^2)/x)

A stacionary solution x =const., i.e. (OverDot[x])=(Overscript[x  , ¨])=0, represents the equilibrium situation such that the end of the fibre is situated at a constant height

(Subscript[x, c])=(F/(η g))

The problem thus seems to be reduced to finding the solution of equation (3). However, this equation only describes the upward motion of the fibre - it is only correct for a growing x, i.e. (OverDot[x])>0

the reason is that only in such a case the resultant force F-m g is responsible for the rate of change of the total momentum of the fibre

(OverDot[p])=m (Overscript[x, ¨])+(OverDot[m]) (OverDot[x])

In the case of a downward oriented motion ((OverDot[x])<0, i.e. decreasing x) the situation is different:

(Overscript[x, ¨])=(F/(η x))-g

Upward Motion

In the case (OverDot[x])>0, integrating the equation (2)

(1/2)(((x(OverDot[x])))^2)=(F/(2η))(x^2)-(g/3) (x^3)+C,

C is an integration constant dependent on initial conditions  (Subscript[x, 0]), (Subscript[(OverDot[x]), 0]):

C=(1/2)((((Subscript[x, 0])(Subscript[(OverDot[x]), 0])))^2)-(F/(2η))((Subscript[x, 0])^2)+(g/3) ((Subscript[x, 0])^3)

Equation (7) can be rewritten in the form

(1/2)((OverDot[x])^2)+(Subscript[V, ef]^↑)(x)=(F/(2η)),      (Subscript[V, ef]^↑)(x)=(g/3) x-(C/(x^2))

(Subscript[V, ef]^↑)(x) ... effective potential for the case of upward motion

Substituting (9) to (3) we obtain the equation

(Overscript[x, ¨])=-(g/3)-((2C)/(x^3))=(Subscript[F, ef])(x)

(Subscript[F, ef])(x) ... effective force

Manipulate[
xm=If[C1<0,x/.FindRoot[(10/3)x-(C1/(x^2))==20,{x,0.1}],-1];
xM=x/.FindRoot[(10/3)x-(C1/(x^2))==20,{x,6}];
Plot[{(10/3)x-(C1/(x^2)),20},{x,0,10},AxesLabel->{x,(Subscript[V, ef])},PlotRange->{-40,40},Ticks->{{{xm,(Subscript[x, m])},{xM,(Subscript[x, M])}},{{20,(F/(2η))}}}],
{{C1,0,"C"},-5,5}]

Conclusions:

For C⩾0 we observe that (Subscript[F, ef])<0 and the initial velocity must necesarily be (Subscript[(OverDot[x]), 0])>0

The initial condition (Subscript[(OverDot[x]), 0])=0 is only possible when (Subscript[F, ef])⩾0

For any value of the exerted force F there exists a maximal height (Subscript[x, M])

the solution (Subscript[x, m]) has the meaning of the starting point for the upward motion with zero initial velocity

Downward Motion

In the case (OverDot[x])<0, integrating the equation (6)

(1/2)((OverDot[x])^2)+(Subscript[V, ef]^↓)(x)=(Subscript[V, ef]^↓)((Subscript[x, 1])),      (Subscript[V, ef]^↓)(x)=g x-(F/η)log x

x2=FindRoot[10x-55Log[x]==-20,{x,1}];
x1=FindRoot[10x-55Log[x]==-20,{x,12}];
xc=FindMinimum[10x-55Log[x],{x,5}];
Plot[{10x-55Log[x],-20},{x,0,20},AxesLabel->{x,(Subscript[V, ef])},Ticks->{{{x/.x2,(Subscript[x, 2])},{x/.x1,(Subscript[x, 1])},{x/.xc[[2]],(Subscript[x, c])}},None}]

Conclusions:

For every potential level higher than (Subscript[V, ef]^↓)((Subscript[x, c])) there exist two turning points such that 0<(Subscript[x, 2])<(Subscript[x, c])<(Subscript[x, 1])

The coordinates of these points fulfill the condition

(Subscript[x, 1])-(Subscript[x, 2])=(Subscript[x, c])log((Subscript[x, 1])/(Subscript[x, 2]))

Necesarily  (Subscript[x, 2])>0, therefore the upper end of the fibre can never fall down back on the plane

The Complete Motion

If the downward motion follows the upward motion, the turning point (Subscript[x, 1]) is identical to the maximal value  (Subscript[x, M]) introduced in the case of upward motion

Similarly, the turning point (Subscript[x, 2]) becomes a new minimal value (Subscript[x, m]) for a subsequent upward motion

To describe the complete motion of the fibre we thus have to alternaly join solutions of equation (9) valid for (OverDot[x])>0 with solutions of (11) for (OverDot[x])<0

In our model, the equations of motion (3) and (6) are rewritten in unified form

(Overscript[x, ¨])=g(((Subscript[x, c])/x)-1)-(1/2)(1+sgn(OverDot[x]))(((OverDot[x])^2)/x)

η=5*^-3; (* lineární hustota *)
F=0.1; (*síla působící svisle vzůru *)
g=9.81;
(Subscript[y, c])=F/η/g; (* rovnovážná poloha *)
(Subscript[y, p])= 1; (* počáteční poloha *)
(Subscript[t, m]) = 20; (* maximální čas, do kterého bude probíhat numerická integrace *)
Konst = F/2/η*(Subscript[y, p])+g/3*(Subscript[y, p])^3;

Manipulate[Module[{sol=NDSolve[{y''[t]== g*((Subscript[y, c])/y[t]-1)-0.5*(1+Sign[y'[t]])*y'[t]^2/y[t],y[0]==p[[1]],y'[0]==p[[2]]},y,{t,0,T}]},Plot[{Evaluate[y[t]/.sol],(Subscript[y, c])},{t,0,T},PlotRange->{Automatic,{0,3.5}},Filling->{1->{2}},BaseStyle->{FontFamily->"Times",FontSize-> 14},AxesLabel->{"t","x"},GridLines->Automatic,PlotStyle->{Thickness[0.005]},PlotLabel->Evaluate[y[t]/.sol/.t->T][[1]],Epilog->{PointSize[Large],Darker[Red,0.3],Point[{T,Evaluate[y[t]/.sol/.t->T][[1]]}]}]],
Grid[{
{Control[{{p,{1,0}},Slider2D}],Dynamic[Style["x[0] = "<>ToString[p[[1]]],15]]},
{SpanFromAbove,Dynamic[Style["(OverDot[x])[0] = "<>ToString[p[[2]]],15]]}
},Spacings->{2,1}],{{T,.01},0.00001,40,Appearance->"Labeled"}]

THANK YOU FOR YOUR ATTENTION