15
$\begingroup$

Torque-free Euler equations experiment seen in low gravity of Russian spacecraft is modelled here with a view to see its tumbling motion around the intermediate axis $\omega_2$ rotation. However its reversal is not observed here. The boundary conditions do play a role, varying them did not much change the sine behaviour towards interfering periodic flips.

Due to easy demonstration possibility here I posted this hopefully interesting problem although strictly speaking it is a physics problem.

{I1, I2, I3} = {8, 4, 0.4};
Dzhanibekov = {I1 TH1''[t] == (I2 - I3) TH2'[t] TH3'[t], 
   I2 TH2''[t] == (I3 - I1) TH3'[t] TH1'[t], 
   I3 TH3''[t] == (I1 - I2) TH1'[t] TH2'[t], TH1'[0] == -0.4, 
   TH2'[0] == 0.08, TH3'[0] == 0.65, TH1[0] == 0.75, TH2[0] == -0.85, 
   TH3[0] == 0.2};
NDSolve[Dzhanibekov, {TH1, TH2, TH3}, {t, 0, 15.}];
{th1[u_], th2[u_], th3[u_]} = {TH1[u], TH2[u], TH3[u]} /. First[%];
Plot[Tooltip[{th1'[t], th2'[t], th3'[t]}], {t, 0, 15}, 
 GridLines -> Automatic]

Please help choose better initial conditions for getting a jump around $\theta_2$ axis. Thanks in advance.

EDIT1:

ICs updated per MichaelE2's suggestion showing the effect on the three angular velocity variations.

The flip frequency is surprisingly dependent on choice of ICs. Is it possible to determine the common frequency analytically?

Wing Nut Flips

Wiki Ref

enter image description here

Improve this question
$\endgroup$
3
  • 1
    $\begingroup$ Use ICs from the physics of the phenomenon you wish to illustrate: TH1'[0] == -0.01, TH2'[0] == 1, TH3'[0] == 0.01, TH1[0] == 0.0, TH2[0] == -0.0, TH3[0] == 0.0 $\endgroup$
    – Michael E2
    Commented Oct 17, 2020 at 0:51
  • $\begingroup$ Thanks soo much. But what beats me is the phenomenon. Time period depending upon the IC !! Did n't read much into the referred articles though. You can edit the code if so you wish to.. $\endgroup$
    – Narasimham
    Commented Oct 17, 2020 at 5:43
  • $\begingroup$ You would expect some dependency to be predictable, such as halving the period if you double the TH derivatives, no? $\endgroup$
    – Michael E2
    Commented Oct 17, 2020 at 18:08

2 Answers 2

Reset to default
24
$\begingroup$

If this is a physical problem then choice of I1,I2,I3 depends on the form of the body we are tested. To make animation we first make a body as, for example,

Graphics3D[{Cone[{{0, 0, 0}, {0, 0, 3}}, 1/2], 
  Cuboid[{-0.2, -1, 0}, {0.2, 1, .7}]}, Boxed -> False]


G3D = RegionUnion[Cone[{{0, 0, 0}, {0, 0, 3}}, 1/2], 
   Cuboid[{-0.3, -1, 0}, {0.3, 1, 1}]];

c = RegionCentroid[G3D];

Figure1
Then we calculate moment of inertia and define equations

J3 = NIntegrate[x^2 + y^2, {x, y, z} \[Element] G3D];
J2 = NIntegrate[x^2 + (z - c[[3]])^2, {x, y, z} \[Element] G3D];
J1 = NIntegrate[y^2 + (z - c[[3]])^2, {x, y, z} \[Element] G3D];



 eq1 = {\[CapitalOmega]1[
     t] == \[CurlyPhi]'[t]*Sin[\[Theta][t]]*
      Sin[\[Psi][t]] + \[Theta]'[t]*Cos[\[Psi][t]], \[CapitalOmega]2[
     t] == \[CurlyPhi]'[t]*Sin[\[Theta][t]]*
      Cos[\[Psi][t]] - \[Theta]'[t]*Sin[\[Psi][t]], \[CapitalOmega]3[
     t] == \[CurlyPhi]'[t]*Cos[\[Theta][t]] + \[Psi]'[t]};
eq2 = {J1*\[CapitalOmega]1'[t] + (J3 - J2)*\[CapitalOmega]2[
       t]*\[CapitalOmega]3[t] == 0, 
   J2*\[CapitalOmega]2'[t] + (J1 - J3)*\[CapitalOmega]1[
       t]*\[CapitalOmega]3[t] == 0, 
   J1*\[CapitalOmega]3'[t] + (J2 - J1)*\[CapitalOmega]2[
       t]*\[CapitalOmega]1[t] == 0};
eq3 = {\[CurlyPhi][0] == .001, \[Theta][0] == .001, \[Psi][
     0] == .001, \[CapitalOmega]3[0] == 
    10, \[CapitalOmega]1[0] == .0, \[CapitalOmega]2[0] == .025};

Finally we export gif file

Export["C:\\Users\\...\\Desktop\\J0.gif", 
 Table[Graphics3D[{Cuboid[{5, 5, -3}, {5.2, 5.2, 5}], 
    Cuboid[{-5, -5, -3.1}, {5, 5, -3}], 
    GeometricTransformation[{Cone[{{0, 0, 0}, {0, 0, 3}}, 1/2], 
      Cuboid[{-0.2, -1, 0}, {0.2, 1, .7}]}, 
     EulerMatrix[{NDSolveValue[{eq1, eq2, eq3}, \[CurlyPhi][tn], {t, 
         0, tn}], 
       NDSolveValue[{eq1, eq2, eq3}, \[Theta][tn], {t, 0, tn}], 
       NDSolveValue[{eq1, eq2, eq3}, \[Psi][tn], {t, 0, tn}]}]]}, 
   Boxed -> False, Lighting -> {{"Point", Yellow, {10, 3, 3}}}], {tn, 
   0, 11.6, .1}], AnimationRepetitions -> Infinity]

Figure 2

This problem has an analytical solution explained by Landau L.D., Lifshits E.M. in Mechanics. Let put $E$ is energy, $M^2$ is a squared angular momentum, $I_1,I_2,I_3$ are principal moments of inertia, $k^2=\frac{(I_2-I_1)(2EI_3-M^2)}{(I_3-I_2)(M^2-2EI_1)}$ , $sn(\tau,k), cn(\tau,k), dn(tau,k)$ -are Jacobi elliptic functions. Then the solution of the problem can be written in a closed form as
$$\Omega_1=\sqrt {\frac{2EI_3-M^2}{I_1(I_3-I_1)}}cn(\tau,k)$$ $$\Omega_2=\sqrt {\frac{2EI_3-M^2}{I_2(I_3-I_2)}}sn(\tau,k)$$ $$\Omega_3=\sqrt {\frac{-2EI_1+M^2}{I_2(I_2-I_1)}}dn(\tau,k)$$ $$\tau=t\sqrt {\frac{(-2EI_1+M^2)(I_3-I_2)}{I_1I_2I_3}}sn(\tau,k)$$ The dynamics of system is determined by two parameters - the period $T$ and the time of the flip $T_f$, which are related to each other as $T=4K(k)\sqrt{\frac{I_1I_2I_3}{(I_3-I_2)(M^2-2EI_1)}}, T_f=\frac{T}{2K(k)} $ where $K(k)$ is complete elliptic integral of the first kind.

Improve this answer
$\endgroup$
6
  • $\begingroup$ Nice animation on Mathematica. Indeed a new kind of science. Can we get the tumbling frequency as a function of moments of inertia and ICs? $\endgroup$
    – Narasimham
    Commented Oct 17, 2020 at 16:12
  • $\begingroup$ This code dated of March 2016 therefore it created first with v.11. The dynamics of system is determined by two parameters - the period and the time of the flip, which are related to each other. If you interested in this theory I can add some equations to my answer. $\endgroup$
    – Alex Trounev
    Commented Oct 17, 2020 at 16:31
  • $\begingroup$ Sure, the analytic solution is interesting to me, and perhaps to others also. $\endgroup$
    – Narasimham
    Commented Oct 17, 2020 at 16:37
  • 1
    $\begingroup$ @Narasimham See update to my answer. $\endgroup$
    – Alex Trounev
    Commented Oct 17, 2020 at 17:16
  • $\begingroup$ Marvellous, thanks for the answer. I suppose time period $T$ of the rotor is much smaller than flip period $T_f.$ $\endgroup$
    – Narasimham
    Commented Oct 17, 2020 at 21:48
14
$\begingroup$

I wrote something up for class once, using the same physical problem, to show how the implicit midpoint method preserves quadratic invariants. The example also shows that the Euler method, by comparison, allows the invariants to drift. But instead of solving for the angles as in the OP, the system is set up in terms of angular momentum. I hesitated to post it at first, but then the OP expressed an interest in the variability of the "period." One can intuit the dependence on the initial conditions from the phase portrait, as there is a saddle point on the $L_2$ angular momentum axis that corresponds to the middle angular inertia magnitude. An orbit slows down as it approaches the saddle point, and the closer it passes to the saddle point, the longer its period will be.

(* for the specific example orbits *)
Ip = {2, 3, 4};
{I1, I2, I3} = Ip;
icSol = First@
   FindInstance[{L1^2/I1 + L2^2/I2 + L3^2/I3 == 5/16, 
     L1^2 + L2^2 + L3^2 == 1, L2 > 0, L3 > 0}, {L1, L2, L3}];
ic = {L1, L2, L3} /. icSol;
L0p = 1; (* used in Mesh of plot *)

(* phase portrait for |L| == 1 *)
plot = Show[
   ParametricPlot3D[{Sin[v] Sin[u], Cos[v], 
     Cos[u] Sin[v]}, {u, -π, π}, {v, 0, π},
    MeshFunctions -> {Function[{L1, L2, L3, u, v}, 
       L1^2/I1 + L2^2/I2 + L3^2/I3]},
    Mesh -> {Subdivide[L0p^2/I3, L0p^2/I1, I1 I2 I3/2]},
    MeshShading -> 
     ColorData["TemperatureMap"] /@ Subdivide[0., 1., I1 I2 I3/2],
    PlotPoints -> 75,
    AxesLabel -> HoldForm /@ {L1, L2, L3}, Axes -> False, 
    Boxed -> False, ViewPoint -> {3, 3, 3}],
   Graphics3D[{
     Thick, Line[{{0, 0, 0}, 1.3 #}] & /@ IdentityMatrix[3],
     MapThread[
      Text, {{Subscript[L, 1], Subscript[L, 2], Subscript[L, 3]}, 
       1.4 IdentityMatrix[3]}]
     }]
   ];

ImplicitMidpoint = (* define NDSolve method *)
  {"FixedStep", 
   Method -> {"ImplicitRungeKutta", 
     "Coefficients" -> "ImplicitRungeKuttaGaussCoefficients", 
     "DifferenceOrder" -> 2, 
     "ImplicitSolver" -> {"FixedPoint", 
       "AccuracyGoal" -> MachinePrecision, 
       "PrecisionGoal" -> MachinePrecision, 
       "IterationSafetyFactor" -> 1 }}};

sys = {ode, ics} = {
    {L'[t] == -Cross[L[t]/I0[t], L[t]]},
    {L[0] == ic, I0[0] == {I1, I2, I3}}};

(* particular solution:implicit midpoint method *)
{solMP} = NDSolve[
   sys, L, {t, 0, 100},
   Method -> ImplicitMidpoint, StartingStepSize -> 2, 
   DiscreteVariables -> {I0}];

(* particular solution: Euler method *)
{solE} = NDSolve[
   sys, L, {t, 0, 100},
   Method -> {"FixedStep", Method -> {"ExplicitEuler"}}, 
   StartingStepSize -> 1, DiscreteVariables -> {I0}];

(* a separatrix:
 *   from the intersections of the sphere |L| = 1 and the planes
 *   L[[1]] / Sqrt[I1 (I3 - I2)] == ±L[[3]] / Sqrt[I3 (I2 - I1)] *)
{solSep} = NDSolve[
   {ode,
    {L[0] == Normalize@{Sqrt[I1 (I3 - I2)], -450, Sqrt[I3 (I2 - I1)]},
     I0[0] == {I1, I2, I3}}},
   L, {t, 0, 100},
   Method -> ImplicitMidpoint, StartingStepSize -> 2, 
   DiscreteVariables -> {I0}];

Clear[I1, I2, I3];
Legended[
 Show[
  plot,
  paths = Graphics3D[{
     {Red, Sphere[ic, 0.03]},
     Thick,
     MapThread[{#2, Tube[L["ValuesOnGrid"] /. #1, 0.015], 
        Sphere[#, 0.02] & /@ (L["ValuesOnGrid"] /. #1)} &,
      {{solMP, solE, solSep}, 
       First["DefaultPlotStyle" /. (Method /. 
            Charting`ResolvePlotTheme[Automatic, 
             ParametricPlot3D])][[2 ;; 4]]}
      ]
     }],
  PlotRange -> All, ViewPoint -> {6, 5, 8}
  ],
 LineLegend[
  First["DefaultPlotStyle" /.
   (Method /. 
     Charting`ResolvePlotTheme[Automatic, ParametricPlot3D])][[2 ;; 4]],
  {"Implicit Midpoint", "Explicit Euler", "Separatrix"}]
 ]

enter image description here

DSolve can solve the system exactly. It returns several cases, some of which are periodic and correspond to closed, nontrivial orbits.

Clear[I1, I2, I3, L0, L1, L2, L3];
Block[{L, I0},
   L[t_] = {L1[t], L2[t], L3[t]};
   I0 = {I1, I2, I3};
   dsol = DSolve[{L'[t] == -(L[t]/I0) \[Cross] L[t]}, L[t], t]
   ]; // AbsoluteTiming

(*  {3.04307, Null}  *)

I got lucky solving by inspection for the constants of integration {C[1], C[2], C[3]} that give the same initial conditions as the numerical example I coded.

constants = {C[1] -> 3/8, C[2] -> 1/8, C[3] -> 0};
dsol[[-3]] /. Thread[{I1, I2, I3} -> Ip] /. constants /. t -> 0
Values@% == Values@icSol

(*
  {L1[0] -> 0, L2[0] -> Sqrt[3]/2, L3[0] -> 1/2}
  True
*)

From the general solution, which is in terms of Jacobi elliptic functions, we can get a symbolic expression for the period and check it against the numerical example.

period = DeleteDuplicates[
   Cases[dsol, _JacobiSN | _JacobiDN, Infinity], 
   Simplify@Abs[List @@ #1/List @@ #2] == {1, 1} &] /. {_[u_, m_]} :> 
   4/Coefficient[u, t] EllipticK[m]
period0 = % /. Thread[{I1, I2, I3} -> Ip] /. constants // RealAbs
N@%

(*     symbolic period
  (2 Sqrt[2] Sqrt[I1] Sqrt[I2] I3 EllipticK[((-I1 + I2) I3 C[1])/
    (I2 (I1 - I3) C[2])])/(Sqrt[I1 - I3] Sqrt[I2 - I3] Sqrt[C[2]])
*)(*   period for the numerical example
  32 Sqrt[3] EllipticK[-2]
  64.9267
*)

Plot[Table[Indexed[L[t], k] /. solMP, {k, 3}] // Evaluate,
 {t, 0, period0}, PlotLegends -> {L1, L2, L3}, 
 PlotStyle -> Table[ColorData[97][k], {k, 3}]]

enter image description here

One can observe separatrices in the phase portrait. These paths go from one $L_2$ pole to the other in infinite time.

Plot[Table[Indexed[L[t], k] /. solSep, {k, 3}] // Evaluate,
 {t, 0, 100}, PlotLegends -> {L1, L2, L3}, 
 PlotStyle -> Table[ColorData[97][k], {k, 3}]]

enter image description here

One can also integrate around an orbit to obtain the period, the time being the integral of $ds/v$, where $v=|d{\bf L}/dt|$ is the speed along the orbit. The code below integrates over $1/4$ of an orbit and multiplies the integral by $4$.

Block[{L, I0, I1, I2, I3, L0, E0, norm},
 norm = Sqrt[#.#] &;
 L = {L1, L2, L3};
 I0 = {I1, I2, I3} = Ip;
 {L0, E0} = {L, L/I0} /. Thread[L -> ic];
 E0 = L0.E0/2;
 L0 = norm[L0];
 Assuming[
  0 < I1 < I2 < I3 && E0 > 0 && L0 > 0 && L ∈ Reals && 
   And @@ Thread[L0 > Abs[L]] && And @@ Thread[Dt@L > 0] && L1 > 0 && 
   2 E0 I1 < L0^2 < 2 E0 I3,
  integrand =(* T <- integrate ds/speed *)
   norm[Dt@L]/norm[-(L/I0) \[Cross] L] /.
    Last@Normal@Solve[{eqL, eqE}, {L2, L3}, Reals]; 
  expr = 4 Integrate[integrand /. Dt[L1] -> 1, L1];
  (* apply FTC *)
  Limit[expr, ##] & @@@ {
      {L1 -> 0, Direction -> "FromAbove"},
      Append[
       Last@Normal@Solve[{eqL, eqE} /. L2 -> 0, {L1}, {L3}, Reals] /. 
        Power[b_, 1/2] :> Power[1/b, -1/2], 
       Direction -> "FromBelow"]} // Differences // First
  ]]

(*  32 Sqrt[3] EllipticK[-2]  *)
Improve this answer
$\endgroup$
3
  • 1
    $\begingroup$ Very nice answer (+1). It looks like some extension of the tennis rocket theorem with numerical cases. $\endgroup$
    – Alex Trounev
    Commented Oct 19, 2020 at 10:42
  • $\begingroup$ @Michael E2: Thanks for the instructive and very detailed numerical illustration. Reading up the Landau/ Lifshitz article might also help further perhaps. And theoretically at least if required $\omega $,energy/momentum are given, it is possible to calculate back the three moments of inertia... right? $\endgroup$
    – Narasimham
    Commented Oct 20, 2020 at 22:37
  • $\begingroup$ It is a great answer. $\endgroup$
    – A little mouse on the pampas
    Commented Nov 6, 2020 at 5:00

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.