Modeling Air Compression Cylinder

[JeroenM]

TL;DR Summary

Trying to model a piston that compresses air and pushes it out of a nozzle

Hi,

I'm have been trying to model an airsoft replica, as an electrical engineer fluid dynamics, thermodynamics is a complete new flied for me and with with help of Ai i have been learning a lot about this subject. So far I have come a long way in making something my first model sort of worked but the output energy was about 40% higher then actual measurements. So i have been trying several things to see where the losses are and how to model them using Matlab. To make things a bit simpler I removed as many components as possible to make the model in steps.

So the first part I'm now trying to model is a cylinder (with diameter of 23.85mm) with a outlet on one end that is 4.5mm in diameter and 35mm long and on the other end there is a piston which will compress the air and push it out of the smaller hole. The piston is driven by a spring (150mm spring with 692N/m springconstant and pre compressed to a length of 55mm).

So the model looks like this

I measured the pressure inside the chamber using an absolute pressure sensor and the graph looks like this:

The pressure peaks at 350kPa when the piston slams into the end of the cylinder from what I conclude is that the flow gets choked at 250kPa (more rapid rise).

From what I so far know and assume is that all the "losses" is due to the friction force of the oring (to create the seal of the piston) which is pressure depended and the fact that the air needs to be pressed through a hole. For now I modelled this using the orifice formula which i'm not sure if the correct way I also tried the sudden contraction formula but that doesn't seem to fit at all.

So I made the following code in matlab (snippet):

1. Calculate the volume in the cylinder based on piston position (t=0 pos = 0, Lc is the cylinder length of 55mm)
2. Calculate the pressure in the chamber using the air mass temperature and volume
3. Calculate the temperature rise due to rapid compression
4. Calculate the air density based on the new temperature and pressure
5. Assume the air velocity inside the cylinder the same as the piston velocity
6. Check if the flow is subsonic or choked flow based on the pressure of the cylinder and atmospheric pressure (P2 is assumed to be P ATM)
7. Calculate the velocity and mass air flow through the "orifice" with either subsonic flow or choked flow. this part I'm not sure if correct because it is not really an orifice plate but more a sudden contraction to a smaller pipe (of 35mm length) in to to atmosphere but so far it seems to fit the best.
8. Calculate the new air mass inside the cylinder
9. Calculate the force of the spring using hooks law (1/3 of the mass is added to the total mass of the piston)
10. Calculate the friction force of the oring based on the pressure inside the cylinder and P0 (P ATM). I use for mu a value of 2 (which seems to be way to high a value of 0.5 or lower is what I have expected) and the A_seal is the area of the oring which is for now the circumference (33.4mm) of the oring times the width (2.58mm) of the oring (also not sure if this is correct).
11. Calculate the pressure force on the piston based on the area of the piston head and the delta pressure on the piston (P_cylinder - P0)
12. Calculate netto force on the piston.
13. Calculate the acceleration, velocity and position of the piston
14. Check if piston is at the end (position > 55mm)

for n = 2:n_steps
    % Chamber volume
    V_c(n) = Lc*Ap - x_p(n-1) * Ap;

    % Update chamber pressure (Ideal Gas Law)
    P(n) = m_air_chamber(n-1) * R * Temp(n-1) / V_c(n);

    Temp(n)=T*(P(n)/P0)^((gamma-1)/gamma); % Calculate temperature rise due to fast adiabatic compression
 
    % Update air density in chamber
    rho_air_chamber(n) = P(n) / (R * Temp(n));
    v_air_chamber(n) = v_p(n-1); %Air velocity in chamber same as piston

    % Density at the orifice (isentropic relationship)
    P_orifice(n) = P(n) - (0.5 * rho_air_chamber(n) * v_air_chamber(n)^2); % Simplified pressure drop to orifice
    rho_orifice = P_orifice(n) / (R * Temp(n));
  
  
    P2=P0; %P2 is Atmospheric pressure

    % Calculate orifice flow
    delta_P_orifice(n) = P(n) - P2;  % Pressure difference across orifice
    pressure_ratio = P2 / P(n);  % Downstream to upstream pressure ratio


    if pressure_ratio > critical_pressure_ratio  % Subsonic flow
        flow_regime(n) = 1;  % Subsonic
        if delta_P_orifice(n) > 0
            v_air_orifice(n) = sqrt( (2*gamma*R*Temp(n))/(gamma-1) * (1-(P2/P(n))^((gamma-1)/gamma)));
            m_dot_orifice(n) = C_d * An * P(n) * sqrt(2 / (R * Temp(n))* (gamma / (gamma - 1))  * ((P2/P(n))^(2/gamma) - (P2/P(n))^((gamma+1)/gamma)  ) );
            %*C_d_f(n)
        else
            m_dot_orifice(n) = 0;
        end
    else  % Choked flow
        flow_regime(n) = 2;  % Chokedflow_regime(n) = 2;  % Choked
        v_air_orifice(n) = sqrt(gamma*R*Temp(n)*(2/(gamma+1)));
        m_dot_orifice(n) = C_d * An *P(n)*sqrt(gamma/(R*Temp(n))) * (2/(gamma+1))^((gamma+1)/(2*(gamma-1))); 
    end

    % Update air mass in chamber and barrel
    m_air_chamber(n) = m_air_chamber(n-1) - m_dot_orifice(n) * dt;  % Air leaving chamber

    if m_air_chamber(n)<0
        m_air_chamber(n)=0;
    end


    % Spring force (Hooke's Law with pre-load)
    F_spring(n) = k * (x_preload - x_p(n-1)); % Spring force with pre-load
 
    % Friction force (based on seal)
    F_friction(n) = (P(n) - P0) * A_seal * mu;  % Friction force on the piston
  
    %Force due to pressure buildup
    F_pressure(n) = (P(n)- P0)* Ap;

    %Net force on piston
    F_n(n) = F_spring(n) - F_friction(n) - F_pressure(n);

    % Piston dynamics
    a_p(n) = F_n(n) / m_p;
    v_p(n) = v_p(n-1) + a_p(n) * dt;
    x_p(n) = x_p(n-1) + v_p(n) * dt;
  
    % Apply saturation to piston position
    if x_p(n) < 0
        x_p(n) = 0;
        v_p(n) = 0;
    elseif x_p(n) > Lc
        x_p(n) = Lc;
        v_p(n) = 0;
        break;
    end
end

So below is a picture of the simulation vs the measured data of the pressure inside the cylinder

So the simulation is pretty close except that the peak is missing in the simulation, so either the simulation is wrong or the measurement is wrong (pressure is measure through a hole on the side of the cylinder right before the endstop).

I used CFD simulation to see what I could expect and got the following results:

The inlet is a mass air flow of 50gr/s and the outlet is a static presure of 101325 pascal. From the total pressure inside the cylinder and ATM I have a pressure difference of 74kPa which I used to calculate the Coefficient of discharge value which is about 0.85. But the simulation gives different pressure values (static, dynamic, relative and total pressure) for the inlet (on the right) and the outlet (on the left) I'm not sure any more which one to use especially of the left side (pressure average is 101325, total pressure average is 151371 pascal so not sure which one to use). If use the total pressure only the pressure difference between inlet and outlet is only 23888 pascal which seems very low and would not explain why the pressure inside the cylinder is 175kPa. Here is really the lack of knowledge to understand what is happing here.

So to sum it up it up where I have my doubts:

• Doubt if the orifice formula is the correct one, but if so not sure if i should P0 as value inside the formula especially because the choked flow seems to start at 250kPa in the measurement.
• There is friction force inside the system but based on my model is seems to be very high, friction coefficient of 2 seems to be very high and unrealistic. Tried to simulate this using FEM but not much succes there or very low friction numbers.

Who can help me to point where I might be wrong ?

Best regards,

Jeroen

 

[erobz]

Maybe just try to focus in on what is happening to the piston. You have a known force from the spring (at least idealistically), known piston mass. What you don't know is the contact friction force from the seal of the piston from cylinder and the internal pressure acting on the piston face. Unfortunately, I think that pressure is very complex behaviorally. My guess is the gas is basically rapidly compressing before any flow initiates. I think the flow viscosity would be critical to that.

What do you get for the pressure in the nozzle if you compress it adiabatically the initial volume cylinder + nozzle volume down to the volume of the nozzle?

 

[JeroenM]

If the i do P1V1^1.4=P2V2^1.4 would be P1*V1^1.4/V2^1.4= 20MPa.

The part of the piston with regarding spring force and mass are easy to model but without friction of pressure it would be to fast.

 

[erobz]

If the i do P1V1^1.4=P2V2^1.4 would be P1*V1^1.4/V2^1.4= 20MPa.

The part of the piston with regarding spring force and mass are easy to model but without friction of pressure it would be to fast.

Ok, so if that is correctly calculated at 20 MPa, and you are measuring peak 0.4 MPa then there must be some flow initiating prior to complete adiabatic compression of the initial volume. The question is how much because of the high level of nonlinearity...

P.s. it would be good if moving forward ( with whoever you work with) if you could type equations in Latex

Here is the code:

$$ P_1V_1^{1.4}=P_2V_2^{1.4} \implies P_1\left(\frac{V_1}{V_2} \right)^{1.4}$$

That renders as:

$$ P_1V_1^{1.4}=P_2V_2^{1.4} \implies P_2 = P_1\left(\frac{V_1}{V_2} \right)^{1.4}$$

See: LaTeX Guide

 

[erobz]

The Darcy Weisbach Equation is developed under the assumption of steady - fully developed- incompressible flow. I have worries there for applicability to this type of problem as it is literally everything but those assumptions.

The thing is if it is reaching 0.4 MPa viscosity has to be responsible IMHO. If there were absolutely no viscosity, the gas should just slip out the barrel without a pressure gradient. However, I'm not sure to what extent purely accelerating the mass of air in front of the piston out the barrel is responsible for the pressure rise. Trying to apply Bernoulli's ( steady - inviscid flow) here to find a pressure is also a mistake in my estimation for some of the reasons above - very unmatching problem assumptions validate Bernoulli's.

I have my doubts on (certainly) my ability to produce a formula. I hope some fluids experts will chime. just thought hopefully to stimulate some thought while waiting.

 

[JeroenM]

I just found an error in the mass flow rate calculation , I updated the main post. There are small differences but that the big peak is not in de simulation. Could be a measurement error.

So lets assume the calculations are correct then I still am doubting about the high friction coefficient of the oring. Im planning to make a cylinder with movable piston on 2 sides which can be pressurized to measure the friction force while under pressure.

Next step would be to add the barrel to the nozzle and model the losses there which I was planning to do based on the darcy-weisnbach and pipe-loss but as you mention that doesnt work for compressible flows. Maybe have to look more at the fanno flow model.

 

[jrmichler]

There are small differences but that the big peak is not in de simulation. Could be a measurement error.

It's probably not measurement error. The piston is getting close to the end wall, and restricting the air flow through the orifice. The air flow out is restricted, the piston has inertia, so the pressure rises rapidly. It's nice to see that your pressure measurement / data acquisition has enough bandwidth to see the pressure spike.

The mechanical vibration from the piston hitting the end could cause an error in the pressure sensor. You could try isolating the sensor with a very short length of rubber hose, but the increased dead volume would cause other problems.

I did not check your calculations, but the approach you listed in Post #1 looks correct. I once wrote Matlab simulation code for a similar problem. My project was to move a 45 pound mass 9 inches in less than 150 milliseconds without slamming or bouncing. We were able to move that mass smoothly and quietly in as little as 120 msec.

My code output very nicely matched measured pressures and piston accelerations. I modelled the output as an orifice with either subsonic or choked flow. The flow restriction at end of stroke was not a problem because my model was for using an air cylinder to decelerate a moving mass, so the velocity approached zero at end of stroke. I assumed polytropic air compression. Your process happens faster, so your use of adiabatic looks good to me. My model did correctly simulate the piston bouncing off the air cushion when the output restriction was too small, and slamming into the end when the output restriction was too large.

 

[JeroenM]

yeah could be that it isnt an measurement error in the rest of the measurement data you can see small oscillations in pressure due to the piston bouncing, I assume the Cd value drops probably very hard at the end. What pressure value did you use for the mass flow calculation also atmospheric pressure ?

In the end application i don't see that happening at all because the pressure will be moving the bb which slows everything down. Do you known something regarding the Cd value and the flow simulations ?

This is the measurement with a barrel of 300mm connected (diameter is 6.05mm, slightly larger then the 4.5mm of the nozzle)

 

[erobz]

What is the mass of the piston, I'm trying to do some sanity check calculations?

 

[jrmichler]

If you have access to an accelerometer, you can do a sanity check for the end of stroke pressure spike. Measure the cylinder acceleration in the direction of piston motion for a normal stroke. Locate the accelerometer as close as possible to the pressure sensor. Then tap the cylinder with a soft hammer until you can approximately duplicate the magnitude and duration of the end of stroke acceleration, while measuring the output of the pressure sensor. Remove the piston for this test. Any resulting output from the pressure sensor is due to the acceleration.

My restriction was a needle valve. I made no attempt to calculate a $C_v$ for the needle valve. We adjusted the needle valve, measured the acceleration of the system, and adjusted the polytropic coefficient and $C_v$ value in the program until the simulation results matched the measurements. The simulation results included both cylinder pressure and system acceleration. That's my recollection, this work was done 20 years ago, and all notes were left behind when I left the company.

 

[JeroenM]

mass of the piston is 28 gram and the spring is 9 gram.

Already tried the accelerometer but data wasn't very useful or repetitive

So what I think happens the Cd value drops when the piston comes near the end and therefor the pressure rises

Below is a plot of the complete system

 

[erobz]

mass of the piston is 28 gram and the spring is 9 gram.

Already tried the accelerometer but data wasn't very useful or repetitive

View attachment 355105

So what I think happens the Cd value drops when the piston comes near the end and therefor the pressure rises

View attachment 355115

Unfortunately, I don't believe the discharge coefficient $C_d$ derives from first principles with any simplicity (maybe not at all?). So, you might be stuck with empirical results or CFD models.

You seem to know how to use tech to get results, I'm not really feeling the putting much effort there. However if we just ignore compressibility, viscosity, friction, etc.( ) I get the ODE:

$$ \frac{- \pi D_p^2}{8} \rho \left( \frac{D_p^4}{D_b^4} - 1 \right) \dot x^2 + k \left( x_o - x \right) - \left( M_p + \frac{1}{3}m_s \right) \ddot x = 0 $$

Where
$D_p$ Diameter of piston
$D_b$ Diameter of barrel
$M_p$ Mass of piston
$m_s$ Mass of spring
$\rho$ density of air at some point between atmospheric and 0.4MPa

I don't think its analytical, but maybe you can get AI to pump out a result numerically in a few seconds! Maybe you have done this one already - I'm just trying to hack away at it.