Section 4.1Ordinary Differential Equations

In the dynamic modelling of chemical engineering systems, ordinary differential equations (o.d.e.s) occur almost exclusively as unsteady state material balances. We will therefore deal only with the kind of o.d.e.s that occur in this context.

Single O.D.E.s

The simplest single equation of the appropriate form is:

dy(t)/dt = -y

In this equation t, representing time, is called the independent variable. The single unknown or dependent variable is y. The solution is strictly speaking y(t), i.e. y as a continuous function of t. Where the equation is solved numerically this will be a set of discrete values of y at corresponding values of t.

To enable o.d.e.s to be solved it is necessary to have a known fixed point in the solution space. For the present type of problem this will invariably be a starting or initial value of y at an initial time normally designated as t=0, i.e.:

y = y0 at t=0

This equation is one of the few o.d.e.s that has an analytical solution. If we solve it with initial condition:

y = 1 at t=0

We obtain:

y(t) = exp(-t)

To obtain a numerical solution we write:
dy/dt ~ dy/dt

Here dy and dt represent small changes in y and t rather than differentials. Hence we can calculate the change in y for a given change in t or timestep dt approximately as in general:

dy = dy/dt dt

Hence we can calculate a `new' y value since:
y(t + dt ) = y (t) + dy = y (t) + dy(t)/dt dt

For this particular equation therefore:
dy = - y dt

And so:
y(t + dt ) = y(t) - dt y(t)

Starting with a known value of y = y0 = 1 at t = 0:
y1 = y(dt) = - y0 dt + y0
y2 = y(2dt) = - y1 dt + y1
y3 = y(3dt) = - y2 dt + y2
... etc.

This numerical solution is implemented here. Try changing the size of the timestep in cell B7 between say 0.01 and 0.5 and see how the numerical solution compares to the analytical one.

Timeconstants, Timesteps and Stability

The equation used as an example above is a special case, with k = 1:

dy(t)/dt = - k y

Time Constants and Time Steps

Note that k has dimensions of time-1 and so we can alternatively write:

dy(t)/dt = - y / T

T is called the time constant of the equation. The analytical solution of this equation is now:

y(t) = exp(-t / T )

Here is a solver for the modified problem. You can change the time constant. If you experiment with different time constants and time steps, you will observe that the accuracy of the solution depends on the ratio of the time step. Be careful at this stage to keep the time step always less than the timeconstant.

In discussion of lags in `black box' models it was shown that the time constant of a linear o.d.e. is a measure of how fast its solution progresses towards a final steady state. E.g. 63% in a time equal to the time constant and 90% in approximate 2.3 times the time constant. If we wanted to observe 99% of the response towards steady state this would require a time given by:
0.01 = exp (-t/T)

This gives gives a time approximately 4.6 times T.

Stability

Now make the time constant and time step exactly the same, e.g. 1. Note that the solution goes immediately to zero after the first timestep and remains there. The solution behaves in general as it should, i.e. it starts at 1.0 and reduces monotonically to zero, the correct final steady state vale, and stays there.

Now make the time step just larger than the time constant, e.g. 1.1. Observe that the solution now starts to behave in a quite different, and completely incorrect oscillatory manner. If the time step is made exactly twice the time constant the solution oscillates between a maximum and minimum value. If the timestep is made any larger, then the oscillations will grow and eventually become infinite.

Sets of O.D.E.s

It is possible to solve several equations which make up a set of o.d.e.s using the above method applied to each of them. For example to solve:

dx/dt = - x
dy/dt = x - y

Starting from initial values y0 and x0 we just apply repeatedly the two relationships:
x1 = - x0 dt + x0
y1 = ( x0 - y0) dt + y0

`Stiffness'

With sets of equations the choice of timestep becomes a significant issue.

It is clear from our investigation of stability that we must not use a time step greater than the time constant, and for reasonable accuracy it should normally be less than one tenth of the timestep. Since we would usually wish to follow most of the behaviour of the system towards its final steady state, since 4.6 times the time constant takes us 90% of the way there for a single equation a total simulation time of perhaps 5 times the time constant should be adequate,

For a single equation, or a pair of equations like those above which have the same time constant, a total of around 50 time steps would be adequate to provide sufficient accuracy and to cover the range of `interesting' behaviour.

Unfortunately, for a set of equations, the shortest time constant governs the acceptable choice of time step and longest determines the length of time for which the simulation must be run. Experiment with this using the program supplied. With the second time constant of say 5 the system will still be a long way from its final steady state, which is both variable and equal to zero, after a simulation time of 2.5. However, making the first time constant too small, i.e. less than 0.1 with the time step of 0.1, still causes instability.

This phenomenon is call stiffness and becomes a problem whenever time constants differ by more than an order of magnitude. Since in practice these may differ by several factors of 10, more sophisticated numerical methods than the one described above often have to be used to allow larger timesteps to be taken.


Next - Section 4.2: The Systematic Construction of Dynamic Models