Hi, I'm new here, but I'm a friend of btrettel. I'm also a mechanical engineer, and, coincidentally, I do a lot of work on homemade spring powered systems like these. While we tend to shoot homemade Nerf darts, the mechanics are nearly identical (main differences in internal ballistics come from air flow, and then the external and impact ballistics are very different).
I'm interested in the work you have done here - as btrettel mentioned to me, it's the first publicly released spring powered simulation. Particularly, I wanted to look through your algorithm and see the order of your approximations.
I suppose that the first thing that I'd like to mention is that I found no errors in your kinetics or kinematics. Functionally, that part of the algorithm is very solid. However, you did make some mistakes in the pressure calculations, and I'll get into that later.
I do want to note that instead of adding an entire new column to get the initial conditions of the step (which are the final conditions of the previous step), you can just use the final conditions of the previous step (old column, but up one row). That would really clean up a lot of the equations, and make it clearer for debugging. (As an example, Column O is the same as Column AH, but shifted by one row.)
I worked out your algorithm from the excel equations. I'm sure you know it, but for everyone else, I'll post it here.
Code: Select all
Cammyd’s algorithm
Variables: dd is the length of the distance step, nf is the final step number, or total number of steps -1. a, v, d are the acceleration, velocity and position of the plunger as functions of the step number. a2, v2 and d2 are the acceleration, velocity and position of the pellet as functions of the step number.
1. Step with distance dd=(df-di)/nf, 0<=n<=nf
2. Calculate a(n) from force balance
3. Calculate v(n+1)=sqrt(v(n)^2+2*a(n)*dd)
4. Calculate d(n)=n*dd
5. Calculate dt=(v(n+1)-v(n))/a(n)
6. Calculate a2(n) from force balance
7. Calculate v2(n+1)=v2(n)+a2(n)*dt=v(n)+a2(n)*(v(n+1)-v(n))/a(n)
8. Calculate d2(n+1)=d2(n)+ (v2(n+1)+v2(n))/2*dt= v2(n)*dt+a2*dt^2/2
Pressure:
“Actual Pressure Calculator” = APC = Patm*(V_c+V_D)
A. V_c(n)=(V_c(0)-n*dd*A_c)
B. V_b(n)=d2(n)*A_b
C. P(n)=APC/(V_c(n)+V_b(n)+V_d)
In your initial conditions, you don't have units on "Actual Pressure Calculator". They should be Joules, if you weren't sure.
The issue with your pressure calculations comes in here though. Your calculation is APC = P*V = constant. This is shown when you use APC/(total volume) = pressure of each step. This, however, is an isothermal relationship, which means the air temperature never changes. In reality, the temperature will change (you're not removing any heat to keep it the same, heat can't move away fast enough in such a short time frame), so you should use the adiabatic relationship, P*(V^Gamma)=constant, where Gamma=1.4 for air.
To determine the order of your approximations, I looked at the algorithm method. This was more difficult than usual because you start by using distance steps in the plunger, then convert that to time steps in the pellet.
In step 3, you use the Euler method (aka degenerate multistep) for interval dd. This gives your plunger velocity error of O(dd) (which is pretty bad - first order). Interestingly though, it gives your plunger position no error, and the acceleration's error is a function of the pressure error, which we'll see later.
In step 5 when you calculate dt, you use Euler again, and you use the plunger acceleration and velocity to do it. This is slightly worse than O(dd). This one is an easy fix though. You used dv/da=dt, where the change in velocity and the initial acceleration are approximations, but you could use dt=(change in position)/(average velocity). This reduces error because your position is defined - not an approximation. By using the average velocity instead of the initial velocity, you also reduce the error (this is the difference between Euler method and Multistep method). Step 7 only uses that dt and the pressure, so it has about the same amount of error as dt.
Finally, step 8
is where it gets bad. While you tried to use the average velocity (multistep method), because of the way you calculated the velocity of each step, it's the same equation as d=v*t+.5*a*t^2. And since your dt is O(<dd), dt^2 is O(<dd).* This is Euler of an Euler. (first order approximation of a first order approximation of a second order ODE). If it were multistep of a multistep, your error would be O(dd^2).
What this means: Basically, it means that your muzzle velocity results should be relatively OK, which it sounds like it is. However, it means that the position at which the pellet leaves the barrel is probably off by a good amount.
The accuracy of you program is a function of your step size. So you can either use a better function (like I was talking about with the methods), and/or you can use a smaller step size.
Because you're at O(dd), by using 2000 steps instead of 20, you are about one hundred times more accurate. That said, by using O(dd^2), you'd only need 45 steps to get the same accuracy. In general though, simulations should use lots and lots of steps - way more than 20.
I do hope that you don't take this feedback too harshly. I only mean to try to help the development of your model. I don't think that any of the criticisms of the model fall back on you - as you said, it's a work in progress.
I really like the work that you've done, and am glad that you made it public - I'm certainly not comfortable making most of my work public just yet (if ever). I hope that in the future, your SGDT will be a valuable tool to help designers and to compare against the work of others.
Out of curiosity, why did you choose to use the steps to be a function of plunger position, instead of time? I've never seen that before, but it has some interesting benefits. That's one thing I'm going to look at including in my simulations.
*Edit: Btrettel pointed out to me that O(dd^2) is better than O(dd), and he's right about that, so this point was moot. Euler of an Euler is still less accurate than just an Euler though.
Edit 2: Apparently strike-through isn't supported here. I just made the font small instead.
