1. Modelling a Control Volume in Matlab
    1. Preliminaries
    2. Upwind scheme
    3. Boundary and Initial Conditions
    4. The equations
      1. Mass Balance
      2. Energy Balance
    5. Implementing the code in Matlab
    6. References

Modelling a Control Volume in Matlab

In this post, I want to discuss how we can develop an Minimum Viable Model for analysing a control volume in Matlab. This makes for a good example because it can cover several areas of thermofluids: the idea of the upwind scheme, the solving of ODEs, the calculation of partial derivatives and so on.

We will be modelling a rigid pipe with two ports; an inlet port and an outlet port. In addition, we assume that the wall is perfectly adiabatic, that is, no heat exchange occurs between the pipe and the surroundings.

Note that if you want to follow along, you will need to have a copy of Matlab (or maybe Octave will do) and the NIST software Refprop (REFrigerant PROPerties).

Now, where do we begin? It might be helpful to go back and read my article to see the how we develop a general form of these equations. Here, we will start straight with the implementation of these equations.


Ok, so where do we begin? We remember that, for any control volume, we need to know three things: two independent thermodynamic variables, and a flow variable. In the field of thermofluids, we often use Pressure and Specific Enthalpy as the state variables. Why? Well, for one thing, we rely on the pressure-specific enthalpy diagram often to try to make sense of thermodynamic cycles. And a second thing is that, in the two-phase (the 'dome' in the P-h diagram below), knowing the enthalpy tells us how far along we are into the two-phase, aka, how much of the fluid we have evaporated/condenesed.

As for the flow variable, the mass flow rate is often selected, because it is more of an 'objective' thing, in that it is unaffected by the cross section it is flowing through. For steady flow, the mass flow rate will be the same everywhere, in a way that the volumetric flow rate (L/s) or the velocity (m/s) won't be.

So, we now have a 'lumped' volume, and we need to know three things about it: the pressure, the specific enthalpy and the mass flow rate.

Upwind scheme

Now, let's talk about one of my most favourite mathematical techniques ever. The idea of the upwind scheme. We use this technique to be able to cohesively determine what part of the model calculates what. For the current model, we need the technique to determine what are the boundary conditions, vs what is calculated by the model itself.

We could, if we wanted, solve all three equations inside the same control volume. The trouble is this: the three governing equations all are deeply interconnected. Changes in one equation will have impacts on the other equations. We call this having the equations 'coupled'. This is a nightmare situation for the differential equation solvers because it brings them to a crawl.

To avoid this, Courant et al[1], back in 1952, proposed the upwind scheme, which was later covered by Patankar[2] in his brilliant 1980 book. Using this method, we 'decouple' the equations from each other. We solve the mass and energy balance on one grid called the 'thermal grid', and, offset by half a cell width, we solve the momentum balance on the 'momentum grid'.

The cells are the same width, and thus, offsetting by half a cell-width means that the centers of one grid coincide with the boundaries of the other grid. This way, whatever is calculated inside one grid is automatically averaged right at the boundary of the neighbouring grid, and thus no interpolation is necessary.

This concept is pretty tricky to understand, but ultimately, it's sufficient if you remember that, if we're dealing with a capacitive component, the mass flow rate is determined by the upwind scheme, and, if we're dealing with a resistive component, the pressure is determined by the upwind component. The incoming enthalpy, meanwhile, is ALWAYS determined by upstream conditions.

Boundary and Initial Conditions

So, having discussed the upwind scheme, we now take the control volume itself to be a 'capacitive' element. That is, the control volume solves for the pressure and specific enthalpy using the mass and energy balance equations. Then, the inlet and outlet ports must supply the 'resistive' element of the equation, aka the mass flow rate. Then, just before the inlet port would be some momentum grid element, and likewise for just beyond the outlet port, were we to build up a complete system model.

Note, however, that in addition to the mass flow rate, we also need to supply the enthalpy of the incoming fluid, as the upwind scheme dictates that the inlet enthalpy be defined by the component upstream of it. Thus, the system setup ends up looking like:

In addition, when solving differential equations that involve time, we need to provide initial conditions. We need as many initial conditions as there are dynamic variables, two in this case (pressure and enthalpy).

The equations

Now, let's talk about the governing equations. First, the mass balance.

Mass Balance

Our flow is one-dimensional. Thus, our equation will take the form:

\[\frac{\partial{(\rho V)}}{\partial t} = \dot{m}_{in}-\dot{m}_{out}\]

Rather straightforward, right? The equation basically states that any refrigerant flowing into the control volume can do one of two things: 1) increase the mass inside the volume, or 2) flow out of the volume.

Notice, however, that our volume is not deformable. The control volume is rigid, and thus the derivative of Volume with time will be zero. The equation then becomes:

\[V\frac{\partial{(\rho)}}{\partial t} = \dot{m}_{in}-\dot{m}_{out}\]

Next, we notice that the equation above is in terms of density, and not in terms of pressure and specific enthalpy, which are our state variables. So what do we do?

We first remember that in thermodynamics, any variable can be found using two other known variables. Then, \(\rho=f(P,h)\)

If we apply the chain rule of differentiation to this fact, we get

\[V\frac{\partial{\rho}}{\partial t} = \frac{\partial\rho}{\partial P}\frac{dP}{dt}\Bigg|_h + \frac{\partial\rho}{\partial h}\frac{dh}{dt}\Bigg|_P\]

Here, the pipe symbols with the subscripts mean that when we are taking the partial derivative of the pressure, we keep enthalpy constant (and vice versa).

Substituting this into the mass balance equation gives:

\[V\left(\frac{\partial\rho}{\partial P}\frac{dP}{dt}\Bigg|_h + \frac{\partial\rho}{\partial h}\frac{dh}{dt}\Bigg|_P\right) = \dot{m}_{in}-\dot{m}_{out}\]

These partial derivatives will be calculated using Refprop and some fluid property magic in the two-phase region. But now let's do the same for the energy balance.

Energy Balance

The energy balance equation will look like this:

\[\frac{\partial{(\rho uV)}}{\partial t} = \dot{m}h_{in}-\dot{m}h_{out}\]

This equation is along the same lines as the previous ones, stating that if we have a flow coming in (carrying energy), it can either flow out, or raise the internal energy of the fluid inside.

Now, as before, we note that the volume is constant. So, we are left with \(\rho.u\) inside the derivative. Then, we also know that the enthalpy is \(h = u + P/\rho\). Thus, \(u = h - P/\rho\). If we substitute this into the equation above and run the crank of deriving it out, we will get:

\[V\left(\left(h\frac{\partial\rho}{\partial P}-1\right)\frac{dP}{dt}\Bigg|_h + \left(h\frac{\partial\rho}{\partial h}+\rho\right)\frac{dh}{dt}\Bigg|_P\right) = \dot{m}h_{in}-\dot{m}h_{out}\]

Implementing the code in Matlab

Right! Onto the good stuff. How do we implement this in matlab? Let's first try to understand how ODEs are set up in Matlab. Well, ODEs are set up in the form of matrix equations in Matlab. They take the form Ax = B, with x being the vector of state variables, A is the matrix of the partial derivatives, while B contains the constraints, or the boundary conditions in our example.

Then, if we are able to set up a function that returns the values of these partial derivatives, we can send this to our ODE function, and it can integrate the equation over time. Thus, we need a function that outputs x = A\B

Let's look at that code now.

function y = controlVolume(~,x)
% controlVolume Lumped control volume
% Inlet boundary conditions: mass flow rate and enthalpy
% Outlet boundary condns: mDot
% State variables: Pressure (Pa) and enthalpy (J/kg)

global fld Vol mDotIn hin mDotOut

P = x(1); % Pa
h = x(2); % J/kg

[drho_dPh,drho_dhP] = getPartialDers(P,h,fld);
rho = refpropm('D','P',P*1e-3,'H',h,fld);

Hin = mDotIn*hin;
Hout = mDotOut*h;

a = Vol*drho_dPh;
b = Vol*drho_dhP;    
c = Vol*(h*drho_dPh-1);
d = Vol*(h*drho_dhP+rho);

M = [a,b ; c,d];
f = [mDotIn-mDotOut ; Hin-Hout];

y = M\f;

This code should hopefully look pretty understandable to you by now. The matrices A and B are here called M and f to be more consistent with the implicit equation nomenclature that Matlab uses. Don't worry about it. It's just another name.

So, what is this function spitting out? The values of dy/dt, aka [dP/dt ; dh/dt]. The getPartialDerivatives function calculates the partial derivatives, mindful of where we are in the Pressure-enthalpy diagram (single phase or two-phase).

% Pressure (Pa) and specific enthalpy (J/kg) state variables
% Inlet and outlet mass flow rates and enthalpies boundary conditions
% (Note: Pressure in Pa for state var, bar for plots and kPa for refprop)

close all

global fld
fld = 'CO2';

% Time discretization
t = 60;
dt = 0.01;
tspan = 0:dt:t;
nt = length(tspan);

% Geometry
global Vol
D = 0.2; % m
L = 0.32; % m
Vol = pi * (D/2)^2 * L; % m

% Boundary Conditions
global mDotIn hin mDotOut
mDotIn = 0.1;
hin = 100000;
mDotOut = 0.1;

% Initial Conditions
P0 = 20e5; % Pa
h0 = 100000; % J/kg
x0 = [P0,h0]; %

% Solving ODEs
optionSet = odeset('OutputFcn',@odeplot,'RelTol',1e-5);
[t,x] = ode15s(@controlVolume, tspan, x0, optionSet);

% Plotting
P = x(:,1).*1e-5;
h = x(:,2).*1e-3;
plotPh(); % my own function for plotting the pressure-enthalpy diagram

Right, so a couple of things to comment here. First the simulation parameters. The way I set up the simulation above is to have a 10 L vessel, full of a two-phase mixture of CO2. Then, we have 1 kg/s flowing in, and 1 kg/s flowing right out again. The incoming fluid is at very cold enthalpies of 100 kJ/kg (corresponding to < -45 C at 20 bar).

So, what would you expect? Well, since the flow in and flow out is the same, the mass should stay constant. And indeed, that is what we see. Look at the P-h diagram below

The green lines are the isochores, aka the constant density lines. The red line, as you can see, proceed along the constant density line. Slowly, the control volume is also cooling down. All seems to make sense.

So, what have we achieved? Well, if you were able to follow along, you should now be able to create what is the basic component for all thermofluid simulations. By creating successions of such control volumes, we can model very complex systems indeed, and that is something I hope to continue talking about in future posts.


[1] Courant, R., Isaacson, E., & Rees, M. (1952). On the solution of nonlinear hyperbolic differential equations by finite differences. Communications in Applied Numerical Methods, 5(3), 243–255.
[2] Patankar, S. (1980). Numerical heat transfer and fluid flow. Series in Coputational Methods in Mechanics and Thermal Sciences.