Section 4.4.8: Applications of ODE's

Control

In the gravity flow tank we might control the outflow with a flow controller.

Without control:

\begin{displaymath}F_2 = k \sqrt{M} = c \sqrt{p_M} \end{displaymath}

where $ p_M = \frac{Mg}{A} $,
A is the area of the tank.

With control:

\begin{displaymath}F_2 = \sqrt{\frac{p_M}{\frac{1}{c^2} + \frac{1}{c_v^2}}} = \sqrt{\frac{M}{\frac{1}{k^2} + \frac{1}{k_v^2}}} \end{displaymath}

kv is a parameter which varies with valve position:

\begin{displaymath}0 \leq k_v \leq k_{v,max} \end{displaymath}

Flow Controller

Simplest type: proportional controller.

Equations:

error = set point - measured value
output = proportional gain * error
or

\begin{eqnarray*}e = 100 (F_s - F_2) / F_{range} \\ y = K e \end{eqnarray*}


Usually scale signals 0 - 100 % of range.

Assume valve 50% open at zero output signal, and linear response.
Might think

\begin{displaymath}k_v = \frac{(y+50)}{100} k_{v,max} \end{displaymath}

But valve cannot close beyond 0 or open beyond 100%:

\begin{displaymath}k_v = k_{v,1} \equiv \frac{(y+50)}{100} k_{v,max} \end{displaymath}

if 0 < kv,1 < kv,max
kv = 0 if $ k_{v,1} \leq 0 $
kv = kv,max if $ k_{v,1} \geq k_{v,max} $
Model contains discontinuities!

kv depends on the flow F2 through the error e and output y. If a single ODE is used:

\begin{displaymath}\frac{dM}{dt} = F_1 - F_2 (= f(M)) \end{displaymath}

we need to solve an equation iteratively for F2 given M at each function evaluation:

\begin{displaymath}F_2 = \sqrt{\frac{M}{\frac{1}{k^2} + \frac{1}{(k_v(F_2))^2}}} \end{displaymath}

It might be better to leave this as an algebraic equation and treat the system as a DAE (see section 4.4.9).
So F2, kv, e, y etc. would be kept explicitly in the model and an equation would be added to define each extra variable.

Level Controller

Instead of trying to control the flow, try to control the level.(For this example, assume level $ \propto $ holdup.)

Similar equations, but

e = 100 (Ms - M) / Mrange

where Ms is the `holdup setpoint', and Mrange is the range of holdup for 0 - 100 % measured value signal.

Also, y = K e,
$ k_{v,1} = \frac{(y+50)}{100} k_{v,max} $ and
kv = fk(kv,1), as before.

Direct evaluation of F2 given M in this case - no iteration.

Example

k = 0.2 kg0.5/s, kv,max = 1.0 kg0.5/s.
Assume Mrange = 500 kg, Ms = 300 kg, and controller gain K = -2.
Minus because low level $\rightarrow$ positive
error $\rightarrow$ negative output, which is needed to close the outflow valve.

Initial condition: M0 = 100 kg
Feed flow: F1 = 4 kg/s until t = 180 s, then F1 = 2 kg/s

Results

Results are tabulated at the end of this section.

Controller fails to reach set point (300 kg) but does achieve steady state. This is normal for `proportional only' controllers. If the current steady state is not exactly the one designed for, a continuing offset is needed give a non-zero error and maintain the control action.

In this case, zero output gives 50% valve opening: kv = 0.5 kg0.5/s
For 300 kg holdup, this needs an outlet flow of $ \sqrt{\frac{300}{\frac{1}{0.2^2} + \frac{1}{0.5^2}}} = \sqrt{\frac{300}{29}} = 3.216 $ kg/s
Since the feed flow is finally 2 kg/s, we need a smaller final flow, more closed valve (18.2%) and hence positive 'error': set point > measured holdup.

Thermodynamics

The independent variable need not be time!

ODE solving methods are useful in thermodynamics, where the independent variable may be a thermodynamic variable.

Example: Pressure dependence of enthalpy at constant temperature


dh = T ds + v dp

h - specific enthalpy (J kmol-1),
T - temperature (K),
s - specific entropy (J kmol-1 K-1),
v - specific volume (m3 kmol-1),
p - pressure (Nm-2)

Find change in enthalpy from a standard pressure p0 (e.g. 1 atm = 1.013$\times$105 Nm-2), to required pressure p1.

Treat T and p as independent variables:

\begin{displaymath}dh = T \left(\frac{\partial s}{\partial T}\right)_p dT + T \left(\frac{\partial s}{\partial p}\right)_T dp + v dp \end{displaymath}

Thus at constant T = T1,

\begin{displaymath}\left(\frac{\partial h}{\partial p}\right)_T = T_1 \left(\frac{\partial s}{\partial p}\right)_T + v \end{displaymath}

Maxwell relations give

\begin{displaymath}\left(\frac{\partial s}{\partial p}\right)_T = - \left(\frac{\partial v}{\partial T}\right)_p \end{displaymath}

So

\begin{displaymath}\left(\frac{\partial h}{\partial p}\right)_T = - T_1 \left(\frac{\partial v}{\partial T}\right)_p + v \end{displaymath}

Equation of state

We would like to get

\begin{displaymath}v = \Phi(T,p)\end{displaymath}

from the equation of state (EOS), with T = T1. If EOS is in the above form, then no problem - just evaluate whenever volume is needed. But usually EOS is pressure explicit:

\begin{displaymath}p = \Pi(T,v)\end{displaymath}

e.g. Redlich-Kwong EOS

\begin{displaymath}p = \frac{RT}{v-b} - \frac{a}{T^{\frac{1}{2}} v(v+b)} \end{displaymath}

(a,b constants for given chemical).

We could solve iteratively for v at each value of p in the integration. But this means many iterative calculations between p0 and p1.

A smarter way

Use implicit differentiation of the pressure-explicit EOS to get expressions for $(\partial v/\partial p)_T$ and $(\partial v/\partial T)_p$ as functions of specific volume at temperature T1:

\begin{displaymath}(\partial v/\partial p)_T = \pi_1(T_1,v) \end{displaymath}

\begin{displaymath}(\partial v/\partial T)_p = \pi_2(T_1,v) \end{displaymath}

Do just one iterative calculation to get v0 at pressure p0. (h0 is the enthalpy at T1 and p0.) Solve the ODE system

\begin{displaymath}\frac{dh}{dp} = - T_1 \pi_2( T_1,v) + v \end{displaymath}


\begin{displaymath}\frac{dv}{dp} = \pi_1(T_1,v)\end{displaymath}

for h and v as functions of p at constant T1.

Section 4.4.8.1 shows how to get functions $\pi_1$ and $\pi_2$ for Redlich-Kwong.

Summary

Control

Dynamic simulation can show effect of control strategy and tuning constants (gain K).Latter not covered in this section. Proportional-only control leads to offset:set point is not usually achieved.

Thermodynamics

Coupled ODEs in a thermodynamic variable (such as pressure) can reduce the computational effort in solving for thermodynamic functions (e.g. enthalpy).

Table of output: Holdup controller example

step time Holdup Outlet flow % control signal
0 0.00 100.0000 0.0000 0.0000
1 10.00 140.0000 0.0000 0.0000
2 20.00 178.6650 0.1954 1.4660
3 30.00 208.2578 1.5985 13.3031
4 40.00 228.7844 2.2156 21.5138
5 50.00 244.8290 2.5444 27.9316
6 60.00 258.2513 2.7553 33.3005
7 70.00 269.8923 2.9068 37.9569
8 80.00 280.2064 3.0238 42.0825
9 90.00 289.4720 3.1183 45.7888
10 100.00 297.8756 3.1972 49.1502
11 110.00 305.5505 3.2647 52.2202
12 120.00 312.5971 3.3235 55.0388
13 130.00 319.0932 3.3752 57.6373
14 140.00 325.1017 3.4213 60.0407
15 150.00 330.6740 3.4627 62.2696
16 160.00 335.8534 3.5001 64.3414
17 170.00 340.6767 3.5341 66.2707
18 180.00 345.1755 3.5651 68.0702
19 190.00 339.3776 3.5250 65.7510
20 200.00 324.6830 3.4181 59.8732
21 210.00 311.0609 3.3109 54.4243
22 220.00 298.5167 3.2030 49.4067
23 230.00 287.0576 3.0945 44.8230
24 240.00 276.6902 2.9854 40.6761
.. ...... ........ ...... .......
27 270.00 252.1062 2.6644 30.8425
30 300.00 236.6112 2.3888 24.6445
33 330.00 228.1177 2.1995 21.2471
36 360.00 223.9512 2.0938 19.5805
39 390.00 222.0400 2.0421 18.8160
42 420.00 221.1926 2.0185 18.4771
45 450.00 220.8229 2.0080 18.3292
48 480.00 220.6627 2.0035 18.2651
54 540.00 220.5636 2.0006 18.2255
60 600.00 220.5452 2.0001 18.2181
66 660.00 220.5418 2.0000 18.2167
72 720.00 220.5411 2.0000 18.2165


Next - Section 4.4.9: DAE Systems
Return to Section 4.4 Index
Return to Section 4 Index

Course Organiser
Last modified: Thu Aug 6 11:36:50 BST