This is exactly what I suggested David. here is an attempt to "glue" the two together using the trebuchet NDSolve output as the initial conditions of the projectile equations. I attach a preliminary notebook based on my previous demonstration. Still needs some refining!

Attachments:

Perhaps you could make a use of System Modeler too? Take a look at Model Your Own Medieval Catapult Design.

SystemModeler is indeed a better choice for modelling this type of "machinery"

But a real trebuchet is more complicated than is depicted in all the above (including my own demonstrations!) If you look at the paper "Trebuchet Mechanics" by Donald B. Siano and at this design from the paper:

you can see that there is a third link in the kinematic chain that was neg;acted in the above: the part indicated l3 (between {x2,y2} and {x3,y3} laying horizontally on the ground where the projectile is attached. Adding this link complicates the analysis (Siano does it at the end!) but it is essential for the trebuchet deadly efficiency since it will throw the load much further and with a higher speed! One day, when more time available, I will attempt to add tis feature but it would be a nice Mathematica/SystemModeler combined exercise.

Any volunteers?

Hi, Thanks

I was able to generate the equations of motions for the Trebuchet and calculate the Lagrangian.

I tried using the equations of Projectile motion as well. But I am not quite sure as to how to combine the equations from Projectile with Air Drag into the Trebuchet equations to animate the projectile to actually leave the Trebuchet and make an impact against something like a wall at a distance say x='constant'.

Is there any specific way to do that?

The equations I got were:

cn = {L1, L2, L3, L4, m1, m2, mb, rg, rc};


(*Coordinates of the short end of the beam*)     (*Defined 'a'*)
x1[?_] := L1 *Sin[?];
y1[?_] := -L1 *Cos[?];
(*Coordinates of the long end of the beam:*)    (*Defined 'b'*)
x2[?_] := -L2* Sin[?];
y2[?_] := L2 *Cos[?];
(*Coordinates of the center of mass of the Counterweight*)     \
(*Defined 'c'*)
x4[?_, ?_] := 
  L1* Sin[?] - L4 *Sin[? + ?];
y4[?_, ?_] := -L1 *Cos[?] + 
   L4* Cos[? + ?];
(*Coordinates of the Projectile end of the Sling *)    (*Defined 'd'*)


x3[?_, ?_] := -(L3 *Sin[? - ?] + 
     L2* Sin[?]);
y3[?_, ?_] := -(L3 *Cos[? - ?] - 
     L2* Cos[?]);


(* potential energy of the three parts:*)
PE[?_, ?_, ?_] := 
  m1 *g* y4[?, ?] + m2 *g* y3[?, ?] - 
   mb* g* rc *Cos[?];

(* kinetic energy of the three parts:*)
KE[?_, ?_, ?_] := 
  0.5*m1*((Dt[x4[?, ?], t, Constants -> cn])^2 + (Dt[
         y4[?, ?], t, Constants -> cn])^2
     ) + 0.5*
    m2 *((Dt[x3[?, ?], t, Constants -> cn])^2 + (Dt[
         y3[?, ?, Constants -> cn], t])^2) + 
   (mb/2) *(rg^2)* Dt[?, t, Constants -> cn]^2;
(* the lagrangian*)
Lag = KE - PE;
L[?_, ?_, ?_] := 
  KE[?, ?, ?] - PE[?, ?, ?];


(* define what is constant, and the derivatives of the angles, \
denoted by the "d" suffix:*)
EL = L[?, ?, ?] /. {Dt[?, t, 
      Constants -> {L1, L2, L3, L4, m1, m2, mb, rg, rc}] -> ?d,
        Dt[?, t, 
      Constants -> {L1, L2, L3, L4, m1, m2, mb, rg, rc}] -> ?d,
           Dt[?, t, 
      Constants -> {L1, L2, L3, L4, m1, m2, mb, rg, rc}] -> ?d, 
    Dt[?, t] -> ?d, Dt[?, t] -> ?d, 
    Dt[?, t] -> ?d};



eq? = 
  Simplify[ 
   Dt[D[EL, ?d], t] - D[EL, ?] == 0 /. {Dt[L1, t] -> 0, 
     Dt[L2, t] -> 0, Dt[L3, t] -> 0, Dt[L4, t] -> 0, Dt[mb, t] -> 0, 
     Dt[m1, t] -> 0, Dt[m2, t] -> 0, Dt[g, t] -> 0, Dt[rg, t] -> 0, 
     Dt[rc, t] -> 0,
          Dt[?, t] -> ?d, Dt[?, t] -> ?d, 
     Dt[?, t] -> ?d, Dt[?d, t] -> ?dd, 
     Dt[?d, t] -> ?dd, Dt[?d, t] -> ?dd}];

eq? = 
  Simplify[Dt[D[EL, ?d], t] - D[EL, ?] == 
     0 /. {Dt[L1, t] -> 0, Dt[L2, t] -> 0, Dt[L3, t] -> 0, 
     Dt[L4, t] -> 0, Dt[mb, t] -> 0, Dt[m1, t] -> 0, Dt[m2, t] -> 0, 
     Dt[g, t] -> 0, Dt[rg, t] -> 0, Dt[rc, t] -> 0,
          Dt[?, t] -> ?d, Dt[?, t] -> ?d, 
     Dt[?, t] -> ?d, Dt[?d, t] -> ?dd, 
     Dt[?d, t] -> ?dd, Dt[?d, t] -> ?dd}];

eq? = 
  Simplify[Dt[D[EL, ?d], t] - D[EL, ?] == 
     0 /. {Dt[L1, t] -> 0, Dt[L2, t] -> 0, Dt[L3, t] -> 0, 
     Dt[L4, t] -> 0, Dt[mb, t] -> 0, Dt[m1, t] -> 0, Dt[m2, t] -> 0, 
     Dt[g, t] -> 0, Dt[rg, t] -> 0, Dt[rc, t] -> 0,
          Dt[?, t] -> ?d, Dt[?, t] -> ?d, 
     Dt[?, t] -> ?d, Dt[?d, t] -> ?dd, 
     Dt[?d, t] -> ?dd, Dt[?, t] -> ?dd}];


sol = Solve[{eq?, eq?, eq?}, {?dd, ?dd, ?dd}];

Thanks

Hi Varun,

I agree with Erik on this. The trebuchet part of the problem determines the coordinates and velocity at the point where the projectile leaves the pouch. These are the initial conditions for a second set of differential equations which describe ballistic flight. There are two ways to string these together: You can do so manually by executing them in sequence, either in separate cells or in the same cell; or they can be placed in a module, with local variables connecting them, and both solutions used for building the plot. To animate them together, they would need to be all encapsulated in the animation.

Using the impact point as an input is another matter. It is really part of the solution. The likely method here is to build the full simulation into a function which can return the impact point. A solver like FindRoot can then be used to solve for some selected input parameter being used to adjust the impact. In a simple projectile problem that might be the launch angle. This gets you to a way to specify the impact point as an input, and retrieve the parameter setting needed, as well as the full trajectory.

Before version 10, this would require that NDSolve be enclosed in the FindRoot process, which means solving the set of diffeqs many times. However, with V10 we have ParametricNDSolve, which allows us to specify an input parameter and returns a function for a solution to the diffeqs, given a value of the input parameter. This output from ParametricNDSolve can be used as an alternative to wrapping the solution process in the solve.

Best regards, David

Hello, Frederick, I'm glad you like my trebuchet demos. The equations used in the code itself are always the same and the one in the text "should" be the same but there might be some differences due to converting Normal style into Traditional? Anyway, the ones in the code are the correct ones.

As to efficiency, yes this is a purely theoretical model. Adding friction, air drag,...losses is not as simple as adding a factor of some %. In the procedure to get the Lagrangian equations (see attached notebooks), you would have to add additional energy factors in the formulas for kinetic- and potential energy. Air drag is proportional to some power of speed etc...Adding this into the energies will change the derivatives etc. and add a lot of complexity to the equations. But it could be done. In a more academic study, this would be required but for me it is beyond my capabilities. You get 3 equations of motion because there are 3 degrees of freedom: the three angular variables as functions of time: theta, psi and phi.

Extracting the 3 equations of motion out of the Lagrangian is explained in the attached "lagrangiansnb".

It is based on a Blog by Andrew Moylan of Wolfram. See "Stabilised..." attached

Attachments:

Hello,

I tried solving both equations separately, but I face an error everytime saying that the limiting value is not a machine sized number and a certain replacement rule is invalid. I cannot quite figure it where that error exists.

I used the below mentioned Projectile equations as suggested for the Trebuchet, but, cannot get the projectile to leave the Trebuchet.

Manipulate[
 ParametricPlot[{{V Cos[\[Theta]] t, 
    V Sin[\[Theta]] t - 
     9.8 t^2/2}, {V vt Cos[\[Theta]]/9.8 (1 - Exp[-9.8 t/vt]), 
    vt/9.8  (V Sin[\[Theta]] + vt) (1 - Exp[-9.8 t/vt]) - vt t}}, {t, 
   0, tf}, PlotStyle -> {{Thick, Blue}, {Thick, Red}}, 
  AxesOrigin -> {0, 0}, AxesLabel -> {x, z}, 
  PlotRange -> {{0, 400}, {0, 300}}, ImageSize -> 500, 
  Epilog -> {Green, PointSize[.02], 
    Point[{{V Cos[\[Theta]] tf, 
       V Sin[\[Theta]] tf - 
        9.8 tf^2/2}, {V vt Cos[\[Theta]]/9.8 (1 - Exp[-9.8 tf/vt]), 
       vt/9.8  (V Sin[\[Theta]] + vt) (1 - Exp[-9.8 tf/vt]) - 
        vt tf}}]}], {{V, 48.8, "initial velocity"}, 1, 100, 
  Appearance -> "Labeled"}, {{\[Theta], 1.09, "angle"}, .1, \[Pi]/2, 
  Appearance -> "Labeled"}, {{vt, 100, "terminal velocity"}, .01, 300,
   Appearance -> "Labeled"}, {{tf, 8.75, "time"}, .1, 17, 
  Appearance -> "Labeled"}]

How do I separately solve for the Trebuchet and the Projectile, so that I can get an animation where the projectile actually leaves the trebuchet?

Thanks.

The demo is beautiful, Erik. Thank you. I love Lagrangian mechanics, and here we have additional interest in knowing these things brought the age of castles to an end!

Hello David,

Apologies for the delayed response.

Thank you for explaining all the concepts in detail and providing a well defined approach for problem solving and modelling.

It really helped in the understanding of different aspects and the way equations are used and solved in Mathematica.

Thanks,

Varun

Hello Erik,

Thank you very much for this helpping and guide in equations.

I am actually holding a trebuchet competition for young and junior engineers in company. It's start from building a small wood trebuchet, then predict a shooting range, then test and verification. I try to use your demo to predict shooting range. Maybe, modified a little to fit practical condition, I don't know how, so far. The small wood toy is lovely, but it after tested, It's always a low efficiency (about 10% only) to convert potential energy to kinematic. So we are struggle to improve the toy efficiency a little bit higher in the meantime.

The trebuchet competition will be hold in this weekend. So far, I am wondering, if any of our young engineers would make something new and useful. I would keep you informed in more details, as soon as we find something interesting.

@Erik Mahieu,

If I understand you right. You mean the 1st phase of Trebuchet, can be seperated two (1.1 and 1.2 ) in detail. http://demonstrations.wolfram.com/TrebuchetBallistics/

phase 1:

1.0 Start from all stationary state

1.1 Only two masses ( counterweight and main beam ) involved in physical with 2 DOF rotation. Becasue, sling or rope at projectile side is not straight and no force. projectile is still on ground in this phase.

1.2 Three masses ( counterweight, projectile and main beam ) involved in physical with 3 DOF rotation, sling is straight with a force to pull projectile until it released.

phase 2: the projectile is free throw with ballistic trajectory in air, with air drag. ( Air drag is clear.)

I think below factors in phase 1 might need to consider for efficiency.

1. mass of main beam as a rigid bar;

2. damping or friction on three axes;

3. collision at relearing or after releasing (or the remaining energy of the system, where is it ? ) ;

I think 3rd part of energy is consumed by collision when the beam hit something as end-stop. Usually we can use the position of end-stop to control the velocity direction of projectile. See Siano's page: http://algobeautytreb.com/

It appears Siano uses a bar as end-stop to control the velocity direction of projectile. The moment, when main bean hit on a bar, should be the moment the projectile releasing. Maybe you have conside the 3 part in side equations ?

@Erik, Anthor finding, after careful looking at both on your horse image and my da Leonard's drawing. It can be those ancestors, open a groove and use a slot underneath the ground, so that sling at projectile side can have more bigger effective angle.That's a simple way to break that geometry constrain limited by ground.

So that l3 neglect should be acceptable, only if the man cut a slot on ground.

David said >

It is easier to find your way through a problem if you put each calculation in its own cell, and you evaluate and look at the result as you go along.

I had thought Mathematica 11.0 has a feature that debugs functions step by step, it has breakpoints i cannot get it to work yet. Documentation is very thin - despite setting breakpoints it never stops :)

I made the below for an earlier Mathematica (4.0). I'm unsure if it works "for complicated functions" in Mathematica 11.0 or not. It is in jutilspkg part of fNBookForm2 (orig. part of PeriodicTable` pkg). It is supposed to do what you said without all the ; editing. This allows being able to freeform work a cell at a time, then have a function made auto. It works but needs tender care, and has a help page.

Attachments:

Hi Vitaliy,

Erik's Demonstrations help us to understand the physical model and trebuchet dynamics. As a organizer, I draft a Mathematica notebook to introduce trebuchet topic for our young players. I attached notebook format, see attachment. (Sorry, only Chinese Version).

Unfortunately, Within the limited time and energy, we didn't develpe model or implement on code beyond Erik's Demonstrations. Anyway, we learned a lot through competition in this year. We would make it even better in next year, definitely.

Attachments: