Exact simulation of a special first order partial differential equation

196

OF

Annales de l'Association inlernationale pour le Calcul analogique

A SPECIAL

N o 4 - - Octobre 1963

EXACT SIMULATION ORDER PARTIAL DIFFERENTIAL

FIRST

b y A r v i d B. H E M P E L

EQUATION

**

SUMMARY The paper describes a method which makes it possible to simulate a special partial differential equation of the first order on an analog computer without any approximations. The method usually requires little computer equipment, although a distance velocity lag simulator is necessary in most cases. The simulation may be performed with either of the independent variables as continuous for any given value of the other. Examples are given of both ways of simulation. The partial differential equation of the form

(Oy/Ot) + v (Oy/Ox) + ay = F (x, t)

(1)

is often encountered when studying the dynamic behaviour of chemical processes (for instance heat exchangers, and tubular chemical reactors). The usual procedure for simttlating this equation on an analog computer, is to approximate the partial equation by a set of ordinary differential equations by writing one of the derivatives as finite differences. To obtain a satisfactory approximation, this method will require much equipment, and especially if this equation is simulated as part of a larger system, lack of equipment will often lead to poor approximations. It will be shown in the following that the special case where

F(x,t) = t(x)

through some sort of apparatus where the property, y, is influenced by the property f(t) acting over the entire length of the apparatus. An example is a heat exchanger, y will be the temperature, and v the velocity of the liquid, f(t) is equal to az, where z is the temperature of the tube wall. The constant, a, is determined by the heat transfer coefficient and the heat capacity of the liquid. By Laplace transforming equation 2 with respect to t, the following ordinary differential equation is obtained

v(dy/d×) +

(s+a) Y = F(s)

(3)

and with the initial condition Y =

0

for

x =

0

Y = Y (x, s) and F (s) are the Laplace transforms of y = y (x, t) and f(t) respectively, i.e.

or more generally if

~o

N

Y (x, s) = fo

F (x, t) = X g, (x) t, (t)

e-.t y (x, t) dt

l---I

the partial differential equation can be simulated without approximations, provided equipment for simulating a distance velocity lag is available. The simple case where F (x, t) = f(t) will be shown first, followed by a treatment of the general case.

Consider the equation (2)

with initial and boundary conditions y = 0

for

t =

y=O

for x = 0

s+a Y (x, s) =

Simulation with t as the continuous variable.

OY/30 + v (3y/~x) + 'V = /(t)

and similarly for F(S). A thorough treatment of the Laplace transform theory also of partial differential equations, may be found in (1). Solving equation 3 with respect to Y yields

0

[ ' ~ - - - a a (1 - - e

)1 F(s)

(4)

The initial condition in x has been used to specify the solution. For a given value of x = xo the term in brackets is the transfer function between Y (Xo, s) and F(s), and this may be realized on the computer with a distance velocity lag simulator and some linear equipment. The block diagram is shown in figure 1. The concept

v and a are constants. Physically this equation represents the behaviour of some property y in material moving at a velocity, v, * Manuscrit received January 20, 1963. ** Chr. Michelsen Institute, Bergen, Norway.

v

Fig. 1.

A.B. Hempel .' Special first order partial differential equation <> is defined in textbooks on control theory or linear analysis, for instance Brown and Campbell [2]. Methods for simulating systems with a given transfer function may be found in the book by Johnson [3]. This simple structure will also appear in the more general case, as will be shown in the fallowing. General case. The equation is now 1'4

fay/at) + v(Oy/ax) + ay =

g,(x) /,(t)

X

(5)

i=l

with initial and boundary conditions y --- yo(t)

for

t = 0

y = y~(x) for x = 0 Following the same procedure as before, the Laplace transformed equation becomes N

~, (dY/dx) + (s + a) Y = yo(X) + X g, (x) F, (~) (6) I=X

s+a

- - - - X

e

197

N

X F, (s) H, (s, o)

V

i=X

A simulation of the transfer ftmction H (s, x0) and I (s, Xo) as given by equation 9, may seem rather complicated at first. It must be remembered that the expressions in equation 9 are the general expressions, and in practical applications the sums will contain only a few terms, or the integrals appearing in equation 7 may be integrated directly without using partial integration. This means that the sums in equation 9 may be expressed in closed forms. What is achieved by the general expression given in equation 8, is to show that .r+a the exponential term e v always will appear as a factor in the solution of the integral. This is the main point in this method of simulation, because here it appears that the only trancendent part of the transfer function is a distance velocity lag. With equipment to simulate this delay, the expression given in equation 10, may be simulated according to the block diagram in figure 2.

with initial condition Y = Yl(s)

for

x = 0 F'"---L~

The solution of this specified by the initial condition may be written: s+ a

r~s_.L

. * -" ..........

4--

l (h t~)

Y (s, x) = Y, (s) e s+a

- - - - X

+ e

s+a

£x

------X

v

yo(x) e

--,

s+a

*'

+(:/vie

,X

v

dx

N

s+a

x

,:,EF'(O~

--

÷

X

g,(xle *'

,~< (7)

By continued partial integration it is found that 1 /

s+a gl(x) e - ' 7 - - x

Y,Is)

dx Fig. 2.

-s-+ a x oo (_ v)j-~ = e v X ~ g,")(x)

j=~ (s + a)J

(S)

Introducing the symbols oo

H, (s, x) =

X (__~)a-__~ g,"' (~) ~:, (s + a)J

and

(9)

we may write the solution --

Y ( s , x ) ----- Y,(s) e

s+a

s+a

- - I(s,o) e

X

v

~, x] +

+

[I(s,x)

N

X F,(~) I-I~(~,x) l=X

0o)

It is seen that the first term in the solution given in equation 10 originates from the boundary condition a t x -- o. It consists of a pure distance-velocity lag multiplied by the constant e- (a/v) Xo. The second term (the terms in brackets) reflects the influence of the initial condition on the solution. The excitation f u n o tion to this term must be an impuls function, whose Laplace transform equals 1, or the term may be multiplied with s, and given a step function whose Laplace transform is l/s, as input function. The term is usually of little interest, because the equation may be studied around a stationary value, and in this case the initial condition will be zero. In any case, however, the effect of the term in brackets will vanish for t >/ xo/v, which is easily seen by using the superposidon principle and noting that

I,(Xo) = e -~l/°~x°''' I,(0)

(11)

Annales de l'Association internationale pour le Calcul analogique

198

The last two terms gives the solution due to the N

input function ~

gi (x)fi (t), which is the result must

i=l

often sought. As an example is chosen the equation

~y/~t + ~,y/~x + ay = x[(t)

'

y -- 0

for

y--0

for x - -

(12)

N O4 ~

continuous function of x at specified values of time. The procedure is exactly as before, except that the equation is Laplace transformed with respect to x, and the time functions f, (t) now must be specified before the problem is programmed for the computer. There is no reason for repeating the procedure for the general case, instead the following example will be shown. The equation

t -- 0

(ay/~t) + v(~y/~x) + a y -

0

Octobre 1963

sin cot g(x)

(14)

is chosen, with boundary conditions

The Laplace transformed solution of this equation may be written ~. s+a 1 v - ~ x Y(s,x) -[x (l--e v )] (13) s+a s+a Putting x -- xL the simulation is made according to the block diagram shown in figure 3.

Y(% Xo)~

y =

0

for x -

y =

O, for

0

t--

0

Laplace transformation with respect to x gives

(dY'/dt) + (vs + a) Y' = Y' -- 0

for

sin~ot O(s) t -- 0

The function Y' = Y' (t, s) is the Laplace transform of Y (4 x) and G(s) is the Laplace-transformed: of g(x). The solution of equation 15 specified by the initial condition is

(vs + a)sin o~t y ' ( t , . 0 -- [

Fig. 3. The solution giving y (t, xL) for f(t) equal to, a step function is shown in figure 4. The responses at the points A, B and C in the block diagram are also given in figure 4. The values in the simulation are v -- 0,25, a = 0,5 and xr, = 1.

(vs + a)" + ~-

lOv

r l'i-'.... I[

---B 0 lOv

[ i

I

I'

5

,/"

,

I t I 1 1 I I I I

(cos,~t- ,-~ *-'~)l GO) (16)

To simplify the simulation this may be written

Y' (t, s) -

I I I

( v s , + o) ~ + o,=

....

10 sek

vs + a

(vs + a) (vs + a) ~- + ~

[sin ~ot

(co.s.o,t- ,-=~ ,-,'~)] G(.O

(17)

• I,

Ft,.xo~

0

At specified values o.f t -- 4 this may be simulated according to the block-diagram in figure 5. 10 sek

5

GCs)

_ ~

÷ j ~ . [

,,.o.a

I

Y'(ti,s)

5v

--e

...... i q

L

5

--- . 'lI, II,

0

5v

10 sek

' ~

--e

i

-1

I0 sek

IBm

C

5

10.~

Fig. 4.

Simulation with x as continuous variable. The method as described so far gives the continuous time solution at specified values of x. In process investigations for control purposes this is usually what is needed. It is also possible, however, to record y as a

Fig. 5. In this simulation the machine time, r, represents the x-value, so that the solution for r > xL (xL being the length of the equipment) has no significance to a physical problem. This brings out an interesting feature of this simulation. It is observed that the distance velocity time delay is vti. If the solution at a time ti is wanted, where ti > XL/V, then the distance velocity delay is larger than the solution time. This means that if only stationary solutions are of interest, which is usual when periodic solutions are studied, the block containing the distance velocity lag can be omitted.

A.B. Hempel "Special first order partial differential equation In figure 6 is shown the result of this simulation with g(x) -- a, i.e. G(s) -- a/s. The coefficients chosen are actual values from a heat exchanger. They are a - - 0,15, v - - 0,17, xL - - 1. The solution is shown for different values of ~o • o -- 0,53, 1,0'7, 1,60 and 2,14.

199

I

2..O

T0

,

m=, ! raw'flick.

'tO

0

g---

OJ

02

o,3

~

o,s

U

0,7

U

o.7

3

Fig. 6 c.

\

o,1

0

~

o,3

o,~,

a~

X~

U

o,7

U

U

N

Fig. 6 a.

o.1

n,2

For each of these frequencies the following values of

oa

~

o.5

U

ot

1,a

Fig. 6d.

t were chosen" t

-

-

ti --

2 rr. 10 ,-1

i

( i - - 1, 2, 3 , . . . . , 1 0 )

in the frequency response diagrams of a heat exchanger. This is sho.wn both experimentally and theoretically

by Hempel [4]. The curves shown in figure 6. may be looked upon as temperature profiles of the fluid along a heat exchanger having a sinusoidal temperature variation in the wall.

Conclusion.

Each curve is the temperature profile at a specific

It has been shown that in some special cases a 1.

time, t i . It is seen that damped standing temperature waves are obtained in the heat exchanger. At low frequencies fo,r the driving temperature variation in the wall, only

order linear partial differential equation may be simulated exactly on an analog computer with the aid of a distance velocity lag simulator. The special cases are those where the ordinary differential equatio,n which appears when the partial equation is Laplace-transformed with respect to one of the two, independent variables, can be integrated analytically. The solution may be recorded continuously in one of the variables and for any value of the other. If only stationary periodic soiutions at speci4ied time and with the space coordinate continuous are o,f interest, this may be obtained without the use of the distance velocity lag simulator. When this method can be applied, it gives the exact solution, and as rule little equipment is needed.

one node is present (at the inlet) (fig. 6a). At o~ -- 2 rr (v/xL), i.e. 0,0'7 rad/sec one node also occurs at the outlet (fig. 6 b).

This node is not sharp, i.e.

y

z.0

.......................

__.~

5 /

~

l

.

-

'

_

J

RE,FERENCE,S

[i] H.S. CARSLAW and J.C. JAEGER : <(Operational -:0

~ F - - .......... o.z

U

,

o.i

i

~

U x .--l,..

03

oJ

o,s

Lo

Fig. 6 b.

there is a small variation in the temperature. The sharpness is a fu.nction of the quantity ax~/v, being exact when this equals zero. At still higher frequencies one or more nodes will appear in the heat exch'ar~ger (fig. 6 c and 6 d). . This standing wave illustrates from another point of vmw the now well known resonance e~ffect appearing

methods in applied mathematics>>. Oxford University Press, second edition, 1947. [2] G.S. BROWN and D.P. CAMPBELL: <>. John Wiley & Sons, Inc., New York. Chapman & Hall, Ltd., London, 1948. [3] C.L. JOHNSON : <>. McGraw-Hill Book Co., Inc., New York, Toronto, London, 1956. [4] A. HEMPEL : <( On the dynamics of steam liquid heat exchangers>>. Trans. ASME. VoI. 83 (1961), Serie D, n r 2, p. 244.