© Strojni{ki vestnik 48(2002)4,210-217                 © Journal of Mechanical Engineering 48(2002)4,210-217
ISSN 0039-2480                                                                                                                        ISSN 0039-2480
UDK 62-225.864:004.42:519.87                                                UDC 62-225.864:004.42:519.87
Izvirni znanstveni ~lanek (1.01)                                                                             Original scientific paper (1.01)
Ra~unsko  re{evanje  inverznega  problema oblikovanja  nadzvo~ne  {obe
A Numerical Solution to the Inverse Problem of Supersonic-Nozzle Design
Vinko Martinis - Branimir Matija{evi} - @eljko Tukovi}
Računsko oblikovanje nadzvočne sobe je občutljivo glede stabilnosti v področju nadzvočnega toka. Računski model, predstavljen v tem prispevku, se izogiba tej nestabilnosti z uvajanjem analitično določene porazdelitve tlaka na osi osnosimetrične nadzvočne sobe. Parametri toka v inverzno oblikovani sobi so preverjeni s programom FLUENT in prikazujejo enakomerno porazdelitev po prečnih prerezih vzdolž sobe. © 2002 Strojniški vestnik. Vse pravice pridržane. (Ključne besede: šobe nadzvočne, oblikovanje šob, modeli računski, problemi inverzni)
The numerical design of a nozzle is sensitive to stability in the region of supersonic flow. In the numerical algorithm presented in this paper the instability is avoided by the introduction of an analytically set pressure distribution on the axis of the axisymmetrical supersonic nozzle. The flow parameters of the inverse designed nozzle are checked by the application code FLUENT and they show a regular distribution on cross-sections along the nozzle.
© 2002 Journal of Mechanical Engineering. All rights reserved. (Keywords: supersonic nozzle, nozzle design, numerical solutions, inverse problems)
0 UVOD
Računski postopek oblikovanja nadzvočne šobe je še posebej občutljiv glede stabilnosti v področju nadzvočnega toka ([1] do [5]). Čeprav je uporaba računskih metod pogosta v praksi, le redko najdemo ustrezni algoritem v obliki uporabniškega programa za rešitev inverznega problema prenosa toplote in snovi. Problem je inverzen, ker je področje neznano [6].
V rešitvi, prikazani v tem prispevku, je v primeru osnosimetrične šobe določena izvirna analitična porazdelitev tlaka na osi simetrije. Izračun oblike šobe in parametrov toka je izveden po koračnem postopku po [7]. Začetni pogoj je izpeljan posebej. Parametri toka, v tako oblikovani šobi, so preverjeni z uporabo programa FLUENT. Dobljeni rezultati se dobro ujemajo.
1 OPIS MATEMATIČNEGA MODELA
Poleg kontinuitetnih, gibalnih in energijskih enačb, ki popisujejo tok v šobi, uvedemo funkcijo toka «P z izrazom:
0 INTRODUCTION
The numerical algorithm of supersonic-nozzle design is particularly sensitive to stability in the region of supersonic flow ([1] to [5]). Although the application of numerical methods is very common in practice it is very rarely possible to find an appropriate application code for the solution of the inverse heat- and mass-transfer problems. The problem is inverse because the domain is unknown [6].
In the solution presented in this paper, for the case of an axis-symmetrical nozzle, the original analytic pressure distribution on the axes of symmetry is defined. The numerical calculation of the nozzle form and the flow parameters were performed with the marching algorithm ac-cording to [7]. The initial condition was derived separately. The flow parameters in the nozzle designed in this way were checked for closed domain with the application code FLUENT and the results obtained correspond very well.
1 DESCRIPTION OF THE MATHEMATICAL MODEL
In addition to equations of continuity, motion and energy, which describe the flow in the nozzle, the stream function >F is introduced with the expression:
V<P = 0
VBgfFMK
stran 210
Martinis V. - Matija{evi} B. - Tukovi} @.: Ra~unsko re{evanje - A Numerical Solution
oziroma v valjnih koordinatah za osnosimetrični        i.e. in cylindrical coordinates for the asymmetrical
primer:                                                                         case:
d             dY
u  Yr + n^ dx            dr
kjer so: x - koordinata na osi simetrije; r - radialna koordinata pravokotna na os simetrije, u, v -projekcije hitrosti w na osi x in r . Sedaj enačbe toka prikažemo preprosto:
where: x is the coordinate in the axis of symmetry;
r is the radial coordinate perpendicular to the axis of
r symmetry; u , v are the projections of velocity w
on the axes x and r . Now the flow equations are
expressed in the simple manner:
(1) (2) (3) (4) (5),
dp        g dn dY        r dx	
dr2 =   2 dY     ru dr       n dx       u	
p = rg	
k + 1           2      J^ --------   -  --------   p k	-   n 2
1
kjer so: k = c p/cv - konstantni eksponent izentrope (za zrak k = 1,4); c, c - specifična toplota pri konstantnem tlaku oz p prostornini; r, p - gostota oz. tlak snovi v šobi.
Geometrijske količine  in hitrosti  so normirane:
where: k = c p/cv is the constant isentropic exponent (for air k = 1.4); c , c are the specific heat at constant pressure i.e. at constant volume; r, p is the fluid density i.e. the pressure in the nozzle.
The geometrical quantities and velocities are normalised:
(6)
kjer so s črto zgoraj označene dejanske količine; a* = (kp*/r*)1/2 – kritična hitrost širjenja nizkotlačnih motenj; p*, r* - vrednosti tlaka in gostote v kritični točki. Normirana funkcija toka je:
(7),
where the overbar denotes the real quantities; a* = (k p*/r*)1/2 is the critical propagation velocity of the low-pressure disturbance; p* and r* are the values of pressure and density at the critical position. The normalised stream function is:
r*a*r*
(8)
Normiranje prostorninskih koordinat s kritičnim r* ni nujno. V tem primeru je mogoče spremeniti merilo za posamezne količine.
2 ALGORITEM REŠEVANJA
Koračni algoritem. Za rešitev gornjega sistema enačb za vsak Y (Y = Yj) je treba določiti neznanke r, p, r, u in v (kot brezdimenzijske količine) za vsak x od vstopa do izstopa iz šobe. Področje je odprto v smeri pravokotno na os simetrije, problem je hiperboličen. Očitno je, da je za koračni postopek treba določiti začetne pogoje. Na osi simetrije je Y0 = 0, kar ne moremo upoštevati kot začetni pogoj, ki ga bomo določili kasneje.
Če so vrednosti neznank v vozlih pri Y = Yj določene, lahko določimo njihove vrednosti pri
Normalising of the space coordinates with critical r* is not necessary. In this case it is possible to change the measure for particular quantities.
2 ALGORITHM OF THE SOLUTION
Marching algorithm. For the solution of the above system of equations for every Y (Y = Yj) it is necessary to determine the unknowns r, p, r, u and v (as dimensionless quantities) for every x from the nozzle input to the nozzle output. The domain is open in the direction perpendicular to the axis of symmetry, and the problem is hyperbolic. It is evident that for the marching algorithm the initial conditions have to be set On the axis of symmetry y0 = 0, which cannot be taken as an initial condition, and will be determined below
If the values of the unknowns in nodes where Y = Yj have already been determined then their val-
| IgfinHŽslbJlIMlIgiCšD I stran 211
glTMDDC
Martinis V. - Matija{evi} B. - Tukovi} @.: Ra~unsko re{evanje - A Numerical Solution
>F=<Fj+1 prek diskretiziranih enačb drugega reda natančnosti. 1. Z diskretizacijo enačbe (2) dobimo:
ues on Y = Yj+1 are determined by means of discreted equations of the second order of accuracy as follows. 1.By discretisation of equation (2) is obtained:
ij+1
r +
kjer sta: x = i Dx, [l] - številka iteracije. 2. Z diskretizacijo enačbe (1) dobimo:
1 ru
l Pu
[l-1]
i, j+1
1
D<F
(9),
where: x = i Dx; [l] is the number of the iteration. 2.By discretisation of equation (1) is obtained:
[ l ]                     1
p\j+1 = p ij --g
3. Z diskretizacijo enačbe (5) dobimo:
1<9v r dx
1_dv r dx
[l-1]
i, j+1
D<F
(10).
3. By discretisation of equation (5) is obtained:
ij+1
4. Z diskretizacijo enačbe (4) dobimo:
5. Nazadnje je:
[l]
dr dx
k + 1
2      k -1 p k
1 + 11)            k-1
4. By
 1[ l ]
(11).
u ij+1

k + 1

5. Finally:
k-1     k-1
p k   -V
(12).
(13).
i, j+1
Algoritem je koračen po koordinati % vsak W = Wj lahko uporabimo kot obliko kanala šobe. Iterativen postopek je potreben, ker so količine na desni strani z indeksom i, j+1 neznane.
Začetni pogoji. Stabilnost rešitve dobimo z izvirno analitično določeno enačbo spremembe tlaka na osi šobe (^= 0):
The algorithm is marching along coordinate Y, and every Y = Yj can be used as a shape of the nozzle channel. The iterative procedure is necessary because the quanti-ties on the right-hand side with index i, j+1 are unknown.
Initial conditions. The stability of the solu-tion is obtained by an original analytical set equation of pressure change on the nozzle axis (Y = 0):
kjer so:
p(x,0) = ath(-bx) + c;     -2<x<3 where:
(14),
p= p
1
p*

a=-2 ( pin-pout )
in            out
Iz
k + 1
From
2 + (k-1)Ma2
kk-1
(15)
za mejne pogoje Main= 0,1 in Maout= 2,1 dobimo:
out
for boundary conditions Main= 0.1 and Maout= 2.1 are obtained:
pin = 1,880 pout= 0,207
(16).
out
Strmina krivulje za x = 0 preko koeficienta b > 0 lahko izberemo poljubno. Priporočeno je, da vzamemo b < 4. Pri manjšem b je sprememba tlaka v šobi počasnejša. V prejšnjih enačbah indeksi in in out označujejo vhodni oz. izhodni prerez šobe.
The slope of the curve for x = 0 over coefficient b > 0 can be chosen arbitrarily. It is recommended that b < 4 is taken. For smaller b the pressure change in the nozzle is slower. In the previous equations the subscripts in and out denote the input, i.e. the output section of the nozzle.
VBgfFMK
stran 212
Martinis V. - Matija{evi} B. - Tukovi} @.: Ra~unsko re{evanje - A Numerical Solution
Določenih pogojev p(x,0) na osi simetrije ne moremo neposredno upoštevati kot začetne pogoje. Vrednost neznank vzdolž tokovnice Y =Y1 izračunamo s povprečenjem njihovega razvoja v red oblike:
The set conditions p(x,0) on the axis of sym-metry p(x,0) cannot be used directly as initial condi-tions. The values of the unknowns on the stream line Y = Y1 are calculated by means of their development in series of the form:
NN
f(x,Y) = Yjfn(x)Yn +YXfn'(x)Yn
(17),
kjer je f (x, Y1) = r, p, v, r in n na Y = Y 1, ki je blizu osi simetrije.
Posamezno odvisnost spremenljivke za osnosimetrično šobo dobimo: a) Iz enačbe (2) izhaja:
where f (x,Y1) = r, p, v, r and n on Y = Y1, which is close to the axis of symmetry.
The single dependence of the variable for the axisymmetrical nozzle is obtained as follows. a) From equation (2):
— dY
rur— = 1 8Y
(18).
Za spremenljivke r, u, r vzamemo red f (r), f(u), f( r):         For variables r, u, r the order f (r), f(u), f( r) is taken:
1      N                   /N                                            (19).

n=0                                2Y    n=0                                      n=0
Po množenju se izenačijo koeficienti z enakimi                  After multiplication, the coefficients with
eksponenti spremenljivke Y na levi strani s svojimi        equal exponents of variable Y on the left-hand side
dvojniki na desni strani. Dobimo r in r .                        are equalised with their counterparts on the right-
hand side, and r and r are obtained. b) Iz enačbe (1) dobimo koeficienta p   in p   na        b) The coefficients pn and p  are obtained in a simi-
podoben način.
c) Iz enačbe (3) izhaja:
od koder določimo koeficienta v  in v  .
d) Iz enačbe (4) določimo rn in rn .    n’
e) Nazadnje iz enačbe (5) izhaja:
lar way from equation (1). c) From equation (3) follows:
dr dx
(20),
du
ru------- rv
dY         dY        gdY
from where the coefficients v and v  are determined.
d) From equation (4) rn and rn are determined.
e) Finally, from equation (5) follows:
dv        1  dp
Na osi simetrije je Y = 0, v = 0, v0 = 0. Iz enačbe (4) r0 = p01/k in iz enačbe (5) izhaja:
On the axis of symmetry Y = 0, v = 0, v0 = 0. From equation (4) r0 = p01/k and from equation (5) follows:
du           dv        1  dp
ru------- r v-----=------—
dY         dY        gdY
Posebej dobimo: rn(x) = 0 , pn’(x) = 0, vn(x) = 0, rn(x) = 0 in u n’ = 0. Glede na to so:
Specially obtained are: rn(x) = 0, pn’(x) = 0, vn(x) = 0, rn(x) = 0 and un’ = 0. Accordingly:
N
Posamezni koeficienti so:
r(x) =  Y1Žr:(x)Y1n+r0
n=0
N
p(x)=Yupn(x)Y1n
n=0
N
v(x) = Y1Yjv'n(x)Y1n
N
r(x) = Z rn (x)Y1n
n=0
N
u(x) = YJun(x)Y1n
The particular coefficients are:
(21).
gnn^nwiRaieKE
02-4
stran 213         |^BSSIrlMlGC
Martinis V. - Matija{evi} B. - Tukovi} @.: Ra~unsko re{evanje - A Numerical Solution
2
1/
\u0r0 J
,        1   / u1      r1
4    \u0    r0
r0 dx 1     1dv[- 1      r/
x dr'
p1 = -k~i—
r0 dx
p2 =-k~;— - p1—,
2    r, ox     2      r,
Ox
,        or'        or/
v1 = u1    0 + u0    1
dx        dx
(22).
k +1
2         k-1k
p0
2
----------p1 -—
kr0u0        2   u0
y
r0 = Pk
r1 =   p1
r0       kp0
Izračun koeficientov začnemo z znanim p(x,0) iz enačbe (14).
Koeficienti za dvoizmerno šobo so drugačni, dobimo pa jih na podoben način.
Prehod na realne parametre za določene robne pogoje na vhodu in izhodu je preprost.
Na opisani način dobljene krivulje Y = konst. (Y1=0,001 kot začetni pogoj in z iterativnim postopkom Y2= 0,005, Y3= 0,01 in Y4= 0,015) so prikazane sliki 1. V koordinatnem sistemu x, r za Y = 0,015 smo izračunali parametre u, v, r in r po enačbah (6) in (7) za tlak p, določen po enačbi (14). Prikazani so na sliki 2.
The calculation of the coefficients starts with the known p(x,0) from equation (14).
The coefficients for a two-dimensional nozzle are different, and they are obtained in a similar way.
The transition to the real parameters for the set boundary conditions on the input and output is simple.
In the described manner the obtained curves Y = const. (Y1= 0.001 as initial condition and by the iterative procedure Y2= 0.005, Y3= 0.01 and Y4= 0.015) are presented in Figure 1. In the coordinate system x, r for Y= 0.015 the parameters u, v, r and r were calculated by equations (6) and (7) for the set pressure p according to equation (14) and they are presented in Figure 2.
Sl. 1. Polmeri sobe za Y = konst, dobljeni z opisano metodo Fig. 1. The nozzle radii for Y = const. obtained by the described method
VBgfFMK
stran 214
Martinis V. - Matija{evi} B. - Tukovi} @.: Ra~unsko re{evanje - A Numerical Solution
Sl. 2. Normirani parametri u, v, r in r za Y = 0,015 Fig. 2. The normalised parameters u, v, r and r for Y= 0,015
= 0.1	M = 0.2	M =	1.0	M = 2.0	M =2.1
					
				777	
		W\\\	ta"		
-1.5
-0.5
0.5
1.5
Sl.3. Profil šobe in krivulje enakih hitrosti (Machovih števil) Fig. 3. The profile of the nozzle and the curves of constant velocities (i.e. of Mach numbers)
Na sliki 3 so za določene mejne pogoje in začetni pogoj prikazani profil šobe in krivulje enakih hitrosti (oz. Machovih števil).
S slike 2 vidimo, da so vse spremenljivke nespremenljive za -2<x < -1,5 in 2 < x<3 in da ta del šobe ni potreben (ravni del povzroči turbulentno mejno plast). Zato je na sliki 3 prikazana prostorska oblika šobe za -1,5 <x< 2.
Za določene mejne pogoje in dobljeno obliko r = r(x) za Y = 0,015 smo izvedli izračun z najbolj znanim programom FLUENT. Kot rezultat smo dobili normirane (bezdimenzijske) vrednosti spremenljivk p, u, v in r kot funkcije x.
Relativna razlika najvplivnejše spremenljivke p(x), izračunana z uporabo programa FLUENT po določeni enačbi (14), kaže zelo majhno odstopanje spremembe tlaka. Z naslednjim primerom je prikazano, da te razlike ne vplivajo na rezultate toka v kritičnem prerezu in na izhodu iz šobe. S tem je pokazano, da je predstavljena metoda natančna kakor FLUENT. Dejansko natančnost obeh izračunov je treba preveriti z natančnimi meritvami.
Primer. Izberemo: m = 0,5 kg/s; Ttot = 300 K; p ut = 101 300 Pa; Main = 0,1; Mao = 2,1; g = 1,4; R = 287 J/kgK. Z osnovnim izračunom za izentropni tok dobimo:
Na vstopnem prerezu: Ttot in = 299,4 K = Tott; ptot in = 926 400 Pa = pot ;
Finally, for set boundary and initial conditions the profile of the nozzle and the curves of constant velocities (i.e. of Mach numbers) are given in Fig. 3.
From Figure 2 it is evident that all variables are constant for -2 <x < -1.5 and 2 < x< 3 and this part of the nozzle is unnecessary (the flat part generates the turbulent boundary layer). Therefore, the three-dimensional shape of the nozzle is presented in Figure 3 for -1.5 <x<2.
For the set boundary conditions and the obtained shape r = r(x) for Y= 0.015 the calculation with the best-known application code FLUENT is performed. As a result, the normalised (non-dimensional) values of the variables p, u, v and r are obtained as a function of x.
The relative difference of the most influential variable p(x) from the calculation with FLUENT according to the set equation (14) shows a very small deviation of the pressure change. The example below illustrates that these differences do not influence the results of flow in the critical section and on the nozzle outlet. This shows that the proposed method is at the FLUENT accuracy level. The real accuracy of both calculations needs to be checked with a precise measurement.
Example. Known: m = 0.5 kg/s; Ttot = 300 K; p ut = 101 300 Pa; Main = 0.1; Mao = 2.1; g = 1.4; R = 287 J/kgK. By elementary calculation for isentropic flow the following values are obtained. At the input section: Ttot in = 299.4 K = Ttot; ptotin = 926 400 Pa = pot ;
^vmskmsmm
02-4
stran 215         |^BSSITIMIGC
Martinis V. - Matija{evi} B. - Tukovi} @.: Ra~unsko re{evanje - A Numerical Solution
V kritičnem prerezu:
Tthrot = 250 K; pthrot = 489 400 Pa; rthrot = 6,821 kg/m3;
uthrot = 316 938 m/s; rthrot = 8,58 mm;
Na izstopnem prerezu:
Tout = 159,4 K; uout = 531,46 m/s;
Za Y = 0,815 iz r* = rthrot = 0,00858 = 0,05 sta vstopni
in izstopni polmer:
rm
min
1,1735
In the critical section:
Tthrot = 250 K; pthrot = 489 400 Pa; rthrot = 6.821 kg/m3;
uthrot = 316 938 m/s; rthrot = 8.58 mm;
At the outlet section:
Tout = 159.4 K; uout = 531.46 m/s;
rthrot    0.00858 For Y = 0.815 from r* ==              = 0,05 the inlet
and outlet radii are:
rm
min
1.1735
rin =rinr* =0,418-0,05 =20,9 mm; rout= routr* =0,207 -0,05 = 10,35 mm in dolžina šobe: L =(1,5 + 2)0,05 = 175 mm.
3 SKLEP
V tem prispevku je prikazana stabilna inverzna računska metoda za oblikovanje nadzvočne šobe z uporabo izvirnega analitičnega izraza za padec tlaka na osi šobe. Na slikah je vidna pravilnost spremembe značilnih parametrov vzdolž šobe. Z uvedbo analitičnega izraza za spremembo tlaka lahko z izbiro koeficienta b v enačbi (14) oblikujemo tako daljše kakor krajše šobe. Raziskava je uporabna tudi za dvoizmerne in krožne šobe. Rezultati se dobro ujemajo z izračunom z uporabo programa FLUENT za izbrano obliko šobe.
rin = rin r* = 0.418-0.05 =20,9 mm; rout = routr* =0.207- 0.05 = 10,35 mm and the nozzle length: L =(1.5 + 2)0.05 =175mm
3 CONCLUSION
This paper presents a stable numerical design method for a supersonic nozzle by the introduction of an original analytical expression for the pressure drop on the axis of the nozzle. The figures show the regularity of the change of the characteristic parameters along the nozzle. Furthermore, by introducing an analytical expression for the pressure distribution it is easily possible, by selecting a coefficient b in equation (14), to design both a longer and a shorter nozzle. The research is also applicable to two-dimensional and annular nozzles. These results agree well with the calculation using the application code FLUENT for the obtained shape of the nozzle.
4 OZNAKE 4 NOMENCLATURE
približna kritična hitrost
dejanska oz. normirana dolžina šobe
Machovo število
masni pretok
dejanski oz. normirani tlak
posamična plinska konstanta, za zrak 260 J/kgK
dejanska oz. normirana radialna koordinata
kritični polmer
dejanska oz. normirana temperatura
dejanska oz. normirana komponenta hitrosti
w v smeri osi šobe dejanska oz. normirana komponenta hitrosti
w v smeri pravokotno na os šobe hitrost dejanska oz. normirana koordinata vzdolž osi
simetrije šobe dejanska oz. normirana koordinata pravokotno
na os šobe koeficient izentrope, za zrak 1,4 dejanska oz. normirana gostota dejanska oz. normirana funkcija tokovnice
Indeksi:
na vstopnem prerezu šobe
na izstopnem prerezu šobe
skupna vrednost
v vratu šobe
na osi simetrije šobe
a*        m/s      approximate critical velocity
L , L    m        real i.e. normalised length of nozzle
Ma                 Mach number,
m&         kg/s     mass flow rate
p , p    Pa       real, i.e. normalised pressure
R                    unit gas constant, for air 260 J/kgK
r , r     m        real, i.e. normalised radial coordinate
r*         m        critical radius
T , T   K        real, i.e. normalised temperature
s u , u   m/s     real, i.e normalised component of velocity w in
direction of the axis of the nozzle            s
v , v   m/s     real, i.e. normalised component of velocity w s                          perpendicular to the axis of the nozzle
w       m/s      velocity x , x   m       real, i.e. normalised coordinate along the axis of
symmetry of the nozzle y , y   m       real, i.e. normalised coordinate perpendicular to
the axis of the nozzle k                    isentropic exponent, for air 1,4
r , r   kg/m3 real, i.e. normalised density Y , Y kg/s    real, i.e. normalised stream function
Subscripts: in                   at input section of nozzle
out                 at output section of nozzle
tot                  total value
trot                 in throat of nozzle
0                     on axis of symmetry
VBgfFMK
stran 216
Martinis V. - Matija{evi} B. - Tukovi} @.: Ra~unsko re{evanje - A Numerical Solution
[1]
[2]
[3]
[4]
[5] [6] [7]
5 LITERATURA 5 REFERENCES
Minkowycz, W.J., E.M. Sparrow, G.E. Schneider and R.H. Pletcher (1988) Numerical heat transfer, Wiley,
New York,.787-834.
Shimura, K., Y. Asako and J.H. Lee (1994) Numerical-analysis for supersonic flows in a cooled nozzle,
Numerical Heat Transfer 26, 631-641.
Cler, D., M. Lamb, S. Farokhi, R. Taghavi and R. Hazlewood (1996) Application of pressure-sensitive paint
in supersonic nozzle research, Journal of Aircraft 33, 1109-1114.
Demetriades, A. (1996) Stabilization of a nozzle boundary layer by local surface heating, AIAA Journal 34,
2490-2495.
Tsunoda, K., T. Suzuki and T. Asai (2000) , Journal of Fluids Engineering 122, 585-591.
Anderson, J.D. (1995) Computational fluid dynamics, McGraw-Hill, New York, 772.
Pirumov, U.G. and G.S. Rosljakov (1978) Gas flow in nozzles (in Russian), Moskovskij univerzitet, Moskva,
351.
Naslov avtorjev: doc.dr. Vinko Martinis
prof.dr. Branimir Matijaševič
mag.Željko Tukovič
Zavod za energetska postrojenja
Fakultet strojarstva i brodogradnje
Sveučilište u Zagrebu
Ivana Lučiča 5
HR-10000 Zagreb, Hrvaška
Authors’ Address:    Doc.Dr. Vinko Martinis
Prof.Dr. Branimir Matijaševič Mag.Željko Tukovič Department of Power Eng. Faculty of Mechanical Eng. and Naval Architecture University of Zagreb Ivana Lučiča 5 HR-10000 Zagreb, Croatia
Prejeto: Received:
15.3.2001
Sprejeto: Accepted:
23.5.2002