How can I simulate this toggle mechanism?

Question

I've been introduced to Mathematica very recently. Basically, I haven't actually "solved" anything in my Mathematica lifetime, but I've done some simulations. With my fractional knowledge, I tried simulating a toggle mechanism today. Here's what I did...

ω0 = .02; ω = .02; r = 2; r1 = 6; r0 = 1.5;
rA[θ_, r_] := r*{Cos[θ], Sin[θ]};
rAy[θ_, r_] := r*Sin[θ];
rAx[θ_, r_] := r*Cos[θ];
rC[θ0_, r0_] := r0*{Cos[θ0] + 1/r0, Sin[θ0] + 3/r0};
rCy[θ0_, r0_] := r0*Sin[θ0] + 3;
rCx[θ0_, r0_] := r0*Cos[θ0] + 1;
xB[θ_, r_] := r*Cos[θ] + Sqrt[(r1^2 - (r*Sin[θ])^2)];
rB[θ_, r_] := {xB[θ, r], 0};
rBx[θ_, r_] := xB[θ, r];
rBy[θ_, r_] := 0;
Manipulate[θ = ω*t; θ0 = ω0*t; 
 rAO = rA[θ, r]; rBO = rB[θ, r]; rDC = rC[θ0, r0];
 Graphics[{{Thick, Darker[Green], Dashed, Circle[{0, 0}, r]}, {Thick, 
    Dashed, Blue, Circle[{1, 3}, r0]},
   {Red, EdgeForm[{Thick, Black}], White, 
    Rectangle[{rBO[[1]] - 0.5, rBO[[2]] - 0.5}, {rBO[[1]] + 1, 
      rBO[[2]] + 0.5}]},
   {Thick, Darker[Red], Line[{{0, 0}, rAO, rBO}]}, {Thick, 
    Darker[Red], Line[{{1, 3}, rDC, rAO}]},
   {Black, Disk[{0, 0}, 0.1], Disk[{1, 3}, 0.1], 
    Disk[{rAO[[1]], rAO[[2]]}, 0.1], Disk[{rDC[[1]], rDC[[2]]}, 0.1], 
    Disk[{rBO[[1]], rBO[[2]]}, 0.1]}},
  PlotRange -> {{-3, 9}, {-3, 5}}, 
  ImageSize -> {600, 250}], {{t, 20, "Motion"}, 20, N[205*Pi], 
  N[Pi/16]}]

I think you can figure out that there's a big issue, due to my enormous expectation from the "visualizing" without even thinking of "solving"...

Both the cranks have the same angular velocity. Because, that's the way I've defined. I've got no other choice. But, the mechanism should begin with the top crank. That crank makes the bottom one to rotate, by transferring power through the connecting rod.

Given a connecting rod of constant length (especially, less than the sum of radii of the cranks), the bottom crank cannot complete one full rotation. It can go only halfway, and then turn back. The slider motion remains the same, but the crank is causing a big trouble to the mechanism due to the ridiculous "elastic behavior" of the connecting rod.

I don't want to analyze this mechanism. So, I'm cool with any answer that essentially simulates the actual mechanism (even if that's a modification of mine) as I just need the simulation only for a video (at least, for now). At the same time, I'll be very happy with "solving" the equations with constraints, because I guess those solutions may help me in the future for other mechanisms.

Answer 1

One way is to set up a DAE: See tutorial/DSolveExamplesOfDAEs and example/ModelConstrainedSystemsAsDAEs.

The constraint that the driver (bottom rotating link) has a fixed length is taken care of by initial conditions and the DE. There are two possible starting positions for the driven link. One might have to inspect the result of Solve to determine which is desired.

constraints = ComplexExpand@{
    (* SquaredEuclideanDistance[{a1[t], a2[t]}, {0, 0}] == r1^2, not needed *)
    SquaredEuclideanDistance[{b1[t], b2[t]}, {a1[t], a2[t]}] == rod^2,
    SquaredEuclideanDistance[{c1, c2}, {b1[t], b2[t]}] == r2^2};

eqn = {
      (* DE to rotate small link *)
    D[{a1[t], a2[t]}, t] == omega/r1 Cross[{a1[t], a2[t]}],
      (* initial conditions *)
    {a1[0], a2[0]} == {0, r1},
    {b1[0], b2[0]} == {b1[t], b2[t]} /.
      First @ Solve[constraints /. Thread[{a1[t], a2[t]} -> {0, r1}], {b1[t], b2[t]}]};

Example

Block[{omega = 1, r1 = 1, rod = 3, r2 = 2, c1 = 1, c2 = 3},
  sol = First @ NDSolve[{eqn, constraints},
     {a1[t], a2[t], b1[t], b2[t]},
     {t, 0, 2 Pi/omega}]];

With[{c1 = 1, c2 = 3},
 Manipulate[
  Graphics[
   {Line[{{0, 0}, {a1[t], a2[t]}, {b1[t], b2[t]}, {c1, c2}} /. sol /. t -> t0]},
   PlotRange -> {{-2, 6}, {-2, 6}}
   ],
  {t0, 0, Dynamic@sol[[1, 2, 0, 1, 1, 2]]}]
 ]

[The piston can be included via trigonometry or by adding another constraint.]

---

A random linkage:

Module[{sol, x, y, t},
 DynamicModule[{wheels = 6, omega = 1, joints, centers, constraints, 
   variables, initial, eqn, t0 = 0., direction = Forward},

  joints = 
   Sort@RandomReal[1, {wheels, 2}] + Table[{i, 0}, {i, wheels}];
  centers = 
   Sort@RandomReal[1, {wheels, 2}] + Table[{i, (-1)^i}, {i, wheels}];
  variables = 
   Transpose@{Array[x[#][t] &, wheels], Array[y[#][t] &, wheels]};
  constraints = ComplexExpand@{
     MapThread[
      SquaredEuclideanDistance[#1, #2] == 
        SquaredEuclideanDistance[#3, #4] &,
      {Most[variables], Rest[variables], Most[joints], Rest[joints]}
      ],
     MapThread[
      SquaredEuclideanDistance[#1, #3] == 
        SquaredEuclideanDistance[#2, #3] &,
      Rest /@ {variables, joints, centers}
      ]
     };
  initial = Thread /@ Thread[variables == joints /. t -> 0] // Flatten;
  eqn = {D[First@variables, t] == 
     omega/SquaredEuclideanDistance[First@joints, 
        First@centers] Cross[First@variables - First@centers]};

  Check[sol = 
    First@NDSolve[Flatten@{eqn, constraints, initial}, 
      Flatten[variables], {t, -2 Pi/omega, 2 Pi/omega}], 
   direction = ForwardBackward, {NDSolve::ndsz, NDSolve::ndcf}];

  Column[{
    Labeled[
     Manipulator[Dynamic[t0], First[x[1][t] /. sol /. t -> "Domain"], 
      AnimationDirection -> direction], "t", Left],
    Graphics[
     Dynamic@With[{pts = variables /. sol /. t -> t0},
       {{Red, Thin,
         MapThread[
          Circle, {centers, 
           MapThread[EuclideanDistance, {centers, joints}]}]},
        Line[pts], Line[Transpose[{centers, pts}]],
        EdgeForm[Black], Red, Disk[#, 0.1] & /@ pts
        }],
     PlotRange -> {{0, wheels + 2}, {-2, 2}}
     ]
    }]
  ]]

Answer 2

I did a solution with contour tracing on the distance function. It gets pretty unstable sometimes, but it's a fun question to experiment with interactivity.

DynamicModule[{p1 = {0, 2}, p2 = {1, 3}, angles = {0, 0}, distance, 
  grad, tangent}, 
 distance[a1_, a2_] := 
  Norm[{Cos@a1, Sin@a1} - (Norm[p2 - p1] {Cos@a2, Sin@a2} + p1)]; 
 grad = D[distance[a1, a2], {{a1, a2}}] /. Abs' -> Sign; 
 tangent = Cross@grad; 
 Column[{Dynamic[
    angles = 
     Mod[# /. FindRoot[distance @@ # == distance @@ angles, {t, 0}] &[
       grad t + angles + .05 tangent /. Thread[{a1, a2} -> angles]], 
      2 Pi]; Graphics[{{Dashed, Circle[], Circle[p1, Norm[p2 - p1]], 
       Locator@Dynamic@p1, Locator@Dynamic@p2}, Red, 
      Line[{{0, 0}, {Cos@#, Sin@#}, 
          p1 + Norm[p2 - p1] {Cos@#2, Sin@#2}, p1} & @@ angles]}, 
     PlotRange -> 5, Frame -> True]], 
   Dynamic@LocatorPane[Dynamic@angles, 
     ContourPlot[distance[a1, a2], {a1, 0, 2 Pi}, {a2, 0, 2 Pi}, 
      ContourLabels -> True]]}]]

Answer 3

You can do this without a NDSolve by calculating the distance from the follower cranks joint to the end of the driving cranks end. Then use this distance with law of cosines to calculate the deviation angle. This is also pretty easy to implement on ANY hardware capable of doing a ArcCos and Atan2 operation (note that in c atan2 parameters are swapped).

With[{p0 = {0, 0}, p3 = {1, 3}, r1 = 1, r2 = 2, r3 = 3},
 Module[{p1, p2, s, d, a1, a2},
  Animate[
   p1 = r1 {Cos[a], Sin[a]};
   s = p1 - p3;
   d = Norm[s]; 
   a1 = ArcTan[s[[1]], s[[2]]];
   a2 = a1 + ArcCos[( d^2 + r2^2 - r3^2)/(2 d r2)];
   p2 = p3 + r2 {Cos[a2], Sin[a2]};
   Graphics[
    Line[{p0, p1, p2, p3}],
    PlotRange -> {{-2, 6}, {-2, 6}}
    ],
   {a, 0, 2 Pi}]
  ]]

Results in the same thing as the accepted answer:

This solution is much much simpler, but won't work in 3D cases. The accepted answer will, and its a more flexible model. Finally, the position of the piston, which wasn't calculated by other posters as trivial is:

pp = {r1*Cos[a] + r4 + Sin[a]^2/(2 r4), 0};

And you get the following by fiddling a bit with the graphics elements: