Temperature control of a water tank
Author's Data: Housam BINOUS
Department of Chemical Engineering
National Institute of Applied Sciences and Technology
Tunis, TUNISIA
Email: binoushousam@yahoo.com

This problem is solved using Berkeley-Madonna in Chemical Engineering Dynamics by J. Ingham et al. (second edition, Wiley-VCH, 2000). The solution using Mathematica is given in the present notebook. We would like to heat water with a heater. Oulet temperature must be equal to 80°C. Inlet water temperature is equal to 20°C. Water flow is equal to 10 m3/hr. Tank volume is equal to 100 m3. Water density and heat capacity are 1 Kg/m3 and 4.19 KJ/(Kg °C). The heater and the temperature measurement device have a first order delay with a time constant equal to respectively TauQ and TauSens. Initial heat input and tank temperature are equal to 2500 KJ/hr and 20°C.

Proportional control.  Kp is the proportional control constant.

V = 100; F = 10; rho = 1; cp = 4.19; Q0 = 2500; TRset = 80; T0 = 20; 
    TauSens = 2; Kp = 300; TauQ = 1;

QC[t_]=Max[0,Q0+Kp (TRset-Tsens[t])]

sol1=NDSolve[{V rho cp D[TR[t],{t,1}]== F rho cp (T0-TR[t])+Q[t],
      D[Q[t],{t,1}]== (QC[t]-Q[t])/TauQ, 
      D[Tsens[t],{t,1}]== (TR[t]-Tsens[t])/TauSens,TR[0]== 20,
      Tsens[0]== 20,Q[0]== 2500},{Q[t],Tsens[t],TR[t]},{t,0,150}]

There is a delay between the tank temperature and the measured temperature. This is due to the sensor first order lag (TauSens=2).

plt1=Plot[{Evaluate[Tsens[t]]/.sol1,Evaluate[TR[t]]/.sol1},{t,0,150},
    PlotStyle->{RGBColor[1,0,0],RGBColor[0,0,1]},PlotRange->All]

Tank temperature (TR=79.959) is slightly less than the set point = 80°C. It is commun for proportional control to have a small offset error. Thus, we need to use a PI or a PID controller.

Evaluate[TR[t]]/.sol1/.t->150

Heat input to water tank.

Plot[Evaluate[Q[t]]/.sol1,{t,0,150},PlotStyle->RGBColor[0,1,1],
  PlotRange->All]

Proportional-Integral control. Kp and TauI are the proportional and the integral control constants.

V = 100; F = 10; rho = 1; cp = 4.19; Q0 = 2500; TRset = 80; T0 = 20; 
    TauSens = 2; Kp = 300; TauI = 9; TauQ = 1;

QC[t_]=Max[0,Q0+Kp (TRset-Tsens[t])+Kp/TauI Int[t]]

sol2=NDSolve[{V rho cp D[TR[t],{t,1}]== F rho cp (T0-TR[t])+Q[t],
      D[Q[t],{t,1}]== (QC[t]-Q[t])/TauQ, 
      D[Tsens[t],{t,1}]== (TR[t]-Tsens[t])/TauSens,
      D[Int[t],{t,1}]== TRset-Tsens[t],TR[0]== 20,
      Tsens[0]== 20,Int[0]== 0,Q[0]== 2500},{Int[t],Q[t],
      Tsens[t],TR[t]},{t,0,150}]

Tank temperature presents oscillatory behavior. Thus, we need to use a PID controller.

Plot[{Evaluate[Tsens[t]]/.sol2,Evaluate[TR[t]]/.sol2},{t,0,150},
  PlotStyle->{RGBColor[1,0,0],RGBColor[0,0,1]},PlotRange->All]

Tank temperature is exactly equal to the set point = 80°C.

Evaluate[TR[t]]/.sol2/.t->150

Heat input to water tank.

Plot[Evaluate[Q[t]]/.sol2,{t,0,150},PlotStyle->RGBColor[0,1,1],
  PlotRange->All]

Proportional-Integral-Derivative control. Kp, TauI and TauD are the proportional, the integral and the differential control constants.

V = 100; F = 10; rho = 1; cp = 4.19; Q0 = 2500; TRset = 80; T0 = 20; 
    TauSens = 2; Kp = 300; TauI = 9; TauQ = 1; TauD = 2;

QC[t_]:=Max[0,
    Q0+Kp (TRset-Tsens[t])+Kp/TauI Int[t]+Kp TauD D[ (TRset-Tsens[t]),{t,1}]]

sol3=NDSolve[{V rho cp D[TR[t],{t,1}]== F rho cp (T0-TR[t])+Q[t],
      D[Q[t],{t,1}]== (QC[t]-Q[t])/TauQ, 
      D[Tsens[t],{t,1}]== (TR[t]-Tsens[t])/TauSens,
      D[Int[t],{t,1}]== TRset-Tsens[t],TR[0]== 20,
      Tsens[0]== 20,Int[0]== 0,Q[0]== 2500},{Int[t],Q[t],
      Tsens[t],TR[t]},{t,0,150}]

With a PID controller, temperature oscillations are quickly damped.

Plot[{Evaluate[Tsens[t]]/.sol3[[2]],Evaluate[TR[t]]/.sol3[[2]]},{t,0,150},
  PlotStyle->{RGBColor[1,0,0],RGBColor[0,0,1]},PlotRange->All]

Difference between set point and actual tank temperature.

Plot[Evaluate[TRset-TR[t]]/.sol3[[2]],{t,0,150},
  PlotStyle->RGBColor[0,1,0],PlotRange->All]

Heat input to water tank.

Plot[Evaluate[Q[t]]/.sol3[[2]],{t,0,150},PlotStyle->RGBColor[0,1,1],
  PlotRange->{-1000,20000}]

Tank temperature is exactly equal to the set point = 80°C.

Evaluate[TR[t]]/.sol3[[2]]/.t->150