Implicit simulation method for distributed parameter

W. RuDpel and M Ze:tz" D~sh:huted parameter

Ii~PLI[C[T

SIMULATION

17L

NETHOD W

-

DISTRIBUTED

PARANETER

RUPPEL and M ZEITZ (*'1

SUMMARY Thla paper p~esents a hyb, td method of s~mulatmg a class of non-l, nea~ d,s ~vbuted parameter systema. The implicit method based on an integral representation o/ the do'nam~c bound, oy value problem Th s integral equatwn ~s formulated using Green's function, defined by the hnear spatial &lfe,ent:al and boundwy operators The h),b~:d ~mplementation of the :mpkc:t ~:mulatwn ~eq:qres the analog mteg~atmn of ordinw 7 d~ffe~entad equations and, sHnultaneously, algeb,'mc operations on the ditigal pa~t. Due to the mult~ple.x-techmque used, the method zs wel[ suited to the s:mula wn of patlkt/ differential equations w'th a separable no~-hnear source te~m. As a t),pzca! e:,ample a catalytic fixed-bed ~eacto~ :s a/ululated on a 17.6red comp,/te~ by application o / t h e -

FOR

-

:/nphc:t method,

RESUME Cet a~tzcle ptdsente une mgthode h),b~lde pour la s:mala'zon d'une dasse de s),st&nes ,} param~tres contolus non hn~a.,res. Cette mdthode lmpl.clte est basde sm une relrdsentatlon inldgrale d'un pJobl&ne d),namtqae de valems aax limltes. Cette dquatmn /ntdgra/e est formulde au moyen d'une /onc. twn de Green, d@me ,) patter des op&atears :li//drentiel spa'~al et au.x limites. Le montage h),bJlde pour cette sHnulatzon imphcite ndcess te ]'~ntdgration analogJque d'@uatwna dff/drentelles oJdlnak'es el, aimultandment, des opdratmns alggbriques sur la pattie numdrtque. Grdce ,'t la techmque 1,.,,[email protected] utdes&, retie md:hode est parttcuL*Orement adaptde ,) la simulation d'@uat/ons aux d&vvdes partielles ave¢ an terme no J, lmdmre s@wable. Comme exemple d'apphcatmn, un ~&ctem catalytlque ~) be ft:,e est simuld sin" un calculateur h),brule par applications de la mdthode imph. -

-

cite.

1.. I N T R O D U C T I O N The simulation of &stributed parameter systems with extremely steep profiles using classical analog/hybrid or hybr, d methods often poses problems due to the computation time and/or the analog/hybrid equipment required. Catalytic fixed-Eed re'actors with highly exothermic reaction represent typical examples of such systems The mathematical descriptmn of these dynamic boundary value problems leads to non-linear partial differential equations which, in general, can only be solved by computer smmlafion. For th,s reason, the hybrid computer offers advantages concerning interactive simulation and relatively cheap computatmn time. 2. S Y S T E M

v2(1)

L.__.

I×(z,t) 1 Fig l -- One-dm~tnsmnal system (m,. -- 2)

with the time coordinate t and the normalized space coord,nate, z or C. The linear part of the pde. is represented by the dtfferent,fl operators

DESCRIPTION 111

Consider the class of one-dimensional distributed parameter systems in F,g. t, governed by p,utml differential equations (pdes) of the type D,~,(z,t) + D, x(z,t) = S ( u , , )

(l)

D, = 4' /,,,, ( ~ ) m . ~1

~lll

----

,,,.

C {t, 2, 4)

ana

(2,,)

~. z'" (%IL

D, = 2 ,t,, (z, t) . . . .

,

n,, E ( t , 2}

(2b)

z £ (o, 1), t > t , , ': Rewsed and extended v,.q~m, of the paper [l] piesented on 7 AICA-Congms, Prague 197 a, ** Dlpl-lng W Ruppel .tnd Dt-Ing Iv[ Zeitz, ln,~,mt fu, Systrmdynamik und Regelungstedm~k, Professm Dr-ling E D Gdles, University of Smttg,trt, [')-7000 Stuttgart t, Keplersnasse 17, Germany

which enable the description of plug flow, diffusion. heat conduction, electrical t r a n s m i s s i o n o f e l a s t , c syst e m s ]2] Other than the control variable u (z, t), all non-hnear terms of tt~e system variable x (z, t), are summarized m die source function S (u. x) of (1). At the boundarms, x (z. t) satisfies the m, conditions

/tnnale.r de l'Aaaocmt/on n;ternatwna/e poe;i" ]e Ca&u/ amdog:que - - N" 3 - - ]/:/l/et 1975

172

B .... a (z,/) = R,,, (v,,,, x) /;/ = t.2, .,.Ut,,, z E {0, 1}

(3)

defined by the boundary functions R,.(l,,,,, x) and by the linear differential operators B.... having a maximum order of t o o - - i . In addition to the input variables t',,, (t), the R,,, (v,,,, .~) can obtain time derivatives and/or non-hnear expressions of the system variable. "Ibis distribution into hnear and non-hnear parts is possible for rno~t pdes and boundary con&tions belonging to physical or industrial processes 121. The initial state of the system -----

, (z,t)

=

(4)

X,,(z)

,'l ~ 0, [ ..... ,~I~- - [

T h e solution :, (z, l) of the dynamic boundary value problem (1)-(4) is to be found, with dependence on the control variables, //(z, t) and t',. (t), and on the initial profiles X,, (z) by hybrid simulation. PRINCIPLE

For hybrid simulat,on, the pde (t) has to be transformed under consideration of boundary and initial conditions into ordinary differenti;d equations (odes) and/or algebraic equations For this purpose, most simulation methods use :t &scretizatmn of time and/or space coordinates or a transformation into the characteristic or modal coordinates [3]. In contrast to this, the lmphc# method is based on an integral transformation, with the Green's function g (z, LS) of the standardized boundary value problem as kernel. For standardization of pdes (1) with the s~mae order m,~, a proposal in [4] is i m p l e m e n t e d . The spatial differential operators D, of the general form (2a) are substituted by ~, t2t o

D,I~ --

q- `8"',,

GD, x ( z , l ) = - - , ( z , t )

+ GQ(u,x)

+ GQ(u,.x)

a,,,,, ( s ( . , . , ) -

E llk--{}

b,,,-aT) a,,,,, =~ o

+ ,8,,,,.

Ob)

Due to the standardization, the Green's function g (z, ~) as well as operations defined by D,t] and B .... do not depend on the coefficients h,. (z) By means of Green's identity, the pde (t) and boundary conditions (3) can be formally written as file integral equatmn

x (z, l)

= ,£ g(z,O [Q(u, x)--D~., 1 a'~ +/(R,,,) (6)

t- f(R,,)

(8)

and a further deterinln'.ttion of the variable x (z, I) by inversion of the operators, G and D r , which are on the left side of (7) y(z,t) = G - ' y ( z , / ) ., (z, t) =

Dr-'

(9)

{y (z, ,'), X,, (z))

(10)

The time integration (I0) requires the consideration of the initial conditions (4). 4. I M P L I C I T S I M U L A T I O N E Q U A T I O N S The simulation equations (7)-00) are considered, from the point of view of a hybr, d implementation, as being on parallel l,nes &, = t,, + k A t ,

k = 0,

i....

(i~)

m the t-z-plane (Fig 2) Under this assumption, the transformation G is equivalent to the calculation of w (z, t,,) =

G Q (u,x)

+ I (Rm)

D,~ u, (z, t,.) = Q (., ..), 13 .... w (z, t,,) =

(12a)

problem

solving the spatial boundary value

z c (o, t) R,,, (,',,,, .9

m = t , 2 ..... m,,,

q("' ' ) -

(7)

The Hnp//cit simulation principle consists of computanon of the right side of Eq. (7)

(5a)

This standardized operator contains the real parameter ,8 which can be chosen arbitrarily Both the term ,8'% x and the rem~zining derivatives of x are included in the new source function I m,-t ~m ,x

+ f(R,,,)

already yields the basic equation of the slmulatmn method. The operator G stands as an abbreviatmn for the integral m (6) The method is called imphH1 simulation since, dealing with the unknown variable .x (z, t) ,nvolves dealing with an implicit equation.

y (z,t) = --x(z,t)

t = t.

is given by n,, profiles X,, (z).

3. I M P L I C I T S I M U L A T I O N

if the dynamic term D t .x is considered as a source term The additional term f (R,,,) is caused by the non-homogeneous boundary conditions (3). Solving Eq. (6) for the transformed dynamic term

(t2b)

02c)

z E {0, t}

,f 0 0

~k÷1 ~k

unn~

~k-1 |

0

~,..

z

Fig 2 --- Dnpllck sm~ulation in the t-z-plane.

1

IF Ruppel and M

Zeltz ' Dtstjthuted paramete~

173

This can be realized by analog integration of two initial value problems if the integral,on results B .... w (1, tl,_t) are used to convert the given boundary values into reqmred initial values Next, the inversion (9) and the time integration (10) are approximated by algebraic operations for a digital implementatmn. Based on the definition (6), the exact equations of the reverse operator G -~ read D,I ~ y' (z, t,,) = ), (z, t,,)

z < (0,i)

B,,, y~ (~, t,,) = 0 n/= t,2 .... m.,

]m

1

/.-

1

1=

,

+

(L~b)

1,2,

,j,,

(14)

H )j+ (l,,)

( t 5)

D k ' {H y+ (t,,), X,,}

(t6)

into Eq (9) gt~es x (1,,+,) =

I

z < {0, I}

Llsing a statable difference approximatmn on Eq (13), a (/.)<. /0)-matrix H is obtained for the relation between the (/,, × t)-vectors, y (it,) and )'÷ (tt,), the components of which are the values, ), (za, tL,) and 1'+ (z a , tk), respectively. Substitution of )L(t,,) =

r

÷

(1%0

By means of these relationships, the original function y (z, tl,) is to be determined from the transformed function y~ (z, tk). It IS understood that the derivatives with respect to the space coordinate can not be realized m an exact way In addition, the following time integration can only be carried out for constant space coordinates. As a consequence of these facts, we have to limit ourselves during the inversion of G to the computatmn of y (z,, It,) for /. equidistant space coordinates z, -

initiated by the I O-pulses which correspond to the spatial discretlzatmn (14) In the pauses between I/O-operations the vectors, ~ (it,) and :x (IE,,,), 'are computed digitally for each sm~ulatmn cycle.

if Dt is approxm~ated by an explicit difference operator Dk As a result of this, simple algebraic operatlons yield the solution .y_(t~,+~) at the next time point. Finally, a hneat lnterpoiation circuit produces continuous profiles .x (z, t~,+~) which are used to compute the functions Q (,, x) and R,,, (v,,,, x) 5. H Y B R I D I M P L E M E N T A T I O N The unpliclt simulatmn is m~plemented according to the hybrid computation scheme in Fig. 3. Realization of the non-hnear expressions, Q (u, x) and R,,, (r,,, ,:,'), the computation (12) of the function w(z, th) and the linear interpolation occur on the analog part. The initial profiles X,, (z) and the transformed function y+ (a, tl,) are fed into the dlgltal part where the new system vector .x (tL,+t) i~, computed by the algebraic operations, H and D~, ~. The analog, hybrid and dlglt.d simulations steps are performed as shown in Fig 4. The analog integrations connected w~th the transformation and interpolation are carried out m high-speed repet~tive mode (IC-OP) The A/D-inputs and the D/A-outputs are

t_

.I

Fig :, --Hybrid ~theme of unpln'lt simulatmn Fig. 5 contains a sketch of the mformatlon flow of a hybrid unpkczt simulation. A Ioglc circuit in the analog computer controls the ent, re smlulatton : The digltal part can be considered to be an ad&tional component of the analog compute~ The digital program, which must normally be written m machine code, because of the high computation speed required, is built into a hybrid frame-program in interpreter language, to enable convement operatlon The interpreter program facdttates the preparation of the simulation in &:dog mode over the wdeo. The simulation Is operated by the keyboard of tb.e analog computer and the simulation results are displayed on an analog oscilloscope These results, aided by a sufficiently high computation frequency, give the tmpressmn of contmously dlangmg spatml profiles, 6. PROPERTIES OF THE I M P L I C I T S I M U L A T I O N Some properties of the unpttcit simt,lation are briefly summarized m so far as they effect the application of :t method. a) Due to the exphc~t technique concermng the time mtegrabon (16), the s~muLabon loop can, depending upon the time and space increments, At and kz, become numer!c'_dly unstable. It can be shown both theoretically and experimentally that the stabdity conditions, ,?at ~ f(.XT), v,'did for the explicit difference method (for the numerical soiut.on of pdes ]61), have to Ee observed in the t,npDczt simulation [7]. These con&tions can be fulfilled on a hybrid compater due to the relatively simple realization of high computation frequencies b) The accuracy of the lmpl:clt simulatmn depends essentmlly on the spat,a[ &scretization of the mversion (IS), since the ,malog solution of the boundary value problem (12) and the time integration (16) are relatively accurate The number of space increments, 1,,, can not, however, be increased

Amudes de l'Assocu#zo~ mter~aHo~ale pore le Cazcul a,~alogique

1.74

--

N" 3

-

-

l/dlle t 1975

-

...... ]

OP

(k+l)TR

k TR

IC

~-

I I/O

n.~

rLEILrLFL. . . . . . . 21A I

,"" ">IZ ~ ~" ,>" . y ~ .>.

ANALOG

~,

~Xo ~N/ /\ \,\ '' ,

, I',+1 . \ , . "~ ~ ~ " . . . . . ~

~7//y7.

cu ,×

,,,,/,
' " ' \

\

"

\

"

\

"

L\ interpola tlon'-,~,', /

/+~

/

(~L;'I V/..,<

".~bc-'?

HYBRID

,/

/ /

-7/

<,,/v,.,//,,

", .., .., ~-'.., ..,., - , , \ N,

f

DIGITAL

/

/

•

/

,

.

V/,<;"

.

~,,,1 / / " / / . " , "

/

" .....

"'"""""

I I Fig. 4, - -

Seque~lce ol¢

zmpltcrt s ~ m u l a t m n

steps ( O P - O p e r a t e M o d e , IC - I n m a l C o n & t / o n Mode, I / O - I n p u t / O u t p u t - P u l s e s , T . - S m m l a t m n Cycle T u n e )

7

VIDEO

IK'1

INTERPRETER

¸

CL GP]

SL PB

/.V11: " I,-,,,.,,,. ,./y,,,>'///,> MACHINE ~ / / / / / / / / CODE , / . ~ ' / / //

If,,-

41,,-

===a,- DI SPLAY

POT INTERPRETER

DIGITAL Fig. 5 - -

KEYBOARD]

lIC

,.,ilil'-

DAC ADC

~I

DCFG

INTERFACE

ANALOG- LOGIC

[ n f o r m a t i o n flow of a h y b t k l mapInH s i m u l a t i o n (CL -Cont, ol Line, G P [ - I n t m r u p t , SL - Sense Line, PB - Pu.',hb u t t o n , D C F G - Dlg*t,lli)' C u n t m l l e d Functmn G e n m a t o r ) ,

W, Ruppel and M. Zettz' Dtstubated pa;amete;

175

arbitrarily, since the differentiating character of the reversion (iS') leads to new errors as a consequence of the noise of A/D-converted signals 2+ (&) 15] A further upper limit for the number 1,, is given by the hybrid computation time for the I/O-transfers and for the operations, H and Dt,-L if the simulation frequency is prescribed. c) The realization of the non-linear source function Q ( u , x ) and the transformation (12) take place in the continuous space domain. Due to this multiplex-techniques, pdes. with a separable source term can be simulated by the mephc~t method w~thout a large amount of eqmpment

"['he model is based on the simphfled assumption that the complex behavior m the reactor can approximately be concentrated in two p h a s e s ' The fluid phase (F) and the catalyst phase (K), both of which are connected by an intensive exchange of mass and heat TI
7. EXAMPLE

The ~/npl:c:t method is applied to the most important pde. (20) of the two-phase model to simulate Tic (z, t).

As an example from chemical engineering, a catalybc fixed-bed reactor (Fig. 6) is simulatcd by the Hnphcll method In the simplest case, a substance A is converted to a substance B in the reactor with the help of a catalyst K. The temperature and concentration profiles are to be simulated under the influence of both the feed varmbles and the wall temperature

Normahzation of the physical variables and introduction of the standardized operator (5a) yield the pde. of the normalized catalyst temperature, a = a (z, t)

catalys~ K .A

,~

v

-

. . . .

+

fi-'x

+

4- t,_, ( x -

Tw(z,I )

I

Tlz,t}

a,

~ t

-

& ;'~ (C'~,, .~)

T't..) 4- tg-'x

I

Z

(22a) (22b)

:, (z, O) := T'w (z, O)

B

,'"

TwlZ,t) ,,

?z ~

?.x ?.x z (0, I) = --~z ( t , l ) = O,

A

0

---

I

Fig. 6. - - Catalyt,c f,xed-bed reactor

T'v (z, t), T'w (z, t) and C',.. (z, t) represent the normalized temperatures and concentration, respect,rely The source term of (22a), Q (C'p, T ' v , x), contains the modified reaction rate r '~ (C'~,, v) caused by algebraic elimination of C~¢. The transformation equations 0,2) corresponding eo the dynamic boundary value problem (22) read

In [8, 9] the following two-phase model of a fixedbed reactor is proposed: ?, C~,.'~

~ Cv

~,t

~z

T 1"t I (P Cp)V

+

hi,,l; ((.l; -[~V

,. (;, c,,)~. 7 7

?.t

Cv).

c,, (o, t) = c,, (t)

4- "'"~ (T,¢ - - Tv) + ,,wr (Tw - - Tt..)

CK +}

,' (CK, T,,) + b,¢ v (C,.. - - C~,)

0 T.: (t' c,)i.

---?z

~t

--- .\i.

(0, t) --

3z

(18a)

(~sb)

T v(0,l) = T,,(t) 7J

07)

?'-' TI,: ....

~z-'

t- (-- a / , . ) ," (( ,~. T~.) -F a,.,,. (Tv - - T,.}

(2o.)

Tit(z, 0) = T w(z,0)

(20b)

(1, t) = 0,

r (C,¢, T,;) --- C,¢./?,,

09)

exp ( - - E/RT,~)

(2t)

/76

Anmdes

d "a wt &._, - - fi-' u't = du't

,,,, (o, t,,) -

de l'/ls.~ocmlton

Q(C'~.,T'r,.x)

(o, ¢)

=

nttermammde

(23a)

o

dz

d ~ we

---+ ,8 ~ ~'~ = O, dz'-' w~ (t, 6,-,)

,.,,.. (o, t,,) -

,

5. sin ,e

(23b)

du,..,

(o, t,,) =

dz

o

pozn /e Calcul m z a l o g l q u e -

N" 3 --

Parameter of standardization Number of space increments S,mulation cycle time T~me scale coefficient Digital storage

]udlet

1975

fi = 0,75 7r 28 T R = 7 . lO -a sec ,~ = 1.00 12 k t6-bit words. j. =

Some typical simulation results concerning the igmtion and blow-out of reactmn are shown in Fig. 7-8 Whereas the continuous profiles, T v (t, t,,) and C~.~(z, It,), result from the analog integration of quasi-stationary pdes (17)-(t8), the hnear interpolation of _TK (6,) yields polygonal profiles.

if the following setup is u s e d :

..~ (~, t,,) +

w (., t,.) =

,,,,, (~, t,.)

(23c)

A simple difference approximation of the Eq. 03) of the inverse operator G -~

1473 To 737~747°K

T oK

d-' y+ (z, z,.) dz'-' dy +

-dz

q- /3-' y+ (z, 6,) = y (z, 6,)

(0, t,,) --

d>,+

dz

(1, t,,) =

0

(24a)

t

(24b)

8Z .

yields the inversion matrix

2

h, h,

H = --

A Z~

0

h.., h,

h.,

o

/~-' A Z-' l0 t - -

t (24c)

2

'/,: h,

h,

~',- ~ . .

h._. = 2 h a =

t

2~ 0

O) ~

7 --b',-

0,003

The time integration (t0), in which the operator D, = at (~/~ t) is substituted, is realized by an explicit first order difference quotient At (Ii,+,)

--

j (t,,)

+

,: (tt,)

2._._..~_______..._.______

a 1

._,"(o) = T',,,

(25)

which is relatively exact due to tile small tune increments a t The time scale of the simulation can be manipulated by means of the integration factor of the vector equation (25) if the relation, At = c~Tit, between real time increments ~ t and the simulation cycle time T u (Fig. 4) is taken into consideration. The reactor pareazneters chosen for the n , p h c l t simulation of the quasi-stationary two-phase model can be found in [8, 91. Particular parameters of the hybrid simulation -v) are :

) The sunulat,on was implemented on the hvbtM computer system EAt PACER t00-680 of Instltut fiat Systemdynamik und Regelungstechnik

CF kmo[ m~

f

300s

0

0,5

z--,--

Fig. 7 - - Transient ptofdes of the reactor for ignmon of reaction Moreover, the ineplici! method has been used to simulate other phenomena of catalytic hxed-bed reactors known from experimental measurements or digital simulations. Multiple steady states of tile reactor, the influence of catalyst deactivation [9] and limit cycles of tile reactor in a closed loop have been simulated and recorded on film [10].

W. Ruppel and M. ZeJtz .' DJst,vbuted parameter

177

1473 To.

Z.

./7>--'-, S

oK

t"

1 873

II!

I

/"

I

1,1

is a feasible mstrument to both study the dynamic behavior of this distributed parameter system and devise control algorithms also implernentable on the hybr,d computer.

6 0 5 ~ GOO°K

)%.

/

III Ii

,/i

/"

ACKNOWLEDGEMENTS The authors gratefully acknowledge many helpful suggemons and comments made by Prof. E,D, Gilles concerning the research reported here. Th,s research was supported by the Deutsche Forschungsgemeinschaft.

I

//

r,

1/

l____

273

0

0,5

0

0,5

0,003

CF m3

t

z ---.~

1

Fig. 8 - - Transient profiles of the reactor for blow-out of reaction. 8. C O N C L U S I O N The implicit method of simaladng dlsLributed parameter systems governed by pdes. takes full advantage of the hybrid computer facilmes, since fast analog integrations and algebraic operatmns on the digital part are implemented simultaneously. Moreover, the presented hybrid simulation of the chemical reactor

[1] GILLES, E.D/ZEITZ, M : Implizltes Simulat,onsverfahren fiir 6rthch verteilte Systeme, 7.AICA-Congress, Prague 1973, Prepints I, pp. 191-194 [2] GILLES, E.D : Systeme mJt verteihen Parametern, Emfuhrung m die Regelungstheone, R. Oldenbourg-Verlag, MTiinchen.Wlen 1973. [ 3 ] V I C H N E V E T S K Y , R. : Hybrid methods for partial differentia!, equations. 6.AICA-Congress, Mb.nchen 1972 and SIMULATION 16 (1971), pp. 168-180. [4] LOBECK, B . Numerische Lfisungsmethode fiir mchtlineare partmlle modellgleichungen der chemischen Verfahrenstechu*k. Chemical Engineering Science 29 (1974), pp. 1320-t328. [5] RUPPEL, W. : Entwurf eines hybnden Programms zur Slmulatmn yon 6rthch verteilten Systemen. D,plomarbelt, University of Stuttgart 1974 (not published). [6] RICHTMYER, R.D./MORTON, K.W. , Dif[erellce Methods /or Initlal-Value-Ptoblems. Intersclence Publishers, New York-London-Sydney 1967. [7] ZEITZ, lV£, : Ein Verfahren zur hybriden Simulation yon 6rtlich verteilten Systemen Regelungstechnik 23 (1975) (to appear). [ 8 ] EIGENBERGER, G. ' On the dynamic behavlor of the catalytic fixed-bed reactor in the regmn of multiple steady states Chemical Engineering Science 27 (1972), pp. 1909-1924 [9] BLAUM, E : Zur Dynamik des katalytischen Festbettreaktors bm Katalysatordesaktivierung. Chemical Engineering Science 29 (!.974), pp. 2263-2277. [10] LOBECK, B/ZEITZ, M : Dynamik-katalyuscher Festbettreaktoren. 6. INTERKAMA-Congress, DiisseldorE 1974, Proceedings of VDI-Verlag, pp. 177-186; The corresponding t6-mm film is registered by Hochschulfllmreferat of the Free University of Berhn.