__Dynamics of a Heated Tank with PI Temperature Control__

`> `
**restart;**

The energy balance for the stirred tank is

`> `
**EB := Diff(T[t],t)=(W*C[P]*(T[i]-T[t])+q)/(rho*V*C[P]): EB;**

where
is the temperature in the tank. The temperature of the exit stream is given by

`> `
**Teqn:=Diff(T[0],t)=(T[t]-T[0]-tau[d]/2*Diff(T[t],t))*2/tau[d]: Teqn;**

The thermocouple shielding and electronics are modelled by

`> `
**Teqn3:=Diff(T[m],t)=(T[0]-T[m])/tau[m]: Teqn3;**

The energy input to the tank is modelled by

`> `
**qeqn:=q=q[s]+K[c]*(T[r]-T[m])+K[c]/tau[I]*epsilon: qeqn;**

where

`> `
**eeqn:=Diff(epsilon,t)=T[r]-T[m]: eeqn;**

and

`> `
**qseqn:= q[s]=W*C[P]*(T[r]-T[i,s]): qseqn;**

The parameters in the model are

`> `
**params :={V=4000/rho/C[P],W=500/C[P],T[i,s]=60,T[r]=80,tau[d]=1,tau[m]=5,tau[I]=2,K[c]=0};**

The initial values of the differential variables are given as a list of equations.

`> `
**Start := [t=0,T[t]=80,T[0]=80,T[m]=80,epsilon=0]: Start;**

Several different approaches for solving this system of equations can now be followed. A difficulty is that the step change in the inlet temperature occurs not at
, but at
. One possible approach is to integrate the equations for the first 10 minutes, then make the step change and integrate the equations with the new parameters for as long as needed. Maple's own dsolve/numeric will accept a procedure for evaluation of the right hand sides of the differential equations. We can program the step change in that procedure. Here is our code for this particular system.

`> `
**deproc := proc(n,t,y,dy)**

local q,qs,k,epsilon;

global params,T;

T[tank]:=y[1];

T[0]:=y[2];

T[m]:=y[3];

epsilon:=y[4];

if t < 10 then T[i]:=60 else T[i]:=40 fi;

qs:=W*C[P]*(T[r]-T[i,s]);

q:=qs+K[c]*(T[r]-T[m])+K[c]/tau[I]*epsilon;

dy:=vector(n):

dy[1]:=(W*C[P]*(T[i]-T[tank])+q)/(rho*V*C[P]):

dy[2]:=(T[tank]-T[0]-tau[d]/2*dy[1])*2/tau[d]:

dy[3]:=(T[0]-T[m])/tau[m]:

dy[4]:=T[r]-T[m]:

for k from 1 to n do dy[k] := subs(params,dy[k]); od;

RETURN(eval(dy));

end:

The arguments (required by Maple) are the number of equations, the independent variable, the vector of dependent variables and a vector of right hand sides (computed by the procedure). Note where we have programmed the step change in the inlet temperature and that set of model parameters is declared to be a global variable. This facilitates changing the parameters without changing the procdure, as we will have to do for later calculations.

We solve the problem using dsolve/numeric using the following command. Note that by using the optional arguments we can tell dsolve/numeric to use the proceure above for the right hand sides. In this case the first three arguments serve no role except that dsolve/numeric expects them to be there. The last argument tells dsolve/numeric to return an array with the results tabulated as shown

`> `
**result:=dsolve({diff(x(t),t)=0,diff(y(t),t)=0,diff(w(t),t)=0,diff(z(t),t)=0}, {w(t),x(t),y(t),z(t)},type=numeric, method=rkf45, initial=vector([80,80,80,0]),start=0,procedure=deproc,value=array([0,seq(i,i=9..20)]));**

To extract the table, perhaps for plotting purposes we use the op command in the following form

`> `
**Ttable := result[2,1]:**

print(Ttable);

The first column is the time, the next three are the tank, thermocouple, and measured temperatures.

**A. Open Loop Performance**

Open loop performance is simulated by setting
, which is what was done for the above calculations. We need to integrate for a longer time, however, so we repeat the above command but integrate to 60 minutes.

`> `
**result:=dsolve({diff(x(t),t)=0,diff(y(t),t)=0,diff(w(t),t)=0,diff(z(t),t)=0}, {w(t),x(t),y(t),z(t)},type=numeric, method=rkf45, initial=vector([80,80,80,0]),start=0,procedure=deproc,value=array([0,seq(i,i=9..60)]));**

We extract the results table

`> `
**Ttable :=result[2,1]: **

and plot the three temperatures.

`> `
**with(linalg):**

Warning, new definition for norm

Warning, new definition for trace

`> `
**plot({seq([seq([Ttable[i,1],Ttable[i,k]],i=1..rowdim(Ttable))],k=2..4)},color=[red,blue,black],labels=[t,T]);**

**B. Closed loop performance**

This requires a change in
to 50 (it was 0 in the first case).

`> `
**params :={V=4000/rho/C[P],W=500/C[P],T[i,s]=60,T[r]=80,tau[d]=1,tau[m]=5,tau[I]=2,K[c]=50};**

repeat the integration

`> `
**result:=dsolve({diff(x(t),t)=0,diff(y(t),t)=0,diff(w(t),t)=0,diff(z(t),t)=0}, {w(t),x(t),y(t),z(t)},type=numeric, method=rkf45, initial=vector([80,80,80,0]),start=0,procedure=deproc,**

value=array([0,seq(i,i=9..50),seq(i*10,i=6..20)]));

and plot the results.

`> `
**Ttable :=result[2,1]: **

`> `
**plot({seq([seq([Ttable[i,1],Ttable[i,k]],i=1..rowdim(Ttable))],k=2..4)},color=[red,blue,black],labels=[t,T]);**

**C. Closed loop performance for **

We increase
to 500

`> `
**params :={V=4000/rho/C[P],W=500/C[P],T[i,s]=60,T[r]=80,tau[d]=1,tau[m]=5,tau[I]=2,K[c]=500};**

and repeat the integration, this time asking for results to be returned for every minute after the 9th so that we can make a nice looking plot.

`> `
**result:=dsolve({diff(x(t),t)=0,diff(y(t),t)=0,diff(w(t),t)=0,diff(z(t),t)=0}, {w(t),x(t),y(t),z(t)},type=numeric, method=rkf45, initial=vector([80,80,80,0]),start=0,procedure=deproc,**

value=array([0,seq(i,i=9..200)])):

The output table is extracted from the result:

`> `
**Ttable :=result[2,1]: **

and the three temperatures plotted as a function of time.

`> `
**plot({seq([seq([Ttable[i,1],Ttable[i,k]],i=1..rowdim(Ttable))],k=2..4)},color=[red,blue,black],labels=[t,T]);**

**D. Proportional Control**

Proportional only control is simulated by setting the term
= 0. We can accomplish this by giving
a very high value as shown here.

`> `
**params :={V=4000/rho/C[P],W=500/C[P],T[i,s]=60,T[r]=80,tau[d]=1,tau[m]=5,tau[I]=1e99,K[c]=500};**

We carry out the integrations for this situation.

`> `
**result:=dsolve({diff(x(t),t)=0,diff(y(t),t)=0,diff(w(t),t)=0,diff(z(t),t)=0}, {w(t),x(t),y(t),z(t)},type=numeric, method=rkf45, initial=vector([80,80,80,0]),start=0,procedure=deproc,**

value=array([0,seq(i,i=9..60)])):

The output table is extracted from the result:

`> `
**Ttable :=result[2,1]: **

and the three temperatures plotted as a function of time.

`> `
**plot({seq([seq([Ttable[i,1],Ttable[i,k]],i=1..rowdim(Ttable))],k=2..4)},color=[red,blue,black],labels=[t,T]);**

Let us repeat the last example with BESIRK rather than
dsolve
.

`> `
**read `e:/maple/numerics/integ/besirk`:**

`> `
**read `e:/maple/utils/utils.mpl`:**

BESIRK needs three arguments. The first is a list of the differential (and algebraic equations):

`> `
**T:='T':**

`> `
**DEqns := [EB,Teqn,Teqn3,eeqn]: DEqns;**

`> `
**Subs(DEqns,qeqn,qseqn,DEqns,right);**

`> `
**DEqns2:=evalf(subs(params,T[i]=40,%)): DEqns2;**

`> `
**read `e:/maple/thermo/td/tdtools.mpl`:**

`> `
**result := BESIRK(DEqns2,Start,0..60,hmax=1);**

We can plot the results using the datatableplot command from the BESIRK package.

`> `
**datatableplot(result,[[t,T[t]],[t,T[m]],[t,T[0]]],color=[red,blue,green]);**

`> `