Tank Heater Example

Heated, Stirred Tank As a simple process example, consider a heated, stirred tank. Liquid flows continuously in and out of the tank via streams Fin and Fout (F is a volumetric flowrate). Inside the tank, there is a liquid level h. A heating coil can be used to transfer heat to the liquid. A second inlet stream D is also present.

This basic example will be used to demonstrate a number of aspects of control modeling. Different sections will apply at different times during the course:

Some equations in this document are being displayed using MINSE, a browser independent approach to displaying equations on the web. If the equations do not look properly formatted, please RENDER them by selecting the link and invoking Ping's MINSE polymediator. This will run a special program that should cause the equations to be formatted for viewing by your browser.

Model Development (Shell Balances)

Maybe the most fundamental method for developing a differential equation model for use in control uses "shell balances". This approach begins by defining a "shell" of "differential thickness" (which makes more sense when done for position dependence than it does for time!). Maybe this makes most sense as an example.

Mass Balance

Our system will be the "water in a tank". Our model will be the total material balance.

Generic Balance Equation
The accumulation is the amount by which the total mass in the tank changes over a "differential slice of time" ?Delta?t, or
Accum = M;(t = t+?Delta?t) - M;(t= t)
The net flow of mass into the tank can be written as the product of the mass flow rate and the time interval
Flow;In = m;i*?Delta?t
and a similar term written for the mass flow out
Flow;Out = m;o*?Delta?t
Notice that both these terms have mass units. Now if the whole thing is put together
M;(t = t+?Delta?t) - M;(t= t) = m;i*?Delta?t - m;o*?Delta?t
'quot(M;(t = t+?Delta?t) - M;(t= t), ?Delta?t) = m;i - m;o
The next step is CALCULUS and we'll need the definition of a derivative (do you remember?)
Defn. of Derivative
This means that if we take the limit of both sides of our working equation as ?Delta?t .approach 0 we end up with
'deriv(M;tank, t) = m;i - m;o
This is the mass balance for the system. The accumulation term is the differential (rate of change with respect to time) of the total mass in the tank. If the mass flow entering the tank is greater than that leaving, the mass in the tank is increasing, the level is rising, and the derivative is positive. The opposite (derivative < 0) is true when the mass in the tank is decreasing.

The accumulation term in our dynamic models will always be a time derivative; the time derivative of the total amount of whatever we are balancing contained within whatever system we have defined. With some thought, and maybe a little practice, one can write the derivative for the accumulation directly and dispense with the intermediate steps. The work will be tricky enough without that piece of math.

Normally, we don't "know" the total mass in the tank. We're more likely to know the volume in the tank, and the volumetric flow rates in and out. Volume is not conserved so we cannot write a volume balance. The bridge between volume and mass is density, so we typically write our mass balance as:

'deriv(V*?rho?,t)=F;i*?rho?;i-F;o*?rho?;o
If you want to try the shell balances again, it goes like this:
(V*?rho?);(time=t+?Delta?t) - (V*?rho?);(time=t) = F;i*?rho?;i*?Delta?t - F;o*?rho?;o*?Delta?t
'lim('quot((V*?rho?);(time=t+?Delta?t) - (V*?rho?);(time=t)),?Delta?t,t .approach 0) = F;i*?rho?;i - F;o*?rho?;o
'deriv(V*?rho?,t)=F;i*?rho?;i-F;o*?rho?;o

Typically, we won't be directly measuring the volume in the tank any more than the mass; what is measured is the level, or height, of liquid in the tank. If the cross-sectional area of the tank (A) is knowable, the volume is the product of area and height, so

'deriv(A*h*?rho?,t)=F;i*?rho?;i-F;o*?rho?;o

An important simplification can be made if we assume that the fluid in the tank is well-mixed (perfect mixing). This lets us state that the properties of the fluid are the same at any point in the tank, including the exit, and consequently ?rho?;o=?rho?, T=To, etc.

Next, we want to look at the variables inside the time derivative -- which vary with time and which are constant? The tank is almost certainly rigid, so the cross section won't be changing. If we also assume that the fluids are incompressible, density will be constant with respect to time.

These simplifications let us write the total mass balance as:

A*?rho?*'deriv(h,t)=F;i*?rho?;i-F;o*?rho?
which we will often rearrange to the form
'quot(F;i*?rho?;i,A*?rho?)-'quot(F;o,A)
Integrating this equation will give us an expression for the level in the tank as a function of time. (This is called "simulation")

Another common simplification occurs if the density of the fluid entering the tank is approximately the same as that of the fluid leaving the tank (equal densities). You need to watch this one for reactors and mixing tanks, but it isn't a problem if all we're doing is heating water. Since the densities are the same, they cancel out, leaving

'deriv(h,t)='quot(F;i,A) - 'quot(F;o,A)

Energy Balance

If our tank has a coil or jacket so that the temperature is changed, an energy balance will be needed. The enthalpy of the flowing streams will be expressed in terms of the heat capacities and temperatures of the streams to produce

'deriv(V*?rho?*c;p*T,t)=F;i*?rho?;i*c;p;i*T;i-F;o*?rho?;o*c;p;o*T;o+Q
Assume perfect mixing and express the volume as the product of the cross-section and the height
'deriv(A*h*?rho?*c;p*T,t)=F;i*?rho?;i*c;p;i*T;i-F;o*?rho?*c;p*T+Q
Assume incompressible fluids and constant heat capacities:
A*?rho?*c;p*'deriv(h*T,t)=F;i*?rho?;i*c;p;i*T;i-F;o*?rho?*c;p*T+Q
and divide through to isolate the derivative
'deriv(h*T,t)='quot(F;i*?rho?;i*c;p;i*T;i,A*?rho?*c;p)-'quot(F;o*T,A)+'quot(Q,A*?rho?*c;p)
Once again, it is time to remember calculus. We'll apply the chain rule, to break the derivative into two parts
h*'deriv(T,t)+T*'deriv(h,t)='quot(F;i*?rho?;i*c;p;i*T;i,A*?rho?*c;p)-'quot(F;o*T,A)+'quot(Q,A*?rho?*c;p)
This sets things up to substitute the results of the mass balance in and eliminate the height derivative
h*'deriv(T,t)+T*'quot(F;i*?rho?;i,A*?rho?)-T*'quot(F;o,A)='quot(F;i*?rho?;i*c;p;i*T;i,A*?rho?*c;p)-'quot(F;o*T,A)+'quot(Q,A*?rho?*c;p)
Combine the T terms (some will go away) and
h*'deriv(T,t)+'quot(F;i*?rho?;i,A*?rho?)*T='quot(F;i*?rho?;i*c;p;i,A*?rho?*c;p)*T;i+'quot(Q,A*?rho?*c;p)
and this is the simplified energy balance for the system.

Depending on the purpose of our model, we might also need more explanation of the heat transfer term Q, for instance
Q=U*A*?Delta?*T
This is particularly important if the Q term depends directly on either h or T.


The Model

The differential equation model for the heated, stirred tank:

'quot(F;i*?rho?;i,A*?rho?)-'quot(F;o,A)
h*'deriv(T,t)+'quot(F;i*?rho?;i,A*?rho?)*T='quot(F;i*?rho?;i*c;p;i,A*?rho?*c;p)*T;i+'quot(Q,A*?rho?*c;p)

BREAK

Steady State Model

The steady state version of this equation is
F;ss*?rho?*c*T;ss=F;i;ss*?rho?*c*T;i;ss+Q;ss
T;ss=T;i;ss+'quot(Q;ss,F;ss*?rho?*c)
which is sometimes useful in evaluating modeling parameters.

Simulation

With the nonlinear model, it is possible to set up the system for a nonlinear simulation. The simulation will numerically integrate the model equation to produce transient response curve. The equation will need to be rearrranged slightly:

'deriv(T,t) = 'quot(F*T;i,V)-'quot(F*T,V)+'quot(Q,V*?rho?*c)
so that it is convenient to plug it into a numerical integration algorithm, such as the Euler method:
T(t+?Delta?t) = T(t) + 'deriv(T(t),t)*?Delta?t

Linear Model

The linear, deviation variable model of this system is:

'deriv(T,t) + 'quot(F;ss,V)*T = 'quot(F;ss,V)*T;i + 'quot(T;i;ss- T;ss,V)*F + 'quot(1,V*?rho?*c) + Q
It shouldn't be too difficult to obtain this directly, however, some may prefer to break the procedure into intermediate steps.

Linearization

The first step in obtaining this version is to use a Taylor series expansion (truncated after the first order terms) to linearize the model:

V*?rho?*c*'deriv(T,t) + [F;ss*?rho?*c*T;ss + F;ss*?rho?*c*(T-T;ss) + ?rho?*c*T;ss*(F-F;ss)] = [F;ss*?rho?*c*T;i;ss + F;ss*?rho?*c*(T;i-T;i;ss) + ?rho?*c*T;i;ss*(F- F;ss)] + Q
Then by using the steady state equation
V*?rho?*c*'deriv(T;ss,t) + F;ss*?rho?*c*T;ss = F;ss*?rho?*c*T;i;ss + Q;ss
and the definition of deviation variables
T=T;ss+T;dev
T;i=T;i;ss+T;i;dev F=F;ss+F;dev
Q=Q;ss+Q;dev
and quite a bit of substitution and cancellation, we get the linearized deviation variable model above.

Inputs & Outputs

At this point, it might be useful to note that this system has one output, T, and three inputs, F, Ti, and Q. Of the inputs, Q is a manipulation, and F and Ti are disturbances.

Laplace Transformation

Before doing the Laplace transformation, let use the residence time

?tau?;R='quot(V,F)
to "clean up" the differential equation
?tau?;R*'deriv(T,t) + T = T;i + 'quot(T;i;ss-T;ss,F;ss)*F + 'quot(1,F;ss*?rho?*c)*Q
then take the Laplace transform
(?tau?;R*s-T(0)+1)*T = T;i + 'quot(T;i;ss- T;ss,F;ss)*F + 'quot(1,F;ss*?rho?*c)*Q
where T, Ti, F, and Q are now functions of the Laplace variable s and no longer functions of time. Remember that the initial value of T is zero because the equation is in deviation variables. The equation can be rearranged to
T = 'quot(1,?tau?;R*s+1)*T;i + 'quot('quot(T;i;ss- T;ss,F;ss),?tau?;R*s+1)*F + 'quot('quot(1,F;ss*?rho?*c),?tau?;R*s+1)*Q

Transfer Functions

Because the system has one output and three inputs, three transfer functions are needed. The Laplace domain equation in transfer function form is:

T = G;D1*T;i + G;D2*F + G;P*Q
where
G;P = 'quot(1,?tau?;R*s+1)
G;D2 = 'quot('quot(T;i;ss-T;ss,F;ss),?tau?;R*s+1)
G;P = 'quot('quot(1,F;ss*?rho?*c),?tau?;R*s+1)
I tend to use GP for transfer functions relating outputs to manipulated inputs and GD for transfer functions relating outputs to disturbance inputs.

Solving the Linear Model

Because this is a linear model, we only need consider one input at a time. The solution if two inputs change will be the sum of the solutions changing one input at a time. Because we are in deviation variables, any inputs that do not change will be zero. In this example, we will consider an input change in Q only.

The linear model can be solved in two ways. The equation can be arranged into a standard form for a first order linear ODE

?tau?;R*'deriv(T,t)+T='quot(1,F;ss*?rho?*c)*Q
that can be solved using the integrating factor method, or (since this type of equation is so familiar by now) by inspection. For a step change in the heat input:
T(t)='quot(?Delta?Q,F;ss*?rho?*c)*(1- 'exp(e,_t/?tau?;R))

Solution by Integrating Factor

If "inspection" is intimidating, the integrating factor method may be used. The integrating factor for this equation will be:

p(x)='exp(e,'integ('quot(1,?tau?;R),t))='exp(e,'quot(t,?tau?;R))
The differential equation is multiplied by the integrating factor
'exp(e,'quot(t,?tau?;R))*'deriv(T,t) + 'quot(1,?tau?;R)*T*'exp(e,'quot(t,?tau?;R)) = 'quot(1,F;ss*?rho?*c)*'quot(1,?tau?;R)*Q*'exp(e,'quot(t,?tau?;R))
and the left hand side can then be contracted into a single derivative:
'deriv(T*'exp(e,'quot(t,?tau?;R),t) = 'quot(1,F;ss*?rho?*c)*'quot(1,?tau?;R)*Q*'exp(e,'quot(t,?tau?;R))
Integrating (given that Q is a step function) yields
T*'exp(e,'quot(t,?tau?;R) = 'quot(1,F;ss*?rho?*c)*'quot(1,?tau?;R)*?Delta?Q*?tau?;R*'exp(e,'quot(t,?tau?;R) + c;i
or
T = 'quot(1,F;ss*?rho?*c)*'quot(1,?tau?;R)*?Delta?Q*?tau?;R + c;i*'exp(e,'quot(_t,?tau?;R)
Since the problem is in deviation variables, the initial condition gives T=0, and so
c;i = _'quot(1,F;ss*?rho?*c)*?Delta?Q
and so
T = 'quot(1,F;ss*?rho?*c)*?Delta?Q*(1-'exp(e,'quot(_t,?tau?;R)))

Solution by Laplace Transforms

The equation can also be solved using Laplace transforms. Since Q is a step function

Q(s) = 'quot(?Delta?Q,s)
and so
T(t)='exp(L,_1)*[G;P(s)*Q(s)]
T(t)='quot(?Delta?Q,F;ss*?rho?*c)*'exp(L,_1)*('quot(1,s*(?tau?;R*s+1))
T(t)='quot(?Delta?Q,F;ss*?rho?*c)*'exp(L,_1)*('quot(1,s)+'quot(_?tau?;R,((?tau?;R)*s+1))
T(t)='quot(?Delta?Q,F;ss*?rho?*c)*'exp(L,_1)*['quot(1,s)-'quot(1,s+'quot(1,?tau?;R))]
T(t)='quot(?Delta?Q,F;ss*?rho?*c)*(1- 'exp(e,'quot(_t,?tau?;R))

Block Diagram of Process

Control Loop


R.M. Price
Original: 2/10/97
Modified: 5/8/2003

Copyright 1997, 2003 by R.M. Price -- All Rights Reserved