A small step: Modelling tanks

Well, not that kind of tanks!

Today, we shall talk about the modelling of nice, civilized tanks - the kind of which makes it possible to have a nice warm bath every evening (boilers). The other kind of tanks might creep in here in the future, so keep your eyes open.

The problem

The task for today is to model the dynamic behaviour of a boiler of sorts (see Fig. 1). The tank has two inlets controlled by valves Z1 and Z2, through which the fluid (let’s say it’s water) of two different temperatures flows in. The valve Z3 controls the outlet. The fluid is heated with an electric heater. We shall find the system of differential equations that describes the evolution of water level h(t) and the temperature T(t) inside the tank.

We make several assumptions. First, that the tank is perfectly isolated. That is, no heat loss to the outside occurs, and no matter how long we wait, the water is not going to get any colder just by itself.

We also assume that the two incoming streams of water mix instantaneously and perfectly, so that the temperature of the fluid inside the tank is uniform. Similarly, the heat from the heater propagates instantaneously and uniformly across the tank. (And did I forget anything?)

Fig. 1. Model of a tank with two inlets, an outlet and a heater (courtesy of prof. Zbigniew Kulesza).

Parameters

What parameters are known? Let us assume the following information is available:

parameter	description
$D$	tank diameter (in $m$ )
$h_{max}$	top level of the fluid (in $m$ )
$\rho$	fluid density (in $kg/m^3$ )
$C_p$	heat capacity of the fluid (in $\frac{J}{kg \: \cdot \: K}$ )
$T_{max}$	maximum temperature of the fluid (in $K$ )
$p_1$	pressure at the inlet #1 (in $Pa$ )
$\mu_1$	valve #1 flow coefficient (dimensionless)
$k_1$	valve #1 actuator input to opening area coefficient (in $m^2/m$ )
$T_1$	temperature at the inlet #1 (in $^\circ C$ )
$p_2$	pressure at the inlet #2 (in $Pa$ )
$\mu_2$	valve #2 flow coefficient (dimensionless)
$k_2$	valve #2 actuator input to opening area coefficient (in $m^2/m$ )
$T_2$	temperature at the inlet #2 (in $^\circ C$ )
$\mu_3$	valve #3 flow coefficient (dimensionless)
$k_3$	valve #3 actuator input to opening area coefficient (in $m^2/m$ )
$U$	voltage at the heater (in $V$ )
$R$	resistance of the heater (in $\Omega$ )
$p_a$	atmospheric pressure (in $Pa$ )
$g$	gravitational acceleration (in $m/s^2$ )

Balance equations

We start with the balance equations, the general form of which is:

Input - Output = Accumulation

In our tank we accumulate two quantities: mass and the heat. Therefore it would be convenient to treat the tank as two separate accumulators, one for each of the quantities. Of course, the behaviour of those will be definitely entangled.

Let’s see… As for the mass of the fluid inside the tank the inputs are the mass flows q1, and q2, while the output is q3. The two former deposit the fluid into the tank, and the latter removes it. Whatever does not flow down the drain, remains in the tank increasing the fluid level.

As for the heat: the two input flows q1 and q2 carry some heat with them, depending on the rates of flow and the temperatures of incoming fluid. The heater adds some energy W to the tank as well. Some of the heat escapes with the outflowing fluid in the q3 stream.

We can write the basic form of the balance equations as follows:

$\frac{dm}{dt} \; = \; q_1 \; + \; q_2 \; - \; q_3$

$\frac{dQ}{dt} \; = \; W \; + \; dQ_1 \; + \; sQ_2 - \; dQ_3$

That is not quite what we want yet. For one thing, we want to find the equations describing the fluid level h(t) and the temperature T(t), which do not come out explicitly yet. Furthermore, we have to calculate all those pesky mass and heat flows, don’t we?

Let’s start with the first problem. We know that both h and T have to come somewhere into the equations for mass and for the heat. It is quite obvious for mass:

$m \; = \; \ro \: V \; = \; \rho \: A \: h$

where the area is of course $A \; = \; \frac{pi \: D^2}{4}$ , and hence h can be calculated as:

$h \; = \; \frac{m}{\rho \: A}$

Therefore, the rate of change of level can be computed quite easily as:

$\frac{dh}{dt} \; = \; \frac{1}{A \: \rho} \: \frac{dm}{dt} \; = \;\frac{1}{A \: \rho} \: (q_1 \; + \; q_2 \; - \; q_3)$

The temperature is a little more challenging. The heat accumulated in the tank is simply the temperature times the mass times the heat capacity of the fluid (we can already substitute mass here as calculated above):

$Q \; = \; m \: C_p \: T \; = \; A \: \rho \: C_p \: h \: T$

Two of the variables on the right side (h and T) depend on time. In calculating the time derivative of Q we have to take this into account (the equation for time derivative of a product of functions is $(fg)' \: = \: fg' \: + \: gf'$ ):

$\frac{dQ}{dt} \; = \; A \: \rho \: C_p \: \frac{d(hT)}{dt} \; = \; A \: \rho \: C_p \: (h \: \frac{dT}{dt} \: + \: T \: \frac{dh}{dt})$

The time derivative of temperature is then:

$\frac{dT}{dt} \; = \; \frac{1}{A \: \rho \: C_p \: h} \: \frac{dQ}{dt} \; - \; \frac{T}{h} \: \frac{dh}{dt}$

We now know all the things that go on the right side of this equation. As we’ll see later on, it is furthermore going to simplify a bit more.

Valve model

Let’s now calculate the flows of the water and the heat. How does a control valve work? The actuator operates linearly and moves the “mushroom” (or other kind of a shutter) inside the valve. The shutter’s movement creates an opening of a certain area, which in the simplest case is proportional to the product of the actuator’s input and a coefficient k. A simple model for mass flow through a valve is given as:

$q \; = \; \mu \: k \: z \: \sqrt{2 \: \rho \: \Delta p}$

where $q$ is the mass flow in [kg/s], $\mu$ is the mass flow coefficient of the valve (dimensionless), $k$ is the rate of the valve opening area to the $z$ - the actuator input. The density of the fluid $\rho$ and the pressure difference $\Delta p$ also come into that equation.

Calculating flows

Armed with the knowledge about the valves, we can proceed with our calculations. Mass flows are now simple to compute:

$q_1 \; = \; \mu_1 \: k_1 \: z_1 \: \sqrt{2 \: \rho \: (p_1 \: - \: p_a)}$
$q_2 \; = \; \mu_2 \: k_2 \: z_2 \: \sqrt{2 \: \rho \: (p_2 \: - \: p_a)}$

The third mass flow q3 is somewhat of a special case. The pressure at the valve is a hydrostatic one, cause by the height of the water level. The atmospheric pressure plays no part here, since it’s present on both sides of the valve:

$q_3 \; = \; \mu_3 \: k_3 \: z_3 \: \sqrt{2 \: \rho^2 \: g \: h}$

The heat streams carry the heat proportional to the temperatures of fluid and the mass flow:

$dQ_1 \; = \; C_p \: T_1 \: q_1$
$dQ_2 \; = \; C_p \: T_2 \: q_2$
$dQ_3 \; = \; C_p \: T \: q_3$

The heat produced by the heater is:

$W \; = \; \frac{U^2}{R}$

Arriving at the final equations

It is now a matter of putting things together. Recall the equations we’ve came up with from the mass and heat balance:

$\frac{dh}{dt} \; = \; \frac{1}{A \: \rho} \: \frac{dm}{dt} \; = \;\frac{1}{A \: \rho} \: (q_1 \; + \; q_2 \; - \; q_3)$
$\frac{dT}{dt} \; = \; \frac{1}{A \: \rho \: C_p \: h} \: \frac{dQ}{dt} \; - \; \frac{T}{h} \: \frac{dh}{dt}$

I’m not going to write down all the steps here. Plugging in and simplifying, we get:

$\frac{dh}{dt} \; = \; \frac{1}{A \: \rho} \: \frac{dm}{dt} \; = \;\frac{1}{A \: \rho} \: (q_1 \; + \; q_2 \; - \; q_3)$
$\frac{dT}{dt} \; = \; \frac{q_1 \: \: (T_1 \: - \: T) \: + \: q_2 \: (T_2 \: - \: T) \: + \: \frac{U^2}{R \: C_p}}{A \: \rho \: h}$

Our system can be noted down in the following form:

$\mathbf{x} \; = \; [ \: h \; T \: ]^T$
$\mathbf{u} \; = \; [ \: z1 \; z2 \; z3 \; U \: ]^T$
$\mathbf{y} \; = \; [ \: h \; T \: ]^T$

We have two state variables, four inputs and two outputs.

Simulation in Matlab

The simulation of the tank model can be easily done in Matlab. We need two scripts. The master script sets up the stage: defines parameters and the simulation range. It is simple enough:

clc; clear all; close all;

%% Data

% Tank & fluid properties:
D = 0.1;        % tank diameter [m]
ro = 1000;      % fluid density [kg/m^3]
Cp = 4190;      % specific heat capacity [J/(kg*K)]

% Input & output properties:
p1 = 200000;    % pressure at valve #1 [Pa]
mu1 = 0.05;     % flow rate for valve #1 [1]
k1 = 0.001;     % coefficient of valve opening area to actuator position [m^2/m] 
T1 = 10;        % temperature of fluid at valve #1 [deg C]

p2 = 200000;    % pressure at valve #2 [Pa]
mu2 = 0.05;     % flow rate for valve #2 [1]
k2 = 0.001;     % coefficient of valve opening area to actuator position [m^2/m]
T2 = 60;        % temperature of fluid at valve #2 [deg C]

mu3 = 0.05;     % flow rate for valve #3 [1]
k3 = 0.01;      % coefficient of valve opening area to actuator position [m^2/m]

U = 230;        % voltage for the heater [V]
R = 0.0001;     % resistance of the heater [ohm]

% General properties:
pa = 100000;    % atmospheric pressure [Pa]
g = 9.81;       % gravitational acceleration [m/s^2]

% Simulation time:
t = [0 : 0.1 : 60]; % [s]

%% Initial conditions
x0 = [0.5 20];  % [ h T ]

%% Solution
options = odeset('AbsTol', 1e-3, 'RelTol', 1e-3);
[t, x] = ode23s(@tank_odefun, t, x0, options, [D ro Cp p1 mu1 k1 T1 p2 mu2 k2 T2 mu3 k3 U R pa g]);

%% Plot results
figure;
plotyy(t, x(:, 1), t, x(:, 2));
legend('h [m]', 'T [deg C]');
xlabel('t [s]');
ylabel('h [m], T [deg C]');

The tank_odefun function calculates the derivatives of the state variables at any time point. We arrived at these equations in previous sections - it is now only a matter of writing them in a form understandable by Matlab.

function [ dx ] = tank_odefun( t, x, p )
%TANK_ODEFUN Calculates derivatives for the tank example

%% Parameters
D = p(1);
ro = p(2);
Cp = p(3);
p1 = p(4);
mu1 = p(5);
k1 = p(6);
T1 = p(7);
p2 = p(8);
mu2 = p(9);
k2 = p(10);
T2 = p(11);
mu3 = p(12); 
k3 = p(13);
U = p(14);
R = p(15);
pa = p(16);
g = p(17);

%% Inputs & state variables
z1 = 1 * heaviside(t-0);
z2 = 1 * heaviside(t-10) - 1 * heaviside(t-30);
z3 = 1 * heaviside(t-20);
U  = 1 * heaviside(t-0);

h = x(1);
T = x(2);

%% Calculate helper coefficients
A = pi*D^2/4;
q1 = mu1*k1*z1*sqrt(2*ro*(p1-pa));
q2 = mu2*k2*z2*sqrt(2*ro*(p2-pa));
q3 = (h>0)*mu3*k3*z3*ro*sqrt(2*g*h);

%% Calculate derivatives
dx(1, 1) = 1/(A*ro) * (q1 + q2 - q3);
dx(2, 1) = 1/(A*ro*h) * (q1*(T1-T) + q2*(T2-T) + U^2/(R*Cp));

end

The control of the inputs is also encoded in the tank_odefun function. It’s done by employing the heaviside functions, which are simply unit steps enabled at given times. An example simulation of the tank model is shown in Fig. 2. In this example, the Z1 valve is always turned on, the Z2 valve is opened at 10 s mark (and shut at 20 s), and the Z3 is opened after 20 s. The heater is always on.

Fig. 2. Matlab simulation of the tank model.

Interactive simulator

I’ve also made an interactive Javascript implementation. It’s a nice toy to play with, and all the inputs, outputs and parameters are plainly available. The human mind is after all a visual and tinkering processor, and it’s best to lay your hands on a thing you want to work out.

The simulator is available here: SIMULATOR

The simulator layout is presented in Fig. 3. The panel on the left is used to tweak the values of parameters (it is only possible while the simulation is not running). The controls are on the right: you can open and close any valve (in range 0 - 1), and turn the heater on and off as well. The state and the outputs are shown in the bottom right corner.

Some additional quirks are present in this model. When the water reaches the top of the tank, it overflows, such that it cannot rise arbitrarily high. A top bound on the temperature is also implemented (water cannot be heated above 100 degrees).

The state equations are solved numerically using a Runge-Kutta fourth order solver at a step $dt \: = \: 0.01 \: s$ . (I shall write about numerical integration some other time.)

Fig. 3. Tank model simulator.

Try to play with the model: you might for example want to keep the constant water level and the temperature while the output changes. In a real industrial scenario, you would use a controller (perhaps a PID). See how difficult is to control the plant by hand.

Controllers - they work, so we don’t have to.

(Source code - both JS and Matlab is available here: https://github.com/dagothar/modelling-tanks)

Do not hesitate to ask or point out bugs!

2017-04-05

Modelling tanks