# Unsure Regarding Temperature Model for an Open System

Tags:
1. Nov 13, 2017

### John Kenneth

Hi guys,

I am currently working on my master thesis. I am supposed to make a dynamic model for a gas system and have some trouble setting up the energy balance. I am a noob when it comes to uploading pictures, so I dont have a figure for this, but lets consider a general control volume (CV), with one flow in and one flow out.
I am interested in deriving an expression for the change in temperature, dT/dt.

I have attempted to solve the dynamic energy balance with respect to temperature, but this ends up in a pretty nasty expression. I have read a lot on thermodynamics lately to solve this, but I am not sure if my final expression, or approach, is correct. I would be extremely glad if someone with expertise on thermodynamics could look into this and see if it makes sense, or if there are any obvious mistakes or misconceptions!

Okay, so this is my approach for solving the energy balance for an open system:

The dynamic energy balance is given as [1]:

\boxed{
\frac{dU}{dt} = H_{in} - H_{out} + \delta Q - \delta W
}

where:
U = internal energy
W = work
H = enthalpy

Hin and Hout are the enthalpy of the flow entering and leaving the control volume, and are defined as

\begin{split}
&H_{in} = \dot{m}_{in} h_{in} \\
&H_{out} = \dot{m}_{out} h_{out}
\end{split}

where the lower case h is the specific enthalpy, h=H/m.

Introducing U = H - pV, the energy balance can be expressed with enthalpy as:

\frac{dH}{dt} = H_{in} - H_{out} + \delta Q - \delta W + pdV + Vdp

where:
p = pressure
V = volume
The last two terms, p dV + V dp, are derived from the chain rule for d(pV).
Considering the varying mass, dH can be written as

\frac{dH}{dt} = m \frac{dh}{dt} + h \frac{dm}{dt} \Rightarrow m\frac{dh}{dt} = \frac{dH}{dt} - h\frac{dm}{dt}

The mass balance for the control volume, dm, is defined as:

\frac{dm}{dt} = \dot{m}_{in} - \dot{m}_{out}

With some rearrangement, the energy balance in enthalpy form becomes

m_{cv} dh = \dot{m}_{in} h_{in} - \dot{m}_{out} h_{out} + \delta Q - \delta W + pdV + Vdp - h ( \dot{m}_{in} - \dot{m}_{out})

From what I understand, the last enthalpy term, h, is the enthalpy currently in the control volume, and not the enthalpy "following" the flow in and out.

Combining the enthalpy terms gives:

\boxed{
m_{cv} \frac{dh}{dt} = \dot{m}_{in}(h_{in} - h_{cv}) - \dot{m}_{out}(h_{out} - h_{cv}) + \delta Q - \delta W + pdV + Vdp
}

where mcv is the mass in the control volume.

Now, I will replace the enthalpy and dp terms.

In my project, I am using the Redlich-Kwong equation of state (EOS). This is given as:

p = \frac{RT}{V_m-b} - \frac{a}{\sqrt{T} V_m (V_m+b)}

where:
R = specific gas constant
T = temperature
Vm = specific volume, also defined as ρ-1
a = given constant for the EOS
b = given constant for the EOS

In the Redlich-Kwong EOS, the pressure is given as a function of temperature and specific volume. Therefore, the time derivative of the pressure will be

\frac{dp}{dt} = \left(\frac{\partial p}{\partial T}\right) \frac{dT}{dt} + \left(\frac{\partial p}{\partial V_m}\right) \frac{dV_m}{dt}

The "real" enthalpy can be calculated using residual properties [2]. The residual property gives the enthalpy as a deviation from the ideal enthalpy. This means that the enthalpy can be written

h = h_{ideal} + h_{residual}

The ideal enthalpy can be written as [1]

h = \int_{T_{ref}}^{T} c_p dT

with Tref as the temperature at a set reference state.
It can be shown that the residual enthalpy is

h_{residual} = \frac{bRT}{V_m -b} - \frac{a}{\sqrt{T} V_m + b} - \frac{3a}{2 \sqrt{T} b} \ln (\frac{V_m + b}{V_m})

The dh term on the left side in the energy balance is the time derivative of the combined ideal and residual enthalpy. Derivation of these two requires the partial differential of h in regards to T and Vm :

dh = \left(\frac{\partial h}{\partial T}\right) \frac{dT}{dt} + \left(\frac{\partial h}{\partial V_m}\right) \frac{dV_m}{dt}

This gives

dh = \left[c_p + \frac{bR}{V_m-b} + \frac{a}{2 T^{\frac{3}{2}}(V_m+b)} + \frac{3a}{4T^{\frac{3}{2}}b} \ln{\frac{V_m+b}{V_m}}\right] \frac{dT}{dt} + \left[ \frac{bRT}{(V_m-b)^{2}} + \frac{a}{\sqrt{T} (V_m+b)^{2}} + \frac{3a}{2 \sqrt{T} V_m(V_m+b)} \right] \frac{dV_m}{dt}

Now, we have to replace the enthalpy terms and the dp term. This requires a bit of algebraic manipulation which ends up in a nasty expression, but bear with me!

Inserting for dh and dp:

m_{cv} \left(\left[c_p + \frac{bR}{V_m-b} + \frac{a}{2 T^{\frac{3}{2}}(V_m+b)} + \frac{3a}{4T^{\frac{3}{2}}b} \ln{\frac{V_m+b}{V_m}}\right] \frac{dT}{dt} + \left[ \frac{bRT}{(V_m-b)^{2}} + \frac{a}{\sqrt{T} (V_m+b)^{2}} + \frac{3a}{2 \sqrt{T} V_m(V_m+b)} \right] \frac{dV_m}{dt} \right) = \dot{m}_{in}(h_{in} - h_{cv}) - \dot{m}_{out}(h_{out} - h_{cv}) + \delta Q - \delta W + pdV + V \left(\frac{\partial p}{\partial T}\right) \frac{dT}{dt} + V \left(\frac{\partial p}{\partial V_m}\right) \frac{dV_m}{dt}

Combine dT terms on left side:

\left[ m_{cv} \left(c_p + \frac{bR}{V_m-b} + \frac{a}{2 T^{\frac{3}{2}}(V_m+b)} + \frac{3a}{4T^{\frac{3}{2}}b} \ln{\frac{V_m+b}{V_m}}\right) - V\frac{\partial p}{dT} \right] \frac{dT}{dt} = \dot{m}_{in}(h_{in} - h_{cv}) - \dot{m}_{out}(h_{out} - h_{cv}) + \delta Q - \delta W + pdV + m_{cv} \left( V \left(\frac{\partial p}{\partial V_m}\right) \frac{1}{m_{cv}} - \frac{bRT}{(V_m-b)^{2}} + \frac{a}{\sqrt{T} (V_m+b)^{2}} + \frac{3a}{2 \sqrt{T} V_m(V_m+b)} \right) \frac{dV_m}{dt}

Dividing on the left side to get dT alone:

\boxed{
\frac{dT}{dt} = \frac{ \left[ \dot{m}_{in}(h_{in} - h_{cv}) - \dot{m}_{out}(h_{out} - h_{cv}) + \delta Q - \delta W + pdV + m_{cv} \left( V \left(\frac{\partial p}{\partial V_m}\right) \frac{1}{m_{cv}} - \frac{bRT}{(V_m-b)^{2}} + \frac{a}{\sqrt{T} (V_m+b)^{2}} + \frac{3a}{2 \sqrt{T} V_m(V_m+b)} \right) \frac{dV_m}{dt} \right]}{ \left[ m_{cv} \left(c_p + \frac{bR}{V_m-b} + \frac{a}{2 T^{\frac{3}{2}}(V_m+b)} + \frac{3a}{4T^{\frac{3}{2}}b} \ln{\frac{V_m+b}{V_m}}\right) - V\frac{\partial p}{dT} \right]}
}

So... This is my final expression for the temperature. I hope it was possible to follow this somewhat cumbersome derivation. If not, do not hesitate to ask for more details and I will try to elaborate. As i stated in the introduction, I would really appreciate if someone could take a look at this, as it would be of great help for my thesis.

Sources:
[1]. SKOGESTAD, Sigurd. Chemical and energy process engineering. CRC press, 2008.
[2]. SANDIP, Roy. Nptel.ac.in. NPTEL :: Chemical Engineering - Chemical Engineering Thermodynamics. (2017). [online] Available at: http://nptel.ac.in/courses/103101004/1 [Accessed 13 Nov. 2017].

2. Nov 14, 2017

### Nidum

You seem to be making your analysis of this problem more complicated than it needs to be .

If the problem is really as simple as you say then the analysis required is both relatively straightforward and well documented .

Is there some added complexity to this problem that you haven't told us about ?

Last edited: Nov 14, 2017
3. Nov 15, 2017

### John Kenneth

Hi,

I also think that my final expression seems a bit complicated, and I might have been a bit confused about all the different approaches available.

When it comes to added complexity, I can tell you that the gas flow will vary in direction a lot. Therefore, the model should simulate a transient/unsteady state. As this is a transient analysis one can not set

\dot{m}_{in} = \dot{m}_{out}

or

dU = 0

Hence, I need to include the mass flow balance, where flows in and out of the volume are not necessarily equal.

For an ideal gas, the enthalpy is only dependent on temperature. I am using the Redlich-Kwong equation of state (RK EOS), thus indicating "real" gas, and not ideal gas. Therefore, I need to consider the enthalpy as pressure dependent as well, giving it the extra residual term. Since the RK EOS is given in terms of temperature and specific volume, the residual enthalpy will also be a function of these variables. Taking the time derivative of hresidual includes partial derivatives of temperature and specific volumes, which leads to a more complicated expression for dh than if it were considered ideal.

I feel that the "real" enthalpy is making the equations a bit more complicated, especially when this term has to be derived with respect to time. Also, the time derivative term of dp contains partial derivatives for temperature and specific volume, which adds even more complexity. Substitutions of these two terms in equation (7) are the main components which leads to the final boxed expression.

I hope that this made things a bit clearer, and that it might help you understand my thoughts when setting up these equations. Do you think what I have done is correct, or have I missed something along the way?
Any help would be appreciated :)

4. Nov 15, 2017

### Staff: Mentor

Some questions: You are processing a non-ideal gas that is described by the RK equation. How sure are you that the gas can't be described by the ideal gas law? Why did you convert to enthalpy on the left-hand side of the equation? Why not stick with internal energy? Is shaft work really being done on this system? You are treating the system as being lumped parameter, with prefect mixing and no spatial gradients within the system?

Have you solve this for steady state operation yet (to get some practice)?

Last edited: Nov 15, 2017
5. Nov 16, 2017

### John Kenneth

Regarding the RK equation, I am supposed to model the gas via some non-ideal gas equation to investigate the differences between ideal and non-ideal gas. So therefore I have to assume non-ideal behavior.

The reason for converting to enthalpy on the left-hand side was simply because the textbook I was using (source [1]) suggested it, and stated that it was a usual way of conducting the energy balance. I also found a good source (source [2]) which explained how the residual enthalpy is calculated, and also contained a detailed derivation for the expression.

There are no shaft work being done in this system. I see now that this might lead to some confusion, and the work term W could be set to zero for the most part. One of the volumes contain a piston which will do boundary work on the system, therefore I kept the work term (a bit sloppy notation from my side). But for the general control volume considered in this post, the work term could be set to zero.

I am treating the system as being lumped parameter as far as I understand the term (mainly solving ODE's). The system contains nitrogen gas which is assumed to stay in gas form at all times, and I assume perfect mixing for the system. I also consider it to be evenly distributed in the control volumes with no spatial gradients.

I solved the temperature equation for closed system, and compared it to the new expression I got. It behaves similar to that of the closed system, but has greater amplitude/values. This is as expected due to the added energy brought to the volume via the mass flow (the comparison was tested for a tank with only one stream into the volume, thus increasing the mass). I think the results looks fine thus far, but I still have to test it for the whole system. As the derivation of the final expression is quite long, I hoped that someone could look over it and see if it still looks reasonable with the given assumptions.

6. Nov 16, 2017

### Nidum

Hard to reconcile that with a simple lumped parameter model . If directions of the flows matter then doesn't that really mean that you need to look at the detailed fluid mechanics of your system ?

This would all be easier to understand if you posted a diagram of the system that you are working on . You don't need CAD or anything complicated . Do a clear sketch by hand or using simple software like Paint and then upload or copy and paste .

7. Nov 16, 2017

### Staff: Mentor

I assume that you have solved this problem, both steady state and transient, for the case of an ideal gas. Given @Nidum's request for a more detailed description of the problem, can you tell us about your approach and results for steady state and transient for an ideal gas?

Are the perturbations to the system for the transient case typically small (so that the analysis can be linearized) or are they going to be large, such as in the startup of the system?

8. Nov 17, 2017

### John Kenneth

Chestermiller: I have not tested the whole system for steady state and transient for an ideal gas. Along the way, while building the model using matlab, I have compared the results to that of ideal gas for special cases to verify that the results are not too far off. The system is acting like a relatively big air spring, and the flow is therefore oscillating back and forth between two volumes. I therefore assume that the system will never be in a steady state long enough to make a steady state simulation accurate. But maybe it would be a good idea to drive the system slow enough to assume steady state, and then test steady state simulation for both ideal and non-ideal behavior? That seems like a good idea for verification of parts of the equations. I will also make a complete model for the whole system with ideal gas in transient state for comparison after I have set it up for the non-ideal case.

As the flow is oscillating back and forth between two volumes, I guess that the perturbations are going to be relatively large due to change in flow direction .

Nidum: The model I make will be used for real time HIL-testing, and maybe some simulator training of personnel. So the main point is to make it run fast enough for these applications, and try to simulate the real behavior as accurate as possible within these boundaries. I will try to add some non-ideal behavior, and see how efficiently the code will run, and how much it differs from similar ideal models. Optimally, one should go into the details of the fluid mechanics, and maybe even look into CFD analysis. I am afraid that if I go deep into the details, this will lead to too long simulation times for the intended applications. When it comes to modelling, I have therefore chosen one-dimensional flow based on control volume analysis. The system is split into four control volumes: One control volume for each volume (1&2), while the pipe is divided into two control volumes. The plan is that it should be possible to split the pipe into n volumes, and see how it effects the dynamics, but for now we will start with the pipe as two control volumes.
I have tried to make some easy sketches in paint to illustrate this.

This is a quick sketch of the total system. The valve can be neglected for now.

This is the two volumes. Volume 1 will act as a piston accumulator with a movable piston, while volume 2 will be an auxiliary volume without any piston or mechanical devices.

This is the gas pipe. The two control volumes will be equal and divided by the valve.

Last edited: Nov 17, 2017
9. Nov 17, 2017

### Staff: Mentor

I am not able to see the diagrams. Did you upload them?

10. Nov 17, 2017

### John Kenneth

I tried to edit the links to the photos now. Does it work now?

11. Nov 17, 2017

### Staff: Mentor

yes.

12. Nov 17, 2017

### Staff: Mentor

What is the rationale for not treating the entire combination of 2 tanks and an intervening pipe as a closed system at a uniform T and P?

13. Nov 19, 2017

### Nidum

@John_Kenneth .

Can you tell us more about the practical purpose of this system and what the actual cycle of operation is ?