 The magazine of the Melbourne PC User Group Another Way To Skin The Cat Max Chandivert |
 |
Most of the software used on personal computers is ready-made commercial software, but a small program written for a specific purpose is often more effective in terms of results versus resources. Sometimes there is no ready-made solution and the ad-hoc solution is the only possible one.
Suppose for a moment that we are contemplating implementing an industrial process or designing a structure or mechanism involving any sort of dynamics. If we can represent the system's future behaviour on a computer through simulation, we can avoid functional pitfalls, including instability, we can test performance and optimise parameters without using any actual hardware.
The prerequisite for this is to be able to represent the characteristics of the system elements by means of mathematical equations. We select sufficient variables to describe the system state at any time and we write the equations which represent the relations between these variables. This is where the component manufacturer's data and the engineer's experience must combine to provide a realistic input. Without a reliable model to work with, garbage in yields only garbage out.
Unfortunately realistic input data often results in equations which are difficult or impossible to solve by the classical methods. The assumptions on which these methods rely are rarely verified in the field. They refer to theoretical representations of actual system components, and are made for the purpose of solving the resulting equations.
A case in point relates to accounting for energy dissipation in various processes. The classical method needs a form which maintains the linearity of the resulting differential equations ("DEs", from now on). If we want to account for dry friction with a threshold value and a further contribution proportional to a displacement or a power of the displacement and in phase with a velocity, we are in trouble. Add possibly the need to represent an arbitrary impressed force whose magnitude and timing may depend upon the value of any one of the other variables, and the classical methods become inadequate.
Classical solutions owe their unique popularity to the fact that we did not know how to treat more complicated equations. It is difficult to see the logics behind that argument. It is like loosing the door key when coming home in the dark and looking for it under the lamp post because that is the only place where there is light.
For most of those previously intractable or difficult cases, the advent of affordable and powerful computing power has opened possibilities by providing the means for practical numerical solutions. Instead of evaluating a formula resulting from a theoretical solution of linear DEs, values of the system variables and their derivatives are calculated as functions of an independent variable from a given set of initial conditions. These values can then be listed as output data, or plotted for visualisation. The forces included can account for any coupling between the variables, or any unorthodox non-linear form. If the DEs can be numerically integrated and we can input an initial state, we can "simulate" the consequent behaviour. Not only does it renders easy a task formerly difficult or impossible, but changing the values of parameters is instantaneous, and amending the format of force representations is just as simple, a great advantage over the classical solution.
The procedure on which this numerical solution rests is the numerical integration of DEs, a noted example of which is the refinement derived by S.Gill (ReŁ1) to the general method of Kutta, often known as the Runge-Kutta method.
I shall attempt to give a simple description of that method. The independent variable can, of course, b anything but in many cases, and for the sake of this explanation, it will be the time. We choose a small interval of time, say h. The procedure consists in evaluating a number, usually four, of evermore accurate estimations of the values of the variables at the end of the time interval from their values a the beginning. For the first interval we use the system starting conditions, then we proceed step by step, the values of the variables at the end of an interval becoming the starting values for the next one.
Given the time, Y0, and four dependent variables, Y1 to Y4, we must evaluate them over each interval. This is done by calculating increments, K0 to K4, to apply to each of the four estimations within each interval. Listing 1 below lists some Turbo Pascal code which is part of the procedure. The basic relation for evaluating the Ks is:
Ki = h * dYi/dt, i = 0 to 4.
Simplifications will result in most cases, and in particular, if YO is the time, dY0/dt = 1 and KO = h. The size of his easily determined. If h is halved without change in the results, then h size was satisfactory. Furthermore h size can be altered midstream, and can even be adapted to various phases of the process. Of course the smaller the value of h, the more calculations are executed and the slower is the simulation.
The Runge-Kutta method was derived originally for first order DEs, i.e. including no variable derivatives higher than the first. However it is readily extendable to systems of simultaneous first order DEs, and, since a higher order DE can be replaced by a corresponding number of first order DEs by introducing suitable auxiliary variables, the method is widely applicable.
Listing 1 shows code fragments of a Turbo Pascal program used to study the behaviour of a flat plate placed in an airstream for various positions of the centre of gravity and various values of the airspeed. The plate can move in the vertical direction and also rotate about a transverse axis. It is therefore a 2 degree-of-freedom system and the variables with the Runge-Kutta method are:
| Y0 : time, Y1 : up/down displacement of c.g., Y2 : velocity in Y1 motion, Y3 : rotation angle, Y4 : velocity in Y3 motion |
Obviously the derivatives of three of these variables, and therefore the Ks, have very simple forms. They are dY0/dt = 1, dYl/dt = Y2, and dY3/dt = Y4. The remaining two, dY2/dt and dY4/dt, are accelerations, and, when written per unit mass or unit moment of inertia, have a familiar form:
acceleration = -(dissipative force) -(elastic f.) + impressed f.
where each of the righthand side terms can include coupling with the other variables. In the short (6k) program from which
Listing 1 is extracted, the dissipative force in dY2/dt (up/down) includes a threshold value and a term proportional to the modulus of displacement, the elastic terms are coupled with rotation, and the impressed forces include the effect of up/down velocity and that of rotation. In dY4/dt (rotation), dissipation forces have the form of viscous damping, the elastic forces are cross-coupled and the impressed moment results from the lift force in dY2/dt. This is evidently a very crude representation with no practical engineering merit, except to demonstrate the procedure. The program plots two phase plane graphs side by side on the screen, and the user can interrupt the simulation at any time. In spite of the model simplified representation, the behaviour of the system is clearly illustrated. The airspeed is varied from 20 to 100 m/sec, and
"e", the distance of the c.g. behind the rotation axis is varied from 0.2 to 0.5 m. The system is stable through the whole speed range with e = 0.2, but unstable at the higher speed with e = 0.5 or more.
A further remark can be made about h. If the process studied is such that the terminal conditions are known instead of the initial ones, it is equally acceptable to proceed backwards using a negative h. This could be useful for instance to study a trajectory from the end point instead of the starting point, which may be unknown. It works just as well backwards!
I have written several short programs using either the original Runge-Kutta procedure or the Gill refinement on an V20 XT with 8087 co-processor in Turbo Pascal 5.5. Even on a slow machine the simulation speed is quite acceptable. I have coded a problem treated first in 1969 on a Digital Corporation Equipment PDP-10, and obtained completely identical results.
Reference
1: S.Gill, "A process for the step-by-step integration of differential equations in an automatic digital computing machine." Proc. Cambridge Phil. Soc., 47: 96--108 (1951)
Program
Listing - download in text format
|
|
|
|
Reprinted from the September 1993 issue of PC Update, the magazine of Melbourne PC User Group, Australia
