Computer modeling

[btrettel]

I had started writing my own computer model to help me optimize my air guns. I can use GGDT, yes, but I like to do things myself and as an engineer I feel a need to know exactly what's going on. My goal is to maximize the efficiency of a small air gun so I will use less air and lighter components.

CALM is the only available model with viewable source code so I've been comparing my model against that and looking at the source code when I've had problems. Any problems I've had so far have been minor compared to what I have now.

My model does the same process as CALM and probably GGDT. I wrote some differential equations that control flow, pressure, and the forces on the projectile. My model figures out how much air flows through the valve in an interval, moves the projectile, changes the pressures accordingly, and repeats until the projectile leaves. This is rather rudimentary, yes, but it should be plenty adequate for a small air gun like what I intend to build.

There are two problems with my model, and they likely are related:
1) My model deviates from CALM at low volumes.
2) My model can produce impossible efficiencies.

At volumes at or higher than about 5 cubic inches, my model and CALM are within 1% of each other. At lower than 5 cubic inches, they deviate a small to substantial amount. I could post some examples, but you can run my model and CALM with the same data and see what's going on too.

I've also noticed that my model can produce situations with efficiency greater than 100%. I've looked at my code over and over again and can't seem to figure out where I've gone wrong. I think this problem is related to the deviation from CALM at low volumes because the impossible efficiencies occur at low volumes with high pressures and high Cv values.

Does anyone have any ideas about where my model goes wrong? If it would be helpful I can write out which specific equations I used on paper with explanations. Below is the Python code of my model at the moment.

Code: Select all

# internal.py - Internal ballistics
from math import *

# inputs

# Vc = input("Volume of air chamber (in^3) = ")
# L = input("Barrel length (in) = ")
# d = input("Barrel diameter (in) = ")
# P = input("Pressure (psig) = ")
# W = input("Projectile weight (lb) = ")
# Pf = input("Equivalent pressure of friction (psi) = ")
# Vd = input("Dead volume (in^3) = ")
# Cv = input("Cv = ")
# Vp = input("Pilot volume (in^3) =")

Vc = 20
L = 12
d = 0.5
P = 150
W = 0.1/16
Pf = 2
Vd = 0.5
Cv = 3
Vp = 0.3

# program specs
h = 0.00001 # interval in seconds

# constants
T = 530 # atmospheric temperature in degrees Rankine
Patm = 14.7 # atmospheric pressure in psia
g = 32.2 # acceleration due to gravity in ft/s^2

# convert P in gauge pressure to absolute pressure
Pc = Patm + P;

# convert weight W in lb to mass m in slug
m = W/g;

def Vb():
	# Finds current barrel volume from X and d
	return (pi/4) * pow(d,2) * x + Vd

# before looping...
t = 0
x = 0
x_dot = 0
Pb = Patm
PVb = Pb * Vb()
PVc = Pc * Vc

while x < L and t < (0.1 - h):
	# keep accelerating the projectile until it leaves the barrel
	if Pc < Pb:
		Pc = Pb
	PVb_dot = 462.24 * Cv * sqrt((pow(Pc,2) - pow(Pb,2))/T) * Patm
	PVb = PVb + PVb_dot * h
	PVc = PVc - PVb_dot * h
	Pb = PVb/Vb()
	Pc = PVc/Vc
	x1_dot = ((12 * pi * pow(d,2))/(m * 4)) * (Pb - Patm - Pf)
	x_dot = x_dot + h * x1_dot
	x = x + h * x_dot
	t = t + h

print "Vf =", x_dot/12
print "Tf =", t
print "Pf =", Pb
ke = 0.5 * pow(x_dot/12,2) * m
pe = (1.0/12) * (P + Patm) * (Vc + Vp) * log((P + Patm)/Patm)
print "KE =", ke
print "PE =", pe
print "eta =", ke/pe
input("Done.")

I have a few theories as to where the extra energy is coming from, but CALM should make the same errors and as far as I can tell CALM doesn't. I try the same situations in CALM (sometimes after modifying CALM to allow me to enter these numbers) and I don't get any output. This is why I think the impossible efficiency comes from the difference between my model and reality, which seems to exist for my model and CALM too.

While I'm on the subject of efficiency, what are the best ways to improve efficiency? I've varied pressure and volume while keeping muzzle velocity constant (yes, it took a good deal of calculation), and found this interesting curve:



This is a line plotted in 3D space. The height is efficiency, the lower right is pressure in psi, and the lower left is volume in cubic inches. You can see the beginning of the impossible increase in efficiency at low volumes on the left--this occurs in all cases actually, but it's most pronounced at higher Cv values.

Achieving efficiencies greater than about 25% seems to be difficult. I've tried other things like increasing the Cv value, and that doesn't seem to increase efficiency beyond maybe 2 or 3 percent after a certain saturation point. Decreasing the dead volume helps too, but not substantially.

[clide]

Have you been able to determine if the problems come from low chamber volume, high pressure, a combination of the two, or either of the two? It is not immediately clear from reading through you post, but I may have missed something.

I assume this part of the code:
if Pc < Pb:
Pc = Pb
came from the this situation occurring and causing problems with the PVb_dot equation, but it looks like even when this situation is occurring you are doing the rest of the calculations with the higher than possible Pb.

You might run a quick check by adding a counter into that part of the code. I would suspect that the correction would be running more often as you move into the region where your results start to differ from other programs.

So basically you probably need to find a way to check if your barrel pressure is increasing more than possible and if so make a correction before running your x1_dot calculation.

That's my best guess anyway. Hope it helps.

[btrettel]

Thanks for the reply clide.

My post wasn't a great example of clarity in writing, so I've edited it a bit.

From what I can tell, the model has a tendency to overestimate the muzzle velocity of small chambers. This is not due to highest pressures because the model seems closer to CALM at higher pressures and lower volumes from my quick number crunching.

As for the part of code you quoted, that's not the first way I plan to do it. The problem is that the code can produce higher pressures in the barrel than in the air chamber if something doesn't keep them equal once they meet. My original plan was to make it so that if the barrel pressure is higher, the flow starts going back into the air chamber. After examining CALM, I saw that akb simply set the two pressures equal to each other once the barrel pressure exceeded the chamber pressure.

I'm going to try a counter to see how many times it does run that routine. Good idea. I'm also going to try setting the barrel pressure to the chamber pressure rather than the chamber pressure to the barrel pressure because if the barrel pressure is higher then setting the chamber pressure to that value adds some energy to the system. This sounds like the problem.... maybe if this doesn't work I'll simply add my reversal code as originally planned.

Another potential source of free energy is the isothermal decompression I assumed. To keep the temperature constant, there must be heat transfer to the system--isentropic decompression should avoid this. I use isothermal decompression for simplicity, though, it wouldn't be too difficult to try isentropic decompression and it'd be more accurate. I could then model the temperature changes and it should make my model a bit more accurate in another respect.

...hmm... I just ran some numbers that gave the impossible figures earlier and I might have found the efficiency problem. The step size (and consequently method of approximating differential equations) makes a huge difference. In my code here I'm using simple forward Euler, a simple but very error-prone method. In my external ballistic model I used the simple Runge-Kutta method, which is very significantly less prone to error. I'm going to switch to Runge-Kutta and see what else happens...

This shouldn't make my model closer to CALM however because I only changed the step size when doing some efficiency work. To figure out the efficiency at different pressures and volumes, I fixed the muzzle velocity. The model ran at set volumes to find the minimum pressure to achieve that velocity over a wide range of volumes and therefore a wide range of pressures. As this was time consuming, I increased the step size to increase speed. Not a good idea obviously.

[D_Hall]

Something that I see right off the bat.....

Code: Select all

   if Pc < Pb: 
      Pc = Pb 

You're trying to account for discrete time steps allowing too much air in the barrel. To compensate, you bump the chamber pressure to match. But that adds energy to one part the system without removing it anywhere else. Voila, >100% efficiency!

You can't just add energy in one place without subtracting it somewhere else. May I suggest something along the lines of....

Note: I don't know Python so please don't hold me to syntax. 'Tis pseudocode. Also note that GGDT tracks everything via mass flow (pressure is along for the ride rather than a driver) so, no, this isn't a tested idea and there may be a bug... But hopefully it gets the idea across.

Code: Select all

if (Pc < Pb) {
  Pc = (Pb*Vb + Pc*Vc) / (2*(Vb+Vc))
  Pb = Pc
  }

If I did that right, it will lower the pressure in the barrel while raising it in the chamber... But only to the appropriate "steady state" pressure (but I didn't check my math so DO check it before you use it!).



And for what it's worth GGDT uses a 3rd order Euler. Not as accurate as a Runge-Kutta but much more intuitive and - in my experience back in my days of writing flight sims for the military - the errors induced by it are an order of magnitude smaller than errors induced by other unknowns (example: How accurate is that pressure gage?).

Oh, and the methodology used by GGDT is very similar to CALM's methodology. My refusal to release source code has less to do with me being a control freak and more to do with past experience releasing source code... In 2 weeks there are 47 versions and everybody comes pounding on your door wanting to know why this or that isn't working. HTF am I supposed to know why YOUR modifications broke the damned code?[/code]

[btrettel]

Thanks for the reply.

I barely know Python so any psuedo-code is fine. I mainly started using Python because one of my programmer friends really loves it and I'd thought I'd try it.

Your modification should be a simple compromise and with that there won't be any energy change. I'll implement it and see how it works.

...trying out your change, I can get much closer to CALM at lower volumes. This I find odd because CALM makes the same errors I did (I got the idea to do that from CALM!). CALM must have some extra code to check something that I don't.

At this point I suppose I'll have to generate some data to fit my model to.

Any ideas on how to improve efficiency beyond 25% or so? See the end of my first post for my thoughts so far.

Edit: Because this thread title might turn some people off, I made a new thread about efficiency: http://www.spudfiles.com/forums/viewtopic,p,222630.html

[btrettel]

My efficiency thread changed my focus a bit so that I'm now focusing on higher pressures. Because higher pressures introduce non-negligible temperature drops I'm changing my model to assume an isentropic process (no heat transfer, etc.) using constant specific heats and the ideal gas law. This'll let me track the temperature changes without much difficulty.

My question at this point is whether a model as I've described is reasonable for small volumes at high pressures. I don't want to waste time on something that's incorrect. I suppose the way to know would be to try it and see, but I don't have that capability at the moment.

I'll be periodically posting code here for others to look at so please bear with me. I'm familiar with the basics behind all of this, but I'm not particularly confident so a second opinion would be appreciated. Right now I have a system of many equations (some differential) written out on paper but fitting them together coherently into a program will be a pain, so it might be a few days before I have anything to show for this.

In the future I plan on tailoring my model for the specific valve design I'm working on. Right now the valve opens instantly, which is an idealization but should be reasonable on reasonably fast valves. Put a free-body-diagram on the valve and see how it moves... this is what I think GGDT does.

Let me reiterate something again just in case it's necessary. Using GGDT would be simpler, but I want to learn exactly what's going on and GGDT doesn't do exactly what I want to do (optimization, designing for efficiency). This doesn't have to do with the lack of source code--I'd look at it but not really have a good idea of what's going on unless it's heavily commented.

[SEAKING9006]

Hmmm, maybe a bit off topic, but I was hoping a few spudfiles users knew python. I was hoping we might be able to have an old school ro-sham-bo program contest some time.

[btrettel]

Programming contests sound like fun, but my skills are barely adequate for the computer modeling I'm doing, so I'd have to decline.

I rewrote my model to assume an isentropic process and the results are pretty close to my old ones. For some simple figures I put in this model said 141 fps, my old one said 144 fps, and CALM said 152 fps. They're all within 10% of each other, which is good.

The problem is the temperatures. The barrel should heat up and the air chamber should cool down. But the differences here seem enormous--the air chamber will cool down 100 degrees Rankine/Fahrenheit and the barrel will get 35 degrees hotter at 30 psi. Is this reasonable? I put 300 psi in the program just to see what it would say and the barrel temperature increases to about 1000 degrees Rankine (~ 500 degrees Fahrenheit). My guess would be that the 30 psi figures are accurate but the constant specific heats and ideal gas assumptions cause some oddities with the higher pressures. Edit: Hmm... I don't think the temperature in the barrel should increase so perhaps I need something else added to this...

The code is below. If necessary I can give derivations of the equations used.

Code: Select all

(code removed because of errors)

Testing the limits of the model, I've noticed that this can produce impossible efficiencies at low temperatures (like 400 degrees Rankine and below) too. This probably has to do with the ideal gas assumptions and use of constant specific heats, but I don't know for sure. Any suggestions or errors spotted are welcome.

Edit #2: My way of figuring out the new temperatures in this is incredibly wrong. I think I was caught up in it not being isothermal and that there's no heat transfer so I just looked at isentropic/adiabatic process equations. All I need to do is figure out the final temperature for the mix of the two different air masses at two different temperatures.

Learning through trial and error can be annoying but it seems to be working so far.

[btrettel]

After a few days of messing around with my code, I think I finally have the temperature drop in the barrel modeled at least approximately correctly. The code is at the bottom.

I basically treated to the gas in the barrel as a closed system without any heat transfer. Using Cp I could figure out the internal energy of the system and I knew the work done (approximately at least), so this wasn't too difficult. My results are moving towards GGDT's now... earlier my model had predicted about where CALM was but at least on the data I'm using now mine is fairly close to GGDT. Below are the results I'm looking at:

btrettel new: 172 fps
btrettel old: 199 fps
CALM: 203 fps
GGDT: 175 fps

My model predicts the temperature drops more than GGDT says it will, but I blame that on constant specific heats and my way of calculating the work done being approximate.

At higher pressures my model starts to deviate from GGDT significantly because I neglect Mach effects, use constant specific heats, and simply don't have as robust a model. But I'm getting there, slowly but surely. By the time I actually build what I intend to my model should be fairly robust.

I do have one question some here might be able to answer. I use Cp, the specific heat at constant pressure, to calculate the temperature drop. While I did model the system as closed, which should mean I use internal energy, the results weren't close to GGDT if I used Cv to figure out internal energy. I figured that makes sense because the system really isn't closed, so I'm using Cp now and find it agrees with GGDT fairly well. But this doesn't make perfect sense. I think I'll email my old Thermo TA... he was helpful about things like this.

Code: Select all

# bags.py - btrettel's air gun simulation
from math import *

# inputs
Cvv = 6
d = 1.590 # barrel diameter, in
L = 18 # barrel length, in
P = 50 # air chamber pressure, psig
Pf = 2 # minimum pressure to move projectile, psia
# Px = 50 # pressure to close valve, psia
Tc = 530.0 # initial air chamber temperature, R
Vc = 40 # air chamber volume, in^3
Vd = 10 # dead volume, in^3
Vp = 0.1 # pilot volume, in^3
W = 0.0749571691 # weight of projectile, lb

# constants
Cp = 2241.1 # specific heat at constant pressure, in/R
Cv = 1606.1 # specific heat at constant volume, in/R
g = 32.2 # gravitational acceleration, ft/s^2
k = 1.4 # ratio of the specific heats
lbtocg = 45360 # multiply by something in pounds to convert to centigrams
Patm = 14.7 # atmospheric pressure, psia
R = 640.12 # gas constant for air, in/R
Tr = 528 # SCFM reference temperature, R
T = 530.0 # atmospheric temperature, R
SG = 1 # specific gravity of the gas at atmospheric pressure and 60 degrees F

# convert weight to mass
m = W/g

# program specs
h = 1e-6 # interval in seconds

# functions
def Vb(xi):
	# finds the current barrel volume
	return (pi/4) * pow(d,2) * xi + Vd

def Q(Pi, Po, Tci):
	# finds the volumetric flow rate at STP
	return 462.24 * Cvv * sqrt((pow(Pi,2) - pow(Po,2))/(Tci * SG))

def x_2(Pbi):
	# finds the current projectile acceleration
	if x_dot > 0:
		return ((12 * pi * pow(d,2))/(4 * m)) * (Pbi - Patm - Pf)
	else:
		return ((12 * pi * pow(d,2))/(4 * m)) * (Pbi - Patm + Pf)

def Tt(mi,mo,Ti,To,Poi,worki):
	# finds the new temperature with mixing, assuming no heat transfer
	return (mi * Ti * h + mo * To)/(mi * h + mo) - (worki)/(Cp * (mi * h + mo))

def m_dot(Qi):
	# find mass flow rate from Q
	return (Patm * Qi)/(R * Tr)

def Ps(mi, Ti, Vi):
	# find current pressure
	return (mi * R * Ti)/Vi

def E(Pmax,Pmin,Vmin):
	# figures out energy stored in a gas
	return (1.0/12) * (Pmin * pow(Pmax/Pmin,1/k) * Vmin - Pmax * Vmin)/(1 - k)

# initial conditions before looping
t = 0
x = 0
work = 0
x_dot = 0
Pb = Patm
Pc = P + Patm
Tb = T
# Tc = T
mb = (Pb * Vb(x))/(R * Tb)
mc = (Pc * Vc)/(R * Tc)

while x < L and t < (0.1 - h):
	
	Pbo = Pb
	Pco = Pc
	Tbo = Tb
	Tco = Tc
	
	if Pb <= Pc: # if barrel pressure <chamber> flow into the barrel
		m_dota = m_dot(Q(Pc, Pb, Tc))
		Tb = Tt(m_dota,mb,Tc,Tb,Pb,work)
		#Tc = Ts(mc,Vc,Pco,Tco)
		mb = mb + m_dota * h
		mc = mc - m_dota *h
	else: # if barrel pressure > chamber pressure --> flow into the chamber
		m_dota = m_dot(Q(Pb, Pc, Tb))
		Tc = Tt(m_dota,mc,Tb,Tc,0)
		#Tb = Ts(mb,Vb(x),Pbo,Tbo)
		mb = mb - m_dota * h
		mc = mc + m_dota *h
	
	Pb = Ps(mb, Tb, Vb(x))
	Pc = Ps(mc, Tc, Vc)
	
	x0 = x
	
	x_dot = x_dot + h * x_2(Pb)
	x = x + h * x_dot
	
	# work = Pb * (Vb(x) - Vb(x0))
	
	b = Pb - (Vb(x) * (Pb - Pbo))/(Vb(x) - Vb(x0))
	a = (Pb - Pbo)/(Vb(x) - Vb(x0))
	
	work = (a/2) * (pow(Vb(x),2) - pow(Vb(x0),2)) + b * (Vb(x) - Vb(x0))
	
	t = t + h

print "Vf =", x_dot/12
print "tf =", t
mair = lbtocg * (P * (Vc + Vp))/(R * T)
print "Moved air mass, with pilot volume (cg) =", mair
mair3 = lbtocg * (P * Vc)/(R * T)
print "Moved air mass, zero pilot volume (cg) =", mair3
mair2 = lbtocg * ((P - Pc + Patm) * Vc)/(R * T)
print "Moved air mass, new valve (cg) =", mair2
print "Old valve vs. no pilot  =", (1 - mair3/mair)
print "New valve vs. old  =", (1 - mair2/mair)
print "New valve vs. no pilot old  =", (1 - mair2/mair3)
ke = 0.5 * pow(x_dot/12,2) * m
print "KE =", ke
pe = E(P + Patm,Patm,Vc + Vp)
print "PE =", pe
pes = E(Pc,Patm,Vc)
print "PE saved =", pes
print "eta =", ke/pe
pe2 = E(P + Patm,Patm,Vc)
print "eta (no pilot) =", ke/pe2
print "eta2 =", ke/(pe2 - pes)
print "	Barrel		Chamber"
print "T:" + "\t" + str(Tb-459.67) + "\t" + str(Tc-459.67)
print "P:" + "\t" + str(Pb-Patm) + "\t" + str(Pc-Patm)
input("Paused...")

[btrettel]

I've added a bit more to the simulation and am waiting for some suggestions or corrections from my old thermodynamics TA. While the model is still crude in some respects, it is somewhere between CALM and GGDT in terms of accuracy now, which should be adequate for what I intend to do at the moment. Checking against GGDT will ensure it is adequate.

To avoid reposting code here a trillion more times I'll have the latest version available here in the future: http://trettel.org/bags.html

My goal with this is to make an open source and cross platform computer model that's accurate with respect to its limitations. There are too many plans for the model at the moment, but it's safe to assume it'll be released in the future with a GUI, a more accurate theory behind it, and more features. The most major deviation from GGDT will be that I'm designing the model to be easy to try out many different configurations with loops, primarily to optimize pneumatic guns for efficiency.

Let me reiterate something again to: this is not intended to replace GGDT, rather, just to be an alternative primarily for people like me who want to tinker with different configurations or have a better idea of what exactly is going on. GGDT is significantly more robust than my model and that robustness shows in its accuracy.

Any suggestions, comments, etc., I would suggest be posted at trettel.org or emailed to me.

[rp181]

I keep hearing about these runge-katta methods. What are these? Ive seen them used everywhere from this to barry's resistance-inductance-capacitance simulator. I can't find any good information on this. Someone want to point me to one?

[btrettel]

Wikipedia explains the basic 4th order Runge-Kutta method in all the detail you should need.

If you're unfamiliar with differential equations you probably won't have any idea what they're talking about, but if you've taken at least some calculus you should follow the discussion.