Trouble solving for end state of two control volumes in a rigid tank

[calzone41]

TL;DR Summary: Struggling to structure the problem and derive an analytical solution for gas expanding into other gas in a rigid tank. Preferred formulation is fixed control volumes. This is not a homework problem.

The problem:
Two control volumes (A and B) are in a rigid tank filled with air. The initial state of each control volume is completely defined. There is a pressure differential between the control volumes so mass flows from A to B. Mass flows such that the pressure in each control volume is equal at time 2. There is no energy transfer or work to or from the tank; frictional losses are ignored. Mass and energy are conserved. There is no heat transfer between air masses except for mass transfer across the control volume boundary. Assume immediate and perfect mixing as mass moves into control volume B. I need to be able to calculate the final state of the air masses in the control volumes (P, m, T). This is straightforward to do numerically, I would like an analytical solution.

Some Background:
This is not a homework assignment. Ultimately I am trying to build a discrete-time-step Stirling engine simulator that models each air mass separately (imagine a gamma configuration engine with hot and cold air masses in a displacer cylinder and a power air mass in a power cylinder.) At each time step, various air masses gain/lose energy and perform work, so will have different pressures. I need to balance these pressures at the end of each time step. Eventually, I will be modeling three connected air masses, but wanted to understand the two-airmass problem first before adding complexity.

Attempted solutions
Systems of Equations:
Using fixed-size control volumes, there are 5 unknowns at the final state (mA, mB, TA, TB, P2). I find 4 equations: conservation of mass, conservation of enthalpy ( that is; hA1mA1 + hB1mB1 = hA2mA2 + hB2mB2. Is this a valid equation?), and ideal gas law for both volumes at the end state. This would not be enough to solve for 5 unknowns. Is there another obvious one? I considered relationships for polytropic processes, but was not confident they apply for fixed control volumes. I did derive a 5th equation that relates internal energy to mass flow for control volume A (shown below) but I'm hesitant about it. Solving this non linear system of equations is beyond me. Perhaps you can point me in the right direction?

Control volume energy balance, integral with respect to time:
Starting with the general energy balance equation for volume A, I eliminated terms for heat flux, work, and mass flow in. I then differentiated (chain rule) um with respect to t. The m-dot term here makes no physical sense (mass does not just disappear) so is zero. Using the definition of enthalpy, the ideal gas law, the specific heat relationship u = CvT, the definition of k = Cp/Cv, and integrating both sides of this separable differential equation. I arrive at internal energy as a function of mass flow and time. assuming mass flow is constant this can be re-written as internal energy as a function of mass in the control volume. (Note: It is possible to perform this integration as du/dm rather than du/dt. The result is the same.) This is an simple, elegant equation (which I cannot verify externally). However, when using this relationship to assemble a differential equation for control volume B, I get an unwieldy equation I cannot figure out how to integrate.

Combined control mass and control volume:
I can also structure the problem as control masses instead of control volumes. Essentially, mass A pushes into mass B, doing work Imagine a piston between the air masses. Using polytropic relationships for isentropic flow, and the relationship of total volume being constant, I can solve for P2. I can then use this P2 and the control volume relationships to solve for the changes in mass, then find other equations of state. This is fine, but feels messy and kludged.

Questions:
Am I making this problem more complicated than it needs to be? Is there a simpler way to solve the problem of gas expanding into other gas in a rigid tank?

Side question: Is conservation of enthalpy a thing? When using this relationship I can get an equation for P2 that is essentially an arithmetic mean [P2 - (Pa1Va + Pb1Vb)/(Va + Vb)], which differs from the P2 found using control masses and polytropic relationships and seems too simple to be real.

Any help in how to think about this problem would be appreciated. Thanks in advance!

 

[anuttarasammyak]

As a first try I take it easy as follows.
Say B is vacuum, gas in A expands free adiabatically. Temperature does not change from $T_A$ . $P_A$ drops to $\frac{V_A}{V_A+V_B}P_A$.
Say A is vacuum, gas in B expands free adiabatically. Temperature does not change from $T_B$. $P_B$ drops to $\frac{V_B}{V_A+V_B}P_B$.
In total the air has pressure $\frac{V_A}{V_A+V_B}P_A+\frac{V_B}{V_A+V_B}P_B$. The air consists of two temeratures air of $T_A$ and $T_B$. I am not sure such coexistance of high temperature gas and low one is realistic.

 

[calzone41]

As a first try I take it easy as follows.
Say B is vacuum, gas in A expands free adiabatically. Temperature does not change from $T_A$ . $P_A$ drops to $\frac{V_A}{V_A+V_B}P_A$.
Say A is vacuum, gas in B expands free adiabatically. Temperature does not change from $T_B$. $P_B$ drops to $\frac{V_B}{V_A+V_B}P_B$.
In total the air has pressure $\frac{V_A}{V_A+V_B}P_A+\frac{V_B}{V_A+V_B}P_B$. The air consists of two temeratures air of $T_A$ and $T_B$. I am not sure such coexistance of high temperature gas and low one is realistic.

Thank you for your quick response. I wish for such a simple solution, but in my case, the gas in one control volume is expanding into other gas, not a vacuum; work is done on or by another gas and temperature will not be constant between initial and end states. This is truly the root of the complexity of this problem.

 

[Chestermiller]

Are you assuming that there is a valve that admits gas from A to B relatively slowly? If that is the case, then $$\frac{d(u_Am_A)}{dt}=+\dot{m}_Ah_A$$and$$\frac{d(u_Bm_B)}{dt}=-\dot{m}_Ah_A$$ with $\dot{m}_A$ negative.

 

[Chestermiller]

View attachment 323890

The above is incorrect. Using my notation with $\dot{m}_A=\frac{dm_A}{dt}$, we get: $$u_Adm_A+m_Adu_A=h_Adm_A$$then$$m_Adu_A=m_AC_vdT_A=(h_A-u_A)dm=P_Av_Adm=RT_Adm_A$$

 

[calzone41]

Hey folks, sorry for the delay. Since I was last working on this, my wife and I had a baby. As you can imagine, things have been busy. I have finally had a chance to return to this problem and have made some progress but have also run into a roadblock that indicates I've made a mistake somewhere in my analysis. Chestermiller, thanks for your post a year ago; you pointed out an error I made and that helped me get to where I am now.

In summary, I believe I have partially solved this problem analytically, but when I run the analytical solution next to a numerical solution (discrete time step simulation), the results diverge substantially. If anyone can help check my math or point out an error in conceptualizing this problem I would be very grateful.

First, the analytical solution.
Assume that Pa > Pb so mass only flows from Volume A to Volume B.

For Volume A:
$$ \frac d {dt} \left( um \right) = \dot m h; \text{ where } \dot m \text{ is negative} $$
using the chain rule and re-arranging we get:
$$ \dot u m =\dot m \left(h-u\right)$$
with $h=u + \frac {PV} m$, the ideal gas law $PV = mRT$, and using Leibniz'a notation, we can use the "inverse chain rule" to get:
$$\frac {du} {dm} = \frac {RT} m $$
From there, using $du = C_VdT$ and separating the differentials we get:
$$\frac {dT} {T} = \frac {R} {C_V} \frac {dm} m $$
Integrating the definite integral and using the relationship $\frac R {C_V} = k-1$ and the simplifying assumption that $h=C_P*T$ we can get:
$$\frac {T_2} {T_1} = \left( \frac {m_2} {m_1}\right) ^{k-1} = \frac {h_2} {h_1} $$
This is an elegant solution and seems reasonable. Using the ideal gas law PV=mRT we can get a similar equation for pressure:
$$\frac {P_2} {P_1} = \left( \frac {m_2} {m_1}\right) ^{k}$$
Anyone want to check my math on the above?

Solving for Volume B:
$$ \dot u_B m_B =\dot m_B \left(h_A-u_B\right); \text{ where } h_A = f\left(m_B\right) \text{ because } m_B = m_{total} - m_A$$
Attempting to rearrange this equation (after substituting in the the result for $h_A$ found above) results in an equation that is not separable. It is a first order differential equation but doesn't have a straightforward solution. Perhaps Laplace transforms can be used, but it did not seem that standard transforms apply and this is stretching my mathematical abilities until I can dig in a bit more.

I had a thought that perhaps because the system as a whole is closed, enthalpy would be conserved. This allows for an algebraic solution, but I'm not confident that enthalpy is a conserved property. A bit of internet research has not answered this question directly. Anyone have a thought on whether enthalpy is conserved in this situation (or ever)?

Solving for Volume B is the primary roadblock in my search for an analytical solution. Any help in this area would be appreciated.

Moving to a numerical solution.
As above, we have relationships for $\dot u_A=\dot m_A\left(h_A-u_A\right)$ and $\dot u_B=\dot m_B\left(h_A-u_B\right)$ (where $\dot m_A=-\dot m_B$)
By setting $\dot m$ constant, approximating differentials with $\frac {du} {dt}=\frac {\Delta u} {\Delta t} = \frac {u_2 - u_1} {\Delta t}$ (and similar for m) we get the following equations that relate the state in time step 2 to the state in time step 1 (assume that the state in time step 1 is fully defined):
Volume A:
$$u_2 = \Delta t \dot m\left( h_1-u_1\right) + u_1 $$
$$m_2 = \dot m \Delta t + m_1$$$$T_2 = \frac {u_2} {C_V} $$
$$P_2 = \frac {m_2RT_2} V $$
$$h_2 = C_P T_2 $$
Similarly for Volume B starting with:
$$u_2 = \Delta t \dot m \left( h_A1 - u_B1 \right) = u_1$$
... and omitting the other relations for B to preserve a shred of brevity for this post.

I put together a simulation using Euler's method to model the change in pressure over time (or, over $\Delta m = \dot m\Delta t$). The simulation can be found in this google spreadsheet. The chart below shows the change in pressure for each volume. ($\dot m$ is set to a constant, not as a function of pressure, so mass continues to move from A to B, even after the pressure in volume A is less than the pressure in volume B). It also shows the analytical solution to the pressure in volume A. Much to my disappointment, the analytical solution differs significantly from the numerical solution, suggesting that I have made an error somewhere above.

Is this simply an error of mathematics? Or have I conceptualized this problem incorrectly? I'll continue my work to find out, but any help from the community here would be appreciated.

 

[Chestermiller]

Hey folks, sorry for the delay. Since I was last working on this, my wife and I had a baby. As you can imagine, things have been busy. I have finally had a chance to return to this problem and have made some progress but have also run into a roadblock that indicates I've made a mistake somewhere in my analysis. Chestermiller, thanks for your post a year ago; you pointed out an error I made and that helped me get to where I am now.

In summary, I believe I have partially solved this problem analytically, but when I run the analytical solution next to a numerical solution (discrete time step simulation), the results diverge substantially. If anyone can help check my math or point out an error in conceptualizing this problem I would be very grateful.

First, the analytical solution.
Assume that Pa > Pb so mass only flows from Volume A to Volume B.

For Volume A:
$$ \frac d {dt} \left( um \right) = \dot m h; \text{ where } \dot m \text{ is negative} $$
using the chain rule and re-arranging we get:
$$ \dot u m =\dot m \left(h-u\right)$$
with $h=u + \frac {PV} m$, the ideal gas law $PV = mRT$, and using Leibniz'a notation, we can use the "inverse chain rule" to get:
$$\frac {du} {dm} = \frac {RT} m $$
From there, using $du = C_VdT$ and separating the differentials we get:
$$\frac {dT} {T} = \frac {R} {C_V} \frac {dm} m $$
Integrating the definite integral and using the relationship $\frac R {C_V} = k-1$ and the simplifying assumption that $h=C_P*T$ we can get:
$$\frac {T_2} {T_1} = \left( \frac {m_2} {m_1}\right) ^{k-1} = \frac {h_2} {h_1} $$
This is an elegant solution and seems reasonable. Using the ideal gas law PV=mRT we can get a similar equation for pressure:
$$\frac {P_2} {P_1} = \left( \frac {m_2} {m_1}\right) ^{k}$$
Anyone want to check my math on the above?

Solving for Volume B:
$$ \dot u_B m_B =\dot m_B \left(h_A-u_B\right); \text{ where } h_A = f\left(m_B\right) \text{ because } m_B = m_{total} - m_A$$
Attempting to rearrange this equation (after substituting in the the result for $h_A$ found above) results in an equation that is not separable. It is a first order differential equation but doesn't have a straightforward solution. Perhaps Laplace transforms can be used, but it did not seem that standard transforms apply and this is stretching my mathematical abilities until I can dig in a bit more.

I had a thought that perhaps because the system as a whole is closed, enthalpy would be conserved. This allows for an algebraic solution, but I'm not confident that enthalpy is a conserved property. A bit of internet research has not answered this question directly. Anyone have a thought on whether enthalpy is conserved in this situation (or ever)?

The first law of thermodynamics tells us that the total internal energy of the combined system is constant.

Solving for Volume B is the primary roadblock in my search for an analytical solution. Any help in this area would be appreciated.

Moving to a numerical solution.
As above, we have relationships for $\dot u_A=\dot m_A\left(h_A-u_A\right)$ and $\dot u_B=\dot m_B\left(h_A-u_B\right)$ (where $\dot m_A=-\dot m_B$)
By setting $\dot m$ constant, approximating differentials with $\frac {du} {dt}=\frac {\Delta u} {\Delta t} = \frac {u_2 - u_1} {\Delta t}$ (and similar for m) we get the following equations that relate the state in time step 2 to the state in time step 1 (assume that the state in time step 1 is fully defined):
Volume A:
$$u_2 = \Delta t \dot m\left( h_1-u_1\right) + u_1 $$
$$m_2 = \dot m \Delta t + m_1$$$$T_2 = \frac {u_2} {C_V} $$
$$P_2 = \frac {m_2RT_2} V $$
$$h_2 = C_P T_2 $$
Similarly for Volume B starting with:
$$u_2 = \Delta t \dot m \left( h_A1 - u_B1 \right) = u_1$$
... and omitting the other relations for B to preserve a shred of brevity for this post.

I put together a simulation using Euler's method to model the change in pressure over time (or, over $\Delta m = \dot m\Delta t$). The simulation can be found in this google spreadsheet. The chart below shows the change in pressure for each volume. ($\dot m$ is set to a constant, not as a function of pressure, so mass continues to move from A to B, even after the pressure in volume A is less than the pressure in volume B). It also shows the analytical solution to the pressure in volume A. Much to my disappointment, the analytical solution differs significantly from the numerical solution, suggesting that I have made an error somewhere above.
View attachment 342255

Is this simply an error of mathematics? Or have I conceptualized this problem incorrectly? I'll continue my work to find out, but any help from the community here would be appreciated.

What you are looking at is numerical inaccuracy. Either use smaller time steps or use a higher order integration scheme than forward Euler.

 

[calzone41]

The first law of thermodynamics tells us that the total internal energy of the combined system is constant.

What you are looking at is numerical inaccuracy. Either use smaller time steps or use a higher order integration scheme than forward Euler.

Total internal energy being $U_A + U_B$ in this case, right?

I re-ran the simulation with a variety of time steps. (The 0.5s stair step series is because I sample the data set at .1 seconds but the simulation only ran at .5 seconds.) The pressure values all quite agree (they essentially overlap!) relative to the analytical solution.

This indicates I've made a mistake somewhere. I'll re-check my math to see if I can find out where.

 

[Chestermiller]

Total internal energy being $U_A + U_B$ in this case, right?

Yes. There is no work done and no heat transferred.

I re-ran the simulation with a variety of time steps. (The 0.5s stair step series is because I sample the data set at .1 seconds but the simulation only ran at .5 seconds.) The pressure values all quite agree (they essentially overlap!) relative to the analytical solution.

View attachment 342302
This indicates I've made a mistake somewhere. I'll re-check my math to see if I can find out where.

 

[calzone41]

I found the error in my numerical solution above. I messed up my algebra and accidentally left out a term. Rather than: $u_2 = \Delta t \dot m\left( h_1-u_1\right) + u_1$ (for volume A), it should be $u_2 = \Delta t \frac {\dot m}{m}\left( h_1-u_1\right) + u_1$

This causes the numerical solution to agree very well with the analytical solution, even for relatively large time steps.

Interestingly, when solving the analytical solution for Volume B using the conservation of energy equation, the equation is quite messy, especially when setting $P_{A2} = P_{B2}$ and solving for $m_{A2}$ as follows:
$$P_{B2} = \frac {U_TR} {C_VV_B} - \frac {V_A}{V_B}P_{A1}\left(\frac{m_{A2}}{m_{A1}}\right)^k$$
$$m_{A2} = m_{A1}\left(\frac{RU_T}{P_{A1}C_VV_T}\right)^\left(\frac 1 k\right)$$
However, I noticed that the final pressure can be calculated as:
$$P_2=\frac {P_{A1}V_A + P_{B1}V_B} {V_A + V_B}$$
which, in retrospect, can be calculated from the conservation of energy equation $U_{A1} + U_{B1} = U_{A2} + U_{B2}$ by using the simplifying assumption $u=C_VT$ and the ideal gas law $\frac {PV} R = mT$.
I'm still looking for a similarly elegant way to calculate final mass or temperature (though, several posts above I shared the analytical solution for pressure in A as a function of mass, which is likely the best method. If anyone knows a simpler way to m or T that avoids the calculus I did above, please let me know.

I took quite a round-about path to get to an answer with this mini-project, but learned a lot on the way. This community has been helpful, providing gentle direction where needed. Thanks to you all!

My next step in this project is to analyze the three chamber problem as alluded to in the original post. I think I have learned enough to tackle it on my own, but the overall project (discrete time step simulation of a Stirling engine) will likely put up additional roadblocks that I'll seek to overcome here.

 

[Chestermiller]

I found the error in my numerical solution above. I messed up my algebra and accidentally left out a term. Rather than: $u_2 = \Delta t \dot m\left( h_1-u_1\right) + u_1$ (for volume A), it should be $u_2 = \Delta t \frac {\dot m}{m}\left( h_1-u_1\right) + u_1$

This causes the numerical solution to agree very well with the analytical solution, even for relatively large time steps.

Interestingly, when solving the analytical solution for Volume B using the conservation of energy equation, the equation is quite messy, especially when setting $P_{A2} = P_{B2}$ and solving for $m_{A2}$ as follows:
$$P_{B2} = \frac {U_TR} {C_VV_B} - \frac {V_A}{V_B}P_{A1}\left(\frac{m_{A2}}{m_{A1}}\right)^k$$
$$m_{A2} = m_{A1}\left(\frac{RU_T}{P_{A1}C_VV_T}\right)^\left(\frac 1 k\right)$$
However, I noticed that the final pressure can be calculated as:
$$P_2=\frac {P_{A1}V_A + P_{B1}V_B} {V_A + V_B}$$
which, in retrospect, can be calculated from the conservation of energy equation $U_{A1} + U_{B1} = U_{A2} + U_{B2}$ by using the simplifying assumption $u=C_VT$ and the ideal gas law $\frac {PV} R = mT$.
I'm still looking for a similarly elegant way to calculate final mass or temperature (though, several posts above I shared the analytical solution for pressure in A as a function of mass, which is likely the best method. If anyone knows a simpler way to m or T that avoids the calculus I did above, please let me know.

I took quite a round-about path to get to an answer with this mini-project, but learned a lot on the way. This community has been helpful, providing gentle direction where needed. Thanks to you all!

My next step in this project is to analyze the three chamber problem as alluded to in the original post. I think I have learned enough to tackle it on my own, but the overall project (discrete time step simulation of a Stirling engine) will likely put up additional roadblocks that I'll seek to overcome here.

Here's a hint. If you know P2, then you know $m_{A2}$ because $$P_A=Cm_A^{\gamma}$$

 

[calzone41]

Thank you. Yes, this can be found from the indefinite integral of the equation I listed a few posts ago (where I took the definite integral):
$$\frac {dT} {T} = \frac {R} {C_V} \frac {dm} m $$
Does this relationship have a name? It resembles the polytropic relationships.

 

[Chestermiller]

Thank you. Yes, this can be found from the indefinite integral of the equation I listed a few posts ago (where I took the definite integral):
$$\frac {dT} {T} = \frac {R} {C_V} \frac {dm} m $$
Does this relationship have a name? It resembles the polytropic relationships.

It is one of the equations characterizing an adiabatic reversible expansion.