. . . PU.BUCAÇÃO DA ABCM • ASSOCIAÇÃO BRASILEIRA DE CIÊNCIAS MECÂNICAS

... . . V O L. XVIII • No. 4 • DECEMBER 1996 ISSN 0100-7386

)OURNAL OF THE ~RAZILIAN WCIETY OF MECHANICAL ~ClENCE~

RrVJm ~RA~IWRA DE CI[NCJA~ M~CÂNJCA~ EDITOR: Leonardo Goldsteln Júnior UNICAMP- FEM- OETF - C.P &122 13M3 970 Ccmpon~s SP Tet. (0191 239 -JoJOô Fax. (01~;. 239 -ó/22

EDITORES ASSOCIADOS: Agenor óe Toledo Fleury !PT - tustotutc co Pcsqu·sas TiKr.~t~gttas Dtv•sao d~ Mec~111ca e Eteln~•Cade - Agrupamento de SisleiT'as ce Gür•t'o!e Cidade Uouws•tária · C.f' 7141 -)1~~4 -910 Sao Paute ~P ret : 10111 263 -221 1 R<!rr.at 5oJ4 Fa<· (0 111 8ü9-3353

Angela Ourivlo Nleckele Po~titíc ia uni·,ersidaóc CatCiica tlo ~,o de Janw•j Departamento de Engenharia Mecânica ~ua Marquês de S~o llitente. 225 Ga·,e~ 22453-900 Ri:> óJP. Janeio o RJ T~t · (0211 2::\9-0719 Fax (~·m 294-gHe

Carlos Atbeflo Carrasco Allennani UNICAMP • fEM · OE - C .I'. 5122 13083-970 Camporo3S SP Te!· (019) 239-843& Fax (019i 239 3i22

Paulo Eigl Mlyagl Univ~rsldade d~ São Paul3 - E>tola p,):iléc.1oca Depanamento de Eng~nha11~ Mecanica - MecdlrCnic• .A.ven•da P101 Mcllc· MJ1aes. 2231 a~5·JB·900 sao Paute SP fel· (0 11 ) B1e-5580 Fax· 1011 ) 813-5471i8 13- 1886

Walter L. Weingaertner Uni;e•s•dade Federal de Santa Catarina Departamento de Eng~nhada Mecânota - Labcratóll•) d& Mecán•ca de Pr~cisu, Campus Unw!rsnalic Tlindade C.P 476 88·J40-902 Florianópoll; SC Tet . !~42) 231-9395!2:l4-5277 Fa).: (n4R) 234-1519

CORPO EDITORIAL:

AICII ~e Fa•o Oolaudo !Pl!C·RJ) .AJH~)IIIO F1anc1sto For1es 1Ur8·1 Armanjo Albtrtau• Jr. (UFSC) . Alair RIOS Neto (INPE) Be~e-Jito Moraes PurQuoiro <EfSC t,SPi Caio Marb Cos1a (EMBRACC) Ca1los Alberto de Almeida(PUC-RJ) Gartos ~.IMrt (l Manin tlJFSC! Ch:v•s Ra:mundü Mah~ka (UfSC; Ema""et Rochd Wc•stl 1 UNESP-TEtS 1 r ranc1sc~ Emiõi~ Baccaro N•·J'O i!Pl-SP: Francis.:c JJ~ 51m~es (UFPbj Genes•o J))é Mencn !E FEl) Hans lngo Weber (UNICAMPi Hcrulquc Rcz~niO!d !EESC USP) Jair Carlos Outra (UFSC) Jo3u All.i•<o HeoL de Jo•nada (llfAGSi Jose JD3o j e ~s~tndoi<l Wf SC) Ju1andlr 11om Yanag111ara (EP USP) Lirio scnaet~r (UFRGS"t Lounvat 6oens (UFsc·, Luõs Carlos Sanj~vat Goes (!TA) Mar~ to Zh·lani ~ UFMGi Mnjses l•M•Iut. ;ÇOPPE-IJFRJi Ni~>O !]e carvalhO l vbG BIUIII (C:tPPE-Uf~J) Nlvat;Jo Lemo~ Cupim tUN!C.".Mi') Pauto Atcnso de O r;wa SO'IIelO (ITA) Pcgeric Marto~s Sa>oanha ja Gama <LNCC) Va oer Ste:tcn Jr (UFU)

REVISTA FINANCIADA COM RECURSOS 00

Programa de Apoio a Publicações Científicas

MCT @l CNPq [!} FINEP

RBCM - J. of the Braz.. Soe. Mechanlcel Sclenc.s Vol. XVIII- No.4- 1996 - pp. 308-317

ISSN 0100.7386 Printed ln Brazil

Position Control of Mechanical Manipulators Using the Optimal Root Locus

José Jaime da Cruz Carlos Alberto Matuoka Universidade de Slo Paulo Escola Politécnica - c p 61 548 Departamento de Engenharia Elétrica 0542<4-970 São Paulo, SP Brasil

Abstract

The design of a robust joint independenl positioo cootrol scheme for mechaoical manipulators using a pote placement approach based oo the Optimal Root Locus is discussed in lhis paper. lts perfonnance is compared to that of the traditional PD structun:, widely used in industrial applications. in situatioos \Were the last ooe is not recommended. The in1luence ofvarious feedforward schemes is also analysed. Keywords: Position Control, Mechanical Manipulators, Pole Placement, Optimal Root Locus, Linear Quadratic Regulator.

lntroduction

Thc objec1ivc of this papcr is to design a robust position control scheme for mechanical manipulators by using a pole placement technique based on the Optimal Root Locus Method (ORLM). The controller structure is joint lndependent (Spong and Vidyasagar, 1989).

The ORLM is a particular way of designing a Linear Quadratic Regulator (LQR) such that the closed-loop potes of the system can be allocated in pre-specified locstions of the s-plane (Thompson, 1980). The use of the LQR as a basis for lhe design is motivated by its robustness characteristics: large stability margins, robustness in face of nonlinearities and tolerance to delays (Anderson and Moore. 1971).

A comparison of its perfonnance with the traditional PD scherne is dane, partícularly in conditions where the use of this controller is not recommended. namely, when low gear ratlo andlor high speeds are present.

The PD controllcr design is perfonned based on a linear dynamical modcl of the manipulator. An upper bound on the system closed-loop natural frequencies at 50% of the joint structural resonant frequencies is imposed (Luh, Walker and Paul, 1980).

T he flexibility of the actuators axes is included in the ORLM design modcl to a llow a larger bandwidth for the closed-loop system.

The perfonnance comparison mentioned above is done through digital simulation. A planar two degree-<Jf-freedom manipulator with revolute joints is used for this purpose. The simulation model is distinct from the nominal one in that the fonner includes the following aspects : existence of an unk:nown load at the end-effector, Coulomb friction at the joints and random errors in thc parameters of the dynamical model.

The influence of various feedforward schemes in the global perfonnance of Lhe control system is also studied.

Mathematical Model

The simulation model will be discussed next since it contains as a part icular case the nominal (design) model.

The equations of the joint torques have been obtained by the Newton-Euler method (Spong and Vidyasagar, 1989)

Manuscript received: December 1995. Technical Editor: Agenor de Toledo Fleury

309

and

J. of the &az. Soe. Mechanical Sc:iencet - Vol. 18, December 1996

[ 2 2 2

m1d 1 +J 1+m211 +m2 1 1 ~cos(+52 - t, 1 ) +M11 +

Ml 112 cos(+,2 - +, 1) l+aJ+

[Mti +MI112 cos(+,2-+, 1) +m211d2 cos(+52 -+11 ) + 2 ..

m2d2 + J2 I +s2 +

[MI 112sin(+,2-+,1) + m211d2 sin <+,2 -+11 )] <+;, - +~2) +

(m 1d 1 + m211 + Ml1) gcos (+, 1) + (m2 +2M) d2 gcos (t,2)

2 -(Mh + m2d~ + J2) +s2 +

. 2 (MI 112 +m211d2)sin(+12 -+,1l+s1 +

(I)

(2)

where, for i= l,2 -r1 is the torque applied to joint i, +,1 is the angle of link i with the horizontal (Fig. 1) and, for each link i, ffij is its mass, IJ is its length, l j is its moment of inertia wíth respect to the center of mass, dj is the distance from its center of mass to joint i and M is tbe mass of the load at the end­effector. Ali these symbols correspond to the simulation model ("real" manipulator); g is the local gravity acceleration.

Fig. 1 Scheme of the TWo o.gr.e-of-frMdom Manlpullltor

Those terms containing second order time derivatives of the joint angles correspond to inertial torques; those ones with first order derivatives are associated to the centrifugai effects; the terms that depend directly on the joint angles represent the gravitational torques.

Notice that lhe Coriolis torques do not appear in this formulation since lhe angles are measured in an inertial frarne.

Let Ti, i• l,2, denote the torque on joint i, excluding the inertialtorque corresponding to link i, i.e.

T. = t - J .~ . I I I SI (3)

J.J. da Cruz et ai.: Position Control of Mecha nicai Manlpulators Using the ... 310

Each actuator is assumed to be a DC motor coupled to a gear train through a flexible axis, as shown in Fig. 2. The motor torques 1;,. 1, ix l,2 are takc:n as the control variables.

The dynarnical model of joint i can be written as:

(4)

- . Ii+ si+ 8 i+si- (Km/n;) <+mi- +s/ 0 ;) + T; = O (5)

where, for each actuator i, Jmi is the moment of inertia, + . is the angular displacement, Bmi is the viscous friction coefficient, Km; is lhe torsional spring con~t ofthe axis, Cm; is the Coulomb friction torque, n, is the gear ratio and B; is the viscous friction coefficient corresponding to the link i. sgn(.) is the sign function.

DCMotor J,..

( 1\ \ Klrf

li r----

~ 77/77777 ~. 1'1

" ~ c;,..,.,.,.. Blri. Qn

777 8i

l:n

Fig. 2 Sc:heme of ln Ac:tu1tor

Notice that Eqs. 4 and 5 express lhe condilions of ( dynamical) balance of lhe torque corresponding, respectively, to the motor and link sides of joint i.

The nominal model can be obtained from the simulalion model above after elimination of lhe terms tbat contain the mass of the load ai the end-effector (since its is admitted as unknown in lhe design) as well as those which represent the Coulomb friclion effects. Furthermore, lhe nominal variables are denoted with the sarne symbols as above but with an over bar. For example, -ti and m; are respectively tbe nominal torque at joint i and the nominal mass of link i.

ln Lu h, Walker and Paul, 1980, the PD controller design is based on the assumption that the actuator axis isinfmitely rigid (in fact tho flexibility of lhe axis just imposes a constraint on the bandwidth ofthe system). This assumption means that

(6)

Pole Placement Using the ORLM

Consider the LQR that results from the minimization of

(7)

subject to

x(t) = Ax(t) +bu(t), (8)

y (t) = hx (t), (9)

311 J . of the Braz. Soe. Mechanical Scien<:es • Vol. 18, Oecember 1996

where x (t) e 9t", u', y (t) e 9l, being lhe remaining dimensions compatible and r> O. As usual , assume that (A,b) and (h,A) are controllable and observable respectively.

It is well known (see, e.g. , Kwakemaak and Sivan, 1972) that the solution to this problem has the state feedback forro

u(t) = - kx(t) (lO)

A procedure to choose r and h sucb tbat the locations of the closed-loop potes are pre-specifíed is presented in the following.

Let the polynomials p

n (s) = b n (s -z.) p I (li) i • I

and

n

d (s) n (s-p;). ( 12) i • I

where p < n, be defined by

~ = h(si - A) -1b

d (s)

and tet D(s) be defined as

l O ( s) - d ( s) d ( -s) + r n ( s) n ( -s)

(13)

(14)

The closed-toop poJes P.·; (A - bk)) are the left half plane roots of D(s) (Kailath, 1980).

~uation D(s) =O can be rewrinen as

n - p - tn~ (s+z,)(s-z.,) ( - J}' 1 a l = - )

n~ = 1 (s+ p1) (s - p;) ( 15)

Hence for varying r the closed-loop poJes can be obtained using the classical root locus techniques (Ogata, 1982). with the parameter (- I t- P, -

1 playing the role ofthe gain.

ln this way, given n(s)/d(s), we take its poJes and zeros as well as their symmetrical ones with respect to the imaginary axis and then draw the root locus following the classical rules. The left half plane portion ofthe ptot represents the LQR potes.

When r-+ O~ • the roots of D(s) that remain finite tend to the zeros of n(s)n( -s); those ones that go to co , follow a Butterworth pattem (Kailath, 1980)

Since the poJes Pi correspond to the eigenvalues ofthe open-loop system matrix A, they are tixed. On the other hand the zeros z.; dcpend on vector h which can be considered as a design parameter.

ln this way, once lhe p desired closed-loop potes (z, ~ ..... z,) have been chosen, lhe poJe placement procedure using the ORLM can be reduced to lhe determination of the vector h which can be summarized as (Thompson, t980):

• if(A,b) is in the cannonical controllable forro, h is given by

J.J. da Cruz et III.: Position Control cA Mechanical Manipulators Uslng lhe ...

h = [b0 b 1 ... bP].

where b; is the i-th degree coefficient ofn(s);

• otherwise,

[ n - p - 1 J - 1 h = [0 ... 01] x 1 ... xP b Ab ... A b .

where

-I X; = (z;I-A) b

Controller Designs

PD Controller Deslgn

(i = l ,2, .. . ,p)

312

(16)

( 17)

( 18)

The block diagram ofthe control system for the i-tbjoint is shown in Fig. 3 (Luh, Walker and Paul, 1980) where +si is the reference signal, ~ is the gain of the tachometer that measures the angular velocity of the motor axis, K+i and K1; are respectively the proportional and derivative gains of the controller, K1; is the torque constant of the motor, R; is the annature resistance, Kbi is the back e lectromotive force constant and T mi is the feedforward torque.

From Fig. 3 tbe actual motor torque can be written as

K1. [ - · J 2 · Tmi = R;' K+;(+1; - +1;) -(Kb;+K1;K1;)+mi -n;t1-ni B;+m;+T~; (19)

K+i and K 1 i are usually set such that: i) the closed-loop poJes must be such that the natural undarnped frequency must be less than o r equal to 50% of the structural resonant frequency of lhe actuator axis; ü) the closed-loop poJes must be critically damped. These conditions lead to (for details see Luh et ai., 1980):

K .R. K ,. m1 1 .. .;;,--

' 4n;K1; (20)

and

(21)

313 J . of the Braz. Soe. Mechanlcal Sciences • Vol. 18, December 1996

T he feddforward torque T:.,; is evaluated using the link generalized coordinates and is based on pre-calculated estimates ofthe inertia, gravitational, Coriolis, centrifugal and viscous torques (Luh et ai., 1980).

Controller Design Using the ORLM The following state vector has been adopted in this case:

(22)

The last element of the vector represents an integration included to guarantee that the steady-state position of the i-th link coincides with the nominal (desired) one for a step command.

The control variable has been taken as

(i) -u = Tm;-Tm;. (23)

which has the natun:: of a correction torque.

With the above choices and considering T; as a disturbance in the nominal version of Eqs. 4 and 5 it is easy to obtain the followi ng state model :

where

o o o o o

Ái) = - Kmifimi Km/ (n;Ím;) - Bm;li01 ;

Km;/ (n;j;) - 2-- Km;/ (n; J;) o

o I o

o o o

o o -B;I i; o

o o

(24)

o o

1/Ími

o o

(25)

Then a set of p = 4 closed-loop potes is chosen on the negative real axis from 1.5 to 2.0 times faster than those corresponding to the PD controller in order to include the structural resonant frequency of the actuator axis in the bandwidth of the closed-loop system. Finally the methodology described in section 3 is applied to find a gain vector k such that

(26)

Notice that in the proposed scheme (see Fig. 4) the nominal motor torque T mi represents a feedforward control signal. According to the nominal model described, T mi is given by

(27)

J .J . da Cruz et ai.: Posltion Control of Mechanlcal Manlpulators Uslng the ... 314

, G oFORWARo c-1_ I - _j T,.;

x<•>- ----i!Ç .. J~ L;, ~ ~ I - MANIPULATOR llji> x(i)

~ J U~ _j I I I L

Fig. 4 Bloek Olagram of tha ORlM Controller.

ln the position control problem lhe function .si(t) is given. To compute T mi it should be noticed lhat

- n. - - - "' - -tmi = ::2- (Jitsi + B;tsi +Ti) + +si/ni

Km i (28}

and that both tmi and +mi can be evaluated analytically by successively differentiating the last expression wilh respect to time.

ln view of Eqs. 1-2 it is evident that the computation of tiand t; is an extremely laborious task (even in lhe case of a two degree-of-freedom manipulator) (Matuoka, 1993).

For the reasons exposed above thls scheme is called here full analytical feeedforward.

Motivated by tbis structure and by the computational burden to evaluate T mi, tbree altemative schemes have been tested:

• Full numerical feedforward

An efficient way of evaluating ;mi and fmi is lhrough a numerical approach. An approximate differentiatioo scheme based on the Newton interpolation formula (Demidovich and Baron. 1973) has been uscd.

Recall that the functions to be differentiated have not to be measured but just evaluated numerically being thus not affected by significant noise. For this reason the numerical derivatives match closely the analytical ones (Matuoka, 1993).

• Simplified feedforward

ln lhis case the axis is considered to be infinitely rigid. Hence, given +s1(t), we compute - -+mi '"' +s;/ni (29)

and

(30)

• No feedforward

This is the simplest case and corresponds obviously to take

-Tmi = O (31)

315 J . of lhe Braz. Soe. Mechanlcal Sciences - Vol. 18, Decembef 1996

Results

The v alues of the nominal model parameters are shown in Table I .

Table 1 Parameters of the Nominal Model

Para meter Value Para meter V alue

11.2 0.75m d,.2 0.375 m

m,.2 3.0 kg .ltt..z 0.5625 kg.m2

Jm1.2 0.0233 kg.m2 R,,2 0.91 n Kct,2 2.0 Nm/A Bm1.2 8.09E-05 Nms/rad

Kct.2 0.05062 Vslrad Kc.t.2 0.10123Vslrad

B,,2 8.09E-05 Nmslrad K.nt.2 574.9 Nmlrad

n,,2 0.5

The values ofthe parameters ofthe simulation model have been obtained by adding a 5o/o pseudo­random error to the corr~ponding nominal ones (except n1.:J. Furthermore the load at the end-effector has been taken as M=l .O kg and the Coulomb friction magnitude Cmi• 0.5% of the maximum gravitational torque at joint i.

The value of the gear ratio (2: I) has been chosen lmentionally low (in practice, it is usually of the order of 100: 1) to permit the comparison of controller performances in cases where joint independent control is not recommended {Spong and Vidyasagar, 1989). Obviously in a concrete design problem the use of such a gear ratio would require the choice of the actuators to be reviewed. This is not the case here since we are only interested in the comparison of controllers performance.

The PD controller design produced a double pole at -75s·1 for the system of Fig. 3. ln view of this the design parameters z, ofthe ORLM have been chosen as {-100; -110; -120; -IJO}(s-1) . Joint 1 is dynamically more complex and bence harder to control. For this reason only the results corresponding to it are presented in the following.

The simulation has been carried out for a position reference signal composed of a pai r of parabollic ares corresponding to two accelerat ion pulses with duration 0.5 sec and amplitude + n/ 9 and - n/ 9rd/ sec1• respectively. The ORLM controUer has been simulated with no feedforward.

Figures S to 8 sbow the results obtained. lt can be noted that the position errors are of the order of I ()3 times smaller for the ORLM controller, although the control cffort has similar magnitudes for both controllers. Hence the pcrformance of the ORLM controller is significantly better than that of the PD controller.

o, ---­'\ \

· O.~ ·\

I I

U nlc t

1 -o.o• , 1

j \, /-~--...... __ ,--------~ :Lo, ~ l 1

"'----- - • j I i I ---·

\I \../

-002

·0.02~

-0 .03 0 0.2 0.$

-~·-~---' I 4 1.$ 1.1 2

Fig. 5 Posltlon Em>r - PO Conttoller

J.J. da Cruz et ai.: Position Control of Mechanical Manipulators Using the ...

Uni< I 21.--~-~--.---,--~-~--.---,--~--.

2~,L--02--D-~--~,--O-.I--1--1~2-~1.4--1-.6-~1 .1-~ r .... (•ool

Fig. 6 Total Torque (continuou• Une) and Feedfoi'Wllrd Torque (daahed Une) - PD Control .. r

.,l /\/', ,. I '

' I' I - ;J i ur ( I

l i I , o~ I -------~

·W.~.2-0~,.~....,.~.-.~ .. -~.--.-.2-.~.-~ .. .,... ~ .. ~ T'me(,.::)

Fig. 7 Poaltlon Error - ORLM Controller

linlrl

__~-

lO

·.~-.::':.2,---:: ... .,.--=o~ .• -""o .• =---,c---,:'::.2-.,..,1. •. _ u- ,:.-4

T""(•t<l

Fig. 8 Total Torque - ORLM Controller

316

317 J. of ttle Braz. Soe- Mechanieal Soenoes- Vot 18. December 1996

Furthermore the tests showed that a good perfonnance of the PD controllcr required the presence of a feedforward tenn.

Ali feedforward schemcs described in have been tested for the ORLM controller. The results showed that the feedforward terms could be eliminated since they did not improve the performance when compared to the case with no feedforward. This fact has thus verifted in practice the robustness of the LQR Íl1 face of model uncertainties and nonlincarities.

Both full analytical and numcrical computation schemes gave practically identical results for the feedforward tenn. Hcnce, ifthe feedforward scheme had been required. the numerical approach should be takcn in view of its simplicity, particularly in the case of manipulators with a large number of degrees of freedom.

ln the simplified fccdforward scheme the magnitude ofthc fcedforward torque is relatively small when compared to Lhe full feedforward scheme. The axis inftnite rigidity is the characteristic of the simpli fied scheme responsible for this behavior. So under this point of view we can say that the simplificd feedforward schcme is closer to the the case with no feedforward.

Tests were also perfonned for the case where the end-effector describes a spiraJ trajectory. Results obtained confinned the conclusions above (Matuoka, 1993).

Conclusions

The results obtained with thc ORLM controller indicate that it must be considered as an altemative to the solution of tbe position control problem in cases whcre thc PD scheme is not appropriate. Low gear ratio and high speed motion ofthe links are typical situations where this occurs.

Since the ORLM controller has the structure of a state feedback, both position and velocity must be measured for each actuator axis and each link axis. One altemative currently under investi~ation considers the use of a Kalman Filter in such a way as to recover the original robustness properties of the LQR as in the LQG!L TR procedurc (Doyle and Stein, 1981).

With respect to the simulations performcd, it should be emphasized that the test conditions can be considered severe with respect to modeling errors. The simulation model (but not the nominal model) included a load at lhe end-cffcctor, Coulomb friction at the actuators as well as random errors in the model parameters. The tests have clearly confirmed in practice the well-known robustness charactcristics of the LQR.

Acknowledgments

The first author wishcs to express h is gratitude to CNPq -Conselho Nacional de Desenvolvimento Cientifico e Tecnológico (Proc. No. 304071/85-4(NV)) and to FAPESP- Fundação de Amparo à Pesquisa do Estado de Sao Paulo (Proc. No. 91/0508-3).

References

Anderson, B.D.O., and Moore, J.B., 1971 "Linear Oprima! Control~, Prentice-Hall, Englcwood Cliffs, NJ. Demidovich B.P .. and Maron, I.A., 1973, "Computational Mathematics", MLR, Moscow. Doyle, J .C and Stein, G .. feb 1981, ''Multivariable Feedback Design: Concepts for a Classicai/Modem Synthesis"

. IEEE Trans. Auto. Contr .. Vol. AC-26. Kailath, T., 1980, "Linear Systcms", Prentice-Hall, Englewood Cliffs, NJ. Kwakemaak H. and Sivan R .. 1972, "Linear Optimal Control", Wiley, New York, NY. Luh, J.Y.S., Walker, M.W. and Paul, R.P., 1980, "On-Line Computational Scheme for Mechanical Maniputators•.

Trans. ASME J. Dyn. Syst. Meas. Contr. Mai\Joka, C.A., 1993, "Controle de Posiç.Ao de Manipuladores Mecânicos Usando o Mttodo do Lugar das Raizes

Otimo", Tese de Mestrado, EPUSP, Sao Paulo, SP. Ogata, K., 1982, "Engenharia de Controle Moderno", Prentiee Hall do Brasil, Rio de Janeiro, RJ. Paul, R.P., 1981, "Robot Manipulators: Mathematics, Prograrnming and Control", MlT Press, Cambridge, MA. Spong. M.W. and Vidyasagar, M., 1989, "Robot Dynamics and Control", Wiley Jntemational, New York, NV. Thompson, P.M., 1980, "Multivariable Control Systems -The Optimal Root Locus (Lecture Notes)", MIT,

Cambridge, MA.

RBCM - J . of tha Bru. Soe. Mec:hantc.l S<:lanc. VoJ. XVIII - No. 4 - 199tl - pp. 318-329

Dinâmica do Rodeiro Ferroviário Railway Wheelset Dynamics

Roberto Splnola Barbosa IPT- Instituto de Pesquisas Tecnológicas Divlslo de Tecnologia de Transportes- C.P. 7141 01064-970 5ao Paulo, SP Brasil

Alvaro Costa Neto USP • Universidade de São Paulo Escola de Engenharia de São Carlos Departamento de Engenharia Mecânica · C.P. 359 1356D-970 SAo Carlos, SP Brasil

Abstract

ISSN 0100.7386 Prlnted ln Brazil

A velocity parametrized railway wbeelset model was produccd for dynamic behavior analysis. Lateral wbeelset excursion reJated to the track is the first vibration mode. lt was observed in the eigen-properties that the wave length ofthis movement is approximately constant and dependent on the wheel conicity. Modal damping ofthe first mode is inverse ly proportional to velocity and becames negative for high speed, resulting on system instability. At low speed, real, distinct and overdamped second mode eigenvalues are inversely proportional to veloc ity and strongly coupled with contact stiffiless. Dynamic wheelset behavior through a variable track trajetory may be observed during model simulation. Keywords: Railway Wbeelset, Dynamic Behavior Analysis, Simulation

Resumo Foi elaborado um modelo linear do rodeiro ferroviário parametrizado em funçAo da velocidade para análise de propriedades e avaliaçlo do comportamento dinâmico. O primeiro modo de vibrar corresponde ao movimento de passeio latttal do rodeiro em relaçAo a via. Observa-se nas autopropriedades que este modo possui comprimento de onda de movimento no espaço aproximadamente constante e fortemente dependente da conicidade da pista de rolamento da roda. O amortecimento modal é inversamente proporcional à velocidade e fica negativo para velocidades elevadas tomando o sistema instável. O segundo modo de vibrar, sobreamortecido, possui, a baixas velocidades, ralz.es reais distintas inversamente proporcionais à velocidade e fortemente acoplado com a rigidez das forças de contato. A simulação deste modelo permite visualizar o comportamento dinâmico do rodeiro em sua interação com a via férrea para um percurso com trajetória variável. Palavru-<bave: Rodeiro Ferroviário, Comportamento Dinlmico, Modelo linear.

Nomenclatura n oroem elo SIStema; t.t = intervalo de ~lculo; {I} = vetor de força externa; Uy = coordenada de a angulo de kid(; {v} autovetor do sistema 2n:

deslocamento lateral; T = Intervalo de tempo: (M) matriz de massa;

•• c coordenada angular, t.. = tempo no k41ssimo (C) = matriz de amortecimento; m = massa do rodeiro: intervalo de tempo T; [K) matriz de rigidez; El momento de Inercia; c; "' fator de amortecimento; (R) = matriz modal; r .. = Força de Contato; ÔL = comprimento de onda; (AI = matriz dinamica do s. "' Velocidades relativas; ).., = i-ésslmo autovalor, sistema; À = conicidade da pista de (x} = representaJ!o vetorial (I) = matriz identidade;

rolamento da roda; de estado e espaço; [<J>) = r vetor de coordenadas; (x}' = vetor de estado matriz Fundamental;

V o velocidade de transposto; translaçjo; { x } = 1• derivada do vetor de (A) = matriz diaponal da

= forçamento externo; estado; exponenc1al dos = tempo; {u} = vetor de forçamento autovalores;

externo;

Manuaalpt received: January 1996. Technical Editor: Agenor de Toledo Fleury

319 J . or lhe Braz Soe. Mec:Mnlcal Scienoes • Vol. 18, December 1996

Introdução A inscrição de um rodeiro ferroviário em curvas se faz com o auxilio de um inteligente sistema

dinâmico estabelecido pela conicidade da pista de rolamento das rodas que produzem diferentes raios de rolamento para cada roda em função do deslocamento lateral do rodeiro.

Um observador trafegando junto ao rodeiro em uma trajetória retilinea percebe que a via férrea se desloca lateralmente quando passa sobre uma trajet6ria curva. Este deslocamento lateral relativo faz com que o raio de rolamento da roda ex tema, devido à conicidade, seja maior do que o raio da roda interna à curva. Desta forma. como a rotação angular do rodeiro é idêntica para as duas rodas (rodeiro considerado rlgido não havendo, portanto, movimento rotacional diferencial entre as rodas) haverá produção de forças longitudinais diferenciadas entre os pontos de contato de cada roda. Estas forças tenderão a produzir um ângulo de ataque do rodeiro em relação à direçi!.o da via férrea buscando a inscriçl!o do rodeiro na curva.

Este efeito restituidor garante a centralização do rodeiro quando trafegar em trajetória retilfnea pelas irregularidades da via férrea lnduz também a contribuição individual de cada roda na geração das forças laterais (centrlpetas) necessàrias para inscrição de trajet6ria curvilínea (melhorando inclusive a segurança contra o descarrilamento) garantindo a guiagem automática do rodelro em curvas. Esta propriedade, entretanto, resulta num sistema dinâmico com freqUência natural definida, amortecimento modal inversamente proporcional à velocidade e podendo apresentar velocidade critica acima da qual o sistema toma-se instável.

Equacionamento O rodeiro ferroviário pode ser representado pelo sistema mecânico mostrado na Fig. I. O sistema

de referência utilizado está vinculado à estrutura do truque e trafega junto a este a uma velocidade constante V & O rodeiro foi modelado com dois graus (n=2) de liberdade: deslocamento lateral do rode iro em relação à via uy c rotação angular q> ~, na direção z conhecido como ângulo de yaw.

y

-.. -r r----1 - ---1 ' -,: .. ~~r .. ~-

1 '• ·;~ ......

. 10'---'-

. •

.X I

... ---(~· y z

u,~de~

C. ll -do­'1'. - IOp*>do

otoouodo-

Fig. 1 Elemento• do Rodelro Ferrovl6rto

R. S. Barbosa et ai. : Dinâmica do Rodeiro ferroviário 320

A obtenção das equações de movimento deste sistema pode ser feita de várias maneiras. Para modelos simples com poucos graus de liberdade, as equações podem ser escritas manualmente sem muita dificuldade. Entretanto, para modelos mais complexos e extensos, esta atividade torna-se desgastante e passivei de erros. Programas computacionais para Sistemas Multicorpos (MBS) permitem a geração automática das equações de movimento para sistemas complexos e nao lineares com facilidade, rapidez e segurança.

Neste caso as equações são obtidas manualmente, assumindo pequenos deslocamentos e desconsiderando os efeitos inerciais do truque, a partir da aplicação da 2a lei de Newton sobre o rodeiro nas direções dos graus de liberdade, como pode ser observado na Eq. (I).

[m OJ { Üy} cy O { uy} { T 1 + T 2 } { Fy} (I) O e cPz + O cxe~ <Pz. + bo(~xt-;x2) = Tq,

As forças desenvolvidas no contato T xi e T yi entram do lado direito da equação como forçamento externo. Entretanto, devido a mecânica de contato, estas forças são proporcionais à velocidade relativa entre as superflcies de contato, dependente dos graus de liberdade e suas derivadas podendo, portanto, tomarem-se parte integrante do sistema passando para o lado esquerdo da equação geral.

De maneira simplificada, conforme a teoria de mecânico de eontato, pode-se exprimir as forças nas direções longitudinal T xi e lateral T yi como sendo proporcionais às velocidades relativas 9 xi e 9yientre as superflcies de contato rodá/trilho. As constantes de proporcionalidade, kx e kY' que relacionam os micro-escorregamentos entre as superflcies (conhecidos como creep) e as forças são expressas da seguinte forma:

9 9x2 Txl = k 2.! e Tlt2 = k.xy (2)

X VO o

svl 9 T = k~ e T =k....l! (3) yl y v o y2 y v

o

Os valores de micro-escorregamento resultantes das velocidades relativas no contato são função das coordenadas do sistema e suas derivadas, da conicidade da pista de rolamento À. raio nominal de rolamento r 0 , semidistância dos pontos de contato b0 . Como r, = r., + ). u, e rl = r. - A. u, as forças são obtidas pelas expressões a seguir:

T l = -k (~ + <i>zbo) e T 2 = + k (~ + <i>zho) X X rO )/0 X X TO VO

(4)

Tyl = Ty2 = ky(-q>z+ifj (5)

Substituindo estes valores na Eq. (I) e rearranjando os termos de forma matricial, tem-se a seguinte

expr:~J: { Üy } 1 r2~ 0 l { Uy } [ cy -2kyl { u } { F } (6)

lo e ii>1. + Vol O 2~b; <i>z + 2kx'-b/r0

cxe; <P: T: Observa-se que após a inclusão das forças de contato aparece o termo de primeira derivada do

vetar de coordenadas, revelando que o sistema possui amortecimento relacionado com as propriedades de contato e inversamente proporcional à velocidade.

A Equação 6 é valida apenas para trajetórias retillneas e pequenos deslocamentos (angular e lateral). Durante a inscrição de uma curva, a trajetória do rodeiro se altera em função da solicitação externa imposta pela variação da posição da via no plano, combinada com as características de resposta do próprio sistema. Considerando a situaçao de regime da inscrição em uma curva, a força centrlpeta aplicada sobre o rodeiro e o torque proveniente da rotação deste em relação ao centro da curva, correspondente as força externas Fye T q,' e aparecem do lado direito da equação.

321 J. oft~ Sraz. Soe. Mechanícal Sclences • Vol. 18, Oecember 1996

Na situação de regime, durante a inscrição de uma curva de raio R, a solicitação externa sobre o rodeiro é composta de duas parcelas: aceleração centrfpeta expressa por V2/R e o torque produzido pela velocidade angular de regime do rodeiro 't' = V

0/R, conforme Fig. 2.

Como o rodeiro é considerado rfgido, o acréscimo de velocidade relativa para cada roda devido à velocidade angular na curva, é dado por: Sx 1 = V 0 b0 /R e Sx2 = V 0 b0 /R (pois Sx = 'l'b0 ). A força centrípeta e o torque decorrente do contato, normalizadas por V & são introduzidas do lado direito da Eq. 6, resultando na equação de movimento do sistema mecânico em curvas:

[m ol{ üy} 1 r2ky O ]{ Üy} [ cy -2kyj{ F} { mv2 /R} O ej éi>t +V o lO 2kxb~ lPt + 2kxÂ.b

0/r

0 ex e~ T: + 2kx:~ /R

adotando r1

.. { uy <Pt}1

• r matricial, obtém-se:

l (7) { (F Y + mV~ / R) (T lP + 2kxb~ /R)} e resumindo a expressão

[MJ {i'}+ ~o (C) {r}+ [K) {r} { f}

~'---~- ____ R _________ .~

r~r Go - ' ,, ~ : 11< ---1 2~)

/ ;:,· .. ., I ORQUR

'l'lWlSlÇÀO f ' . R,

• -o ~· - ~: -·~~.

. ... s

l G ..!..!. R R0 s0

Fig. 2 Oescriçlo de Curva Cln:ular e de Tl'1n&lçlo

(8)

R. S. Barbosa et ai.: Oinamica do Rodetro Ferrovtãno 322

Solução do Sistema de Equações

A soluçllo do sistema de equações diferenciais de segunda ordem é feita com a redução do ~rau do sistema passando para a representação de estado de espaço com x

1 = { Üx <Pz u q> } , s ua

respectiva derivada x1 = { Üy iPz Úy ~z}1 e u1 .. {f 0 }1 • resultando num SIStema tonÕnuo de

ordem 2n do tipo x = Ax + Bu. Neste caso obtém-se:

{x} = [~o [M)-1 [C] - [M] - 1 [K]l {x} + [ [M]-1] {u} (9)

[ I] [O] [O]

Para sistemas dinâmicos homogêneos com coeficientes constantes a solução temporal continua direta (Gasch e Knothe, 1987), que leva o vetor de estado (x}10 de condições iniciais no tempo t0 , para o tempo t qualquer, é dada por:

(lO)

A forma construtiva da matriz fundamental [ Cl>] para sistemas amortecidos é expressa por seus autovalores arranjados de forma matricial multiplicados pela matriz modal :

-I (Cl>] (1 - t.,} = [ R] [/\] (t - 1.,} [R]

onde a matriz diagonal [i\] é constitulda da exponencial dos autovalores À.i como segue:

Àl (1-t.l o o e

[ /\] o Ã.,(t - 1.) e o

o o e Àla (I - t.)

e a matriz modal [R] constitulda dos 2n autovetores {v} da seguinte forma:

[R) = [ {ri} {r2} ... {rn}] = [À.IVI À.2v2 ... ~n"2n] vi v2 ... v2o

( II)

{12)

( 13)

A solução completa do sistema de equações diferenciais ordinárias não homogéneas pode ser obtida de forma continua no tempo com auxílio da integral de convolução para uma excitação externa {u} qualquer a partir das condições iniciais {xl(to) do s istema por:

{x}<• - •.> = [Cl>}(t - 1.) [x)lt.l + ~[Cl>) (t - t} [B) {u}(tld-r o

(14)

A representação em tempo discreto da Eq. 14, assumindo o interesse apenas nos instantes de tempo Tk = k T de período T = t- t

0 igualmente espaçado com excitação constante, terá a forma:

( 15)

323 J . of the Braz. Soe. Medlanlcat Sciences- Vol. 18. Decembef 1996

resultando na expressão final completa na forma matricial utilizada para o cálculo do vetor de estado no instante ltc+t a partir do estado do instante tk com a matriz dinâmica do sistema (A], matriz de transição [<!>] e matriz [B] de combinação de aplicação do forçamento externo { u } (conforme proposto em Ogata, 1993):

{16)

Análise das Propriedades do Sistema

O cálculo das autopropriedades da matriz dinâmica (A] do sistema, de ordem 2n, permite obter as caracterfsticas naturais e modos de vibrar do sistema parametrizado em função da velocidade V o.

Observa-se na Fig. 3a os autovalores do primei ro modo, pares complexos conjugados, apresentados no plano complexo para diferentes valores de velocidade. Este modo corresponde ao movimento de passeio lateral do rodeiro conhecido como Lacd ou Hunting (ver Fig. 4a).

A freqUência natural deste modo é linearmente proporcional à velocidade, conforme apresentado na Pig 5a, o que corresponde a um comprimento de onda de movimento lateral ôL aproximadamente constante (Fig. 6b • 11 ,5 metros para este caso). O fator de amortecimento ~ possui valores em tomo de 0.24 atenuando-se com o aumento da velocidade. A 60 m/s o autovalor apresenta parte real positiva (fator de amonecimento menor que zero) caracterizando a velocidade critica acima da qual o sistema toma-se instável, conforme apresentado nas Figs. 3a e 5b.

lnspecionando os autovetores deste modo, observa-se que possuem valores complexos, o que representa movimentos relativos com ângulo de fase. O attaso do ângulo de yaw (cp1)em relação ao movimento lateral na direção lateral y ( uy) é de aproximadamente n / 2 (Tabela I - I 04,19° e Fig. 6a ). Portanto, quando o movimento lateral é mâximo, o ângulo de yaw é próximo de zero (ver Fig. 4).

lO ·c: 1! ·c;,

§ Cl ~

~

a. Autovalores do Primeiro Modo (Plano Complexo) 4or

20 1

o-'

-2l -40

,.+I + _.,.. f~ .f - -I t -f; -i

+ i-

_.___ ' ---'- ...J __ __.~

-3.5 -3 -2.5 -2 -1 .5 -1 -0.5 o 0.5

50

Parte Real b. Autovalores do Segundo Modo (Plano Complexo)

~~--------~ -~----. I I I I I I I I t

I I

' I I

I I I '

o L l J . . I

· :Veloc~ades (mls~~) : ; ; : ; ~· : ' o I I jq' I I I I I I

! • I I I I I I I I I ' I

:- ~ -: ~ --~ ~~ ~lW~J&- --- : . ---1 I I I I I I 1 I I f

I +++. '+-+~4

I 1 I I ! : ~, I o I I

o I

-SOL.L......_L-1.__,__. __ _,

-10.

Fig. l AutoVIIOrH do Primeiro (x,+) e Segundo Modoe (x, +)

-1tt Parte Real

:. S. Barbosa et al.: DiMmlca do Rodeiro Ferroviário

L Descrição da Trajetória do Rodeiro

c:s cr-I =-~=====-:t~I~=~--~:L~;~==~l~\~~·~·

Via Principal

Comprimento de Onda (óJ -------

.Movimento de Passeio Lateral (Lacet)

b. Descriçlo do Aparelho de Mudança de Via

I ~ -~- --]..~ ... ,..,.,. I ..

. - I - -L . 1'···- • '·~ .... - . - - . - _ .,.. 1 I 1 ··-~......_

~~-1-- I ~ ---+-=:-t- / -

I ~-~-r-..... ~/ Trecho Reto ! Transiçao / "" ~

/ í '\ / I

CUrva Circular (Raio R)

Fig. 4 OeKrtçlo di Tr1jetóna do Rodelru • do Apar~lho de MudenÇI de VIa

324

325

~ Q)

E

~ o E < ê5 êií u..

Ê l1l 'O c o .9 c Q)

E ·c a. E o (.)

J . o1 the Braz. SOe. Mechanlcal Sciences- Vol. 18, December 1996

6! ,- ;-1

:~1- --: : -": ~;_:_~_I~~~;~:_;~-~:-----:1 ~ o o I

ol ------ . ~ o 10 20 30 40 50 60

Velocidade (mts)

0.3~---~----

1 b I--~---

0.2r --- --- · ~ - - - ~---=~-=-.:..-=--=-- - --~ - --- -- -:---- ---I ~ , -------· ' ' + ' ~--- ' O.fr _ _ _ _ ~ _____ __ ~ ___ ___ -·- __ ____ + __ ~- - · - __ __ _

I • I I --,

o - - - - - - ~ - - - - - - - ~ - - - - : - - - .. - - - ~ - - - - - - -: -~=-=--:-:.--0.1L_ ___ .__ _ _ _ ,_ ___ ,_ _ _

o

11.8,

11.7

11 .6

11.5

11.4

11.3 o

10 20 30 40 50 Velocidade (m/s)

Fig. 5 Proprled1dn do Prtmelro Modo: Frequlncl• JQtun~le F.tor de Amortecimento

' ' -·~~-· --- .. --. !... -- - --- - '--- - -- ---- - -- -- ,_- - -~-.. __ I I I

-- I -------.---.. '=--~~ . - - .. - - - - . - - - -· - - - ,_ - - -

b

' O ........ ....___ I

-- c;,~icld~d~ ~; R~a~ ~;1õ -~-~--- -:- -- -- -- ~ ---- ---··---_.-- ~~ I . - - --- ' - - - :~.::----'-------

--'----10 20 30 40 50

Velocidade (mls)

Fig. S C.111ctenatlcu do Prtmelro Modo: .lngulo de Fue do Primeiro Modo e Comprimento de Onda do Movimento Latel'lll

60

60

R. S. Bartosa et ai.: Dinâmica do Rodeiro Ferroviário

Tabela 1 Autovetorea dos Modos para Velocidade de 40 m/s

Coordenadas

Deslocamento Lateral (Uy)

Ângulo de yam ( q> ) z Ângulo de Fase

Autovetor 1° Modo

-0,0174 + I 0,0353

+ 0,0212 + i 0,0045

104,1go

Autovetor 2'1 Modo

-0,0045 + i 0,0002

-0,0003 + i 0,0027

93,42°

326

O segundo modo possui um par de autovalores reais distintos (sistema sobreamortecido) até 20 m/s conforme apresentado na Fig. 3b. A partir desta velocidade os autovalores passam a ter parte imaginária diferente de zero, formando um par complexo conjugado subamortecido (Fig. 3b). Em baixas velocidades, o módulo dos autovalores reais (correspondente a constantes de tempo) são inversamente proporcionais à velocidade (Fig. 7a). A partir de 20 m/s a freqüência natural amortecida cresce rapidamente para valores idênticos ao primeiro modo (Fig. 7b). O fator de amortecimento deste modo é elevado chegando a 0,97 para velocidade de 60 m/s (Fig. 8a). O ângulo de fase observado na Fig. 8b, cresce rapidamente com a velocidade, sendo 93,42° para velocidade de 40 m/s (Tabela 1)

... !. ..

~ ~ " ~ " o ..

> ::i 2

" <

-;;

~ !. .:? • • ~ z ...

~

"' o • ii:. E <

a 800 ~ : : . I

:::f \ :: : ~ : -: : -: -- ---; -- . -- ' : . ----. --j 200 r -~-- ----.-- _ _ ___ __ I

O- - I I

o , o 20 30 40 50 60 Veloci<la<le (mi$)

61 ' ' ' .,-- ----~- --- --~- ---- -~------_: :-------+-----~~~ 2r------ ~ ------ ~. --7-- --2-~~~~~-;----- - -~ -----. ~

I / I ' 0 ' ' ' I _.

o 1 o 20 30 40 50 60 Veloci<lade (m /a)

Fig. 7 Proptf.ctades do Segundo Modo: Módulo doa Autovalores FreqUência Natural Amortecida

i l : :: ------; -----l -: : -: : : : : : :-T~--:-:<: 0 1 O ~---. I :--J IL E o. 9 7 - - - - • :. - • - - - - i ------ ~ --- --- -: . -- . --. : '

.... (!)

~ c .. .. .. ., IL

<(

0.96 o

_J ________ J_ ______ ~---------

10 20 30 40 50 60 Velocidade (m/s)

-2:l _ --- --~ ------~\_ -_--- -~- . -~- ~ - -_b_ -~ ·40[ : : ·~ ___ : __ ____ :_ - __ : __ ___ _

1 1 1 i i

-60 r - - - - - - r - - - - - · ~ - · - - : ~ - - - - - - -:- - J I ' I ....... ...t.. I '

·80 I - r - - - - - - -i .::-.:-_ __ - "' I • - - - - -•- - - - -

' ~ ----~ ------· -----100~----------------·L-------

o 1 o 20 30 40 50 60 Velocidade (m/s)

Fig. 8 C41racter1atlcaa do Segundo MOdo: Fator de Amortecimento do Segundo Modo Ângulo de Faae entra oa Graus de Uberdade

327 J . ofthe Braz. Soe. Medlanlcal Sdences • Vol 18, December 1996

Avaliação da Sensibilidade A Equação 7 foi parametrizada em V o onde foram analisadas as autopropriedades. Entretanto

outros parâmetros possuem variações, inclusive não Lineares, do sistema mecânico real. Observando o perfil da roda constata-se que a conicidade varia bastante próximo da região do friso. Adicionalmente a rigidez das forças de contato roda/trilho é influenciada, entre outros fatores, pelas proporções da elipse de contato (alb). Toma-se conveniente o conhecimento das tendencias das autopropriedades devido a variação destes valores.

Para um determinada velocidade constata-se, por inspeção da matriz dinâmica do sistema (A], que para o primeiro modo, devido ao aumento da conicidade da roda ),. ou aumento da rigidez do contato (~e ky):

FreqOência natural deste modo aumenta; Velocidade critica diminui; Fator de amortecimento diminui e Comprimento de onda diminui.

O comprimento de onda não se altera com a variação da rigidez. Entretanto, o módulo da raiz do segundo modo aumentou bastante especialmente em baixas velocidades (diminuiu a constante de tempo). A conicidade afeta pouco as características do segundo modo.

Simulação A partir do modelo construido é possível simular, utilizando algoritmo apropriado (Barbosa,

1993), o comportamento do rodeiro na inscrição de uma curva ou em um aparelho de mudanças de via (AMV) com variaçllo da trajetória da via. Esta é a situação onde se encontra maior adversidade para inscrição do rodeiro e local de elevada incidência de acidentes.

Um AMV é constituldo basicamente por três elementos primários de geometria: trechos de reta, curva de transiçllo e curva circular. Ao inscrever o AMV o rodeiro se depara inicialmente com uma mudança angular repentina da direçllo da via (ângulo de kiclc). Este trecho corresponde à agulha, que nos modelos mais simples de AMV é um segmento reto de trilho que desloca a direção da via lateralmente iniciando o desvio. Após este trecho inicia-se a curva de transição com variação linear de curvatura (clotóide- ver Fig. 2), variando de zero a IIR e, finalmente, o trecho em curva circular de raio R. A Figura 4b apresenta de forma esquemática os elementos do AMV.

Uma aproximaçllo simplificada da entrada da agulha do AMV (trecho reto com ângulo ex de kick) pode ser feita por um pequeno trecho de grande curvatura, resultando num pulso de curvatura de valor R ~ V o ~tI a (onde ~t =intervalo de cálculo).

O resultado da simulação do comportamento do rodeiro entrando em um AMV pode ser visualizado na Fig. 9. São apresentados o deslocamento lateraJ do rodciro em relação ao centro da via (Y) e o ângulo de ataque do rode iro em relação à direção da via (ângulo de yaw).

201 I lnsa1cao 11m AMo/

,J f • - ~ • • • • : • ~~ • • - I ~~ ~~~~ •- • -

I\ : : I \ I I '

I I I · f \0, lrkioet.w 1rm

I ' L \ I _.:;,;;.-x--~ _ __.

~ 101 ~ \ - : • • • r/ -- T : =~, ~~ ---

i J ~\ \ ~ -(\}-I ~ --- -~- ; ~ J\ \ , '\ , ~2m I I I :(\ \ : : > o t - \' ~ 1- \ f- - \ ~"-..::..>--

l I ! I ' V · lleioc:lcledl: 20 ,.. I v· , . .

-51 · ~~ ·· ··- ·- · ----·101'-----'-

0 10 20 30 4() 50 eo ~(m)

Fig. fi ReauiUidos da Slmul~o da Dlnlmlca do Rodalro

R. S. Barbosa et ai.: Dinâmica do Rodeiro Ferroviário 328

Observa-se durante o inlcio da inscrição do AMV (localizado à 2 metros da posição inicial) a mudança repentina do ângulo de ataque devido a agulha e do deslocamento lateral. Em seguida o rodeiro reverte o deslocamento lateral buscando a centralização, atingindo deslocamento lateral de 14,8 mm no trecho reto da agulha a 4,5 metros do inicio da simulação.

A partir de 17 metros inicia-se a curva de transição de 2 metros em seguida a curva circular. O valor máximo de 15,5 mm de deslocamento lateral ocorre na posição de 24,7 metros, estabilizando posteriormente a oscilação em tomo do valor de equillbrio de 10,95 mm de deslocamento lateral na curva circular de 300 metros de raio.

Conclusões

Foram apresentados os elementos mecânicos para elaboração do modelo linear do rodeiro ferroviário. As forças no contato roda/trilho, fundamentais para a representatividade do sistema, foram modeladas e incluldas nas equações. As autopropriedades do sistema parametrizado em função da velocidade V 0 foram calculadas e analisadas.

O primeiro modo de movimento com freqUência baixa corresponde ao passeio lateral do rodeiro em relação à via. Este modo possui comprimento de onda de movimento no espaço aproximadamente constante para a faixa de velocidades investigada, sendo fortemente influenciado pela conicidade À. da pista de rolamento da roda. O amortecimento modal apresenta valores decrescentes em função da velocidade, tornando-se menor que zero para velocidades acima de 60 rnls, o que corresponde a um sistema instável, definindo-se uma velocidade critica.

O segundo modo, sobreamortecido, apresenta raizes reais distintas com módulo inversamente proporcional à velocidade e fortemente dependente da rigidez k.,. e ky do contato. A partir de 20 rnls as raízes tomam-se complexas conjugadas (subamortecido) e com valor de freqüência natural amortecida que se aproxima do primeiro modo em alta velocidade (Figs. 5a e 7b).

A simulação da inscrição do rodeiro em via com trajetória variável, permite observar o comportamento dinâmico do rodeiro. O cálculo dos movimentos permite identificar os locais de maior deslocamento lateral e ângulos máximos. Apesar de simplificado, este modelo permite visualizar o comportamento rodeiro na sua interação com a via para percurso com trajetórias variáveis como é o caso de AMV's. Desta forma é posslvel investigar o comportamento dinâmico do rodeiro em vias curvas gerando subsidio para analise, concepção e projeto de AMV's.

A análise do sistema linear permite identificar o substrato das caracteristicas do sistema dinâmico como preparação para implementação do modelo tridimensional completo não linear utilizando a técnica e programas de Sistemas Multicorpos (MBS).

-- Dados do Sistema

Conicidade da Roda Raio Nominal da Roda (36") Semilarg. Rodeiro (bitola via) Semidistãncia da Susp. Primária Massa do Rodeiro Momento de Inércia do Rodeiro Rigidez lat. Susp. Primáría Rigidez long. Contato Rigidez Lat. Contato Rigidez Tore. Susp Prim. Rigidez Long. Susp. Primaria

Referências

À = 0,10 T0 = 0.4572 m b0 = 0,7175 m 8 0 = 0,61 m m = 1751 kg e =soo kg m2 cy = 1,00 x 102 N!m kx ::6,20 X 1o& N ky = 6,50 x 1 o& N kzz = 1,70 x 1Q8 N mlrad. ex = 4,5687 x 10e N/m

Barbosa, R. S., Setembro, 1993, "Estudo da Dinâmica Longitudinal do Trem", Dissertação de Mestrado, Universidade-Estadual de Campinas (Unicamp), Campinas, SR!> Paulo, Brasil.

Costa, A., 1994, u Application of Multibody Systems (MBS) Techniques 10 Vehicle Modelling" S.A.E. São Paulo.

329 J . af lhe Braz. Soe. Mechanat S<:iences - Vol. 111, December 1 ~

Costa, A~ Barbosa R. S., Weber, H. I. e Pascal J. P., Jwlho 1994, uum Programa Modular de Sistemas Multicorpos (MBS) para a Simulaçlo de Velculos e ComposiçOes Ferroviárias", Seminário Franco Brasileiro de Transportes Pôblicos, Slo Paulo.

Gasch, R. e Knotbe, K. , 1987, wStrukturdynamik Diskrete Systeme", Springer-Verlag, Vol. I. Ogata, K., 1993, "Discrcte-Time Systems", Prentice-Hall, 2" Ed.

RBCM • J. ot the Braz. Soe. Mechenlcel Sclences Vol. XVIII· No.4 • 1996- pp. 330..337

ISSN 0100-7386 Printed ln Brazil

Cogeneration System Design Optimization José Antônio Perrella Balestleri UNESP - Universidade Estadual Paulista Faculdade de Engenharia de Guaratínguetá Departamento de Energia 12500.000 Guaratlnguetá. SP Brasil

Paulo de Barros Correia UNICAMP - Universidade Estadual de Campinas Faculdade de Engenharia Mecênica Departamento de Energia 13083-970 Campinas, SP Brasil

Abstract Cogeneration system design deals with severa! parameters in lhe synlhesis phase. where not only a lhermal cycle must be indicated but lhe general arrangement, type, capacity and number of machines oeed to be defined. This problem is not trivial because many parameters are considered as goals in lhe project. An optimization technique th8t considers costs and revenues, reli8bility, poUutant emissions and exergetic efficiency as goals to be reached in lhe synthesis phase of 8 cogeneration system design process is presented. A discussion of appropri8ted values and lhe results for 8 plllp and paper plant integration to 8 cogeneration system are shown in order to illustrate lhe proposed melhodology. Keywords: Cogeneration, Optimization. Multiple Programming, Design, Efficient Solutions.

lntroduction

Desígn can be defíned as any activity developed in a synthesis phase, where a new or a modified concept is created to satísfy certain demands. and in an analysis phase, that involves testing the concepts proposed in the previous step in its capacity to satisl)t such demands.

Cogeneration system planning and design can be developed according to analytical or optimization models that will permit the definition ofthe equipment and operational procedures.

Some important questions must be answered in the cogeneration system planning:

How is the thennal cycle to be adopted? Which equipment must be chosen?

How many units of each equipment is indicated? How is the individual capacity of each equipment?

Which enthalpy levei drop must be explored?

Both analytical and optimization models can be used to generate differcnt technical solutions to a specific case-study; ali the proposed solutions can be grouped on a rectangle, as proposed by Bohem, 1987, and represented in Fig. I.

Nomenclatura EC02 = emission factor of FIX = equipment capital cost. pu fuel costs. US$/kg fuel

co2. kg co2;k.g tuel 106 US$ rhs steam demand in the ES02 = emission factor of s feasible region of a process, 106 kg/year

so2• kg so21kg tuer problem X exergetic loss factor, ENOx = emission factor of m = flow rate. 1 OS kg/year xi = 1- T0(ós/óh)

NOx, kg Noxlkg fuel P(ij = power generated, WNV F = capital recovery factor, nK number of ares in a óh = enthalpy drop, kJ/kg

year·1 netwof1( ST steam turbina M = operating/malntenance pe = earning for selling GT = gas turbine

costs electricity, US$/kg HRSG = heat recovery steam steam generator

Presenteei at the 1995 ASME lntemational Mechanícal Engineering Congress and Exposition- San Francisco. CA. November 12-17. 1995. ASME permlssion to reprint the paper is acknow1edged. Associate Technical Editor: Carlos Alberto Carraaco Altemani.

331 J ol the Btaz. Soe. Mechanical Scien<:es- Vol. 18, Oecember 1996

Ali Proposed Solutions

Efficient SOl utlons Feasíble SOiutions

Fig. 1 SOiutlon• Spac• (Bo.hm, 1987)

Feasible solutions are those that satisfy ali constrains of a problem, considering thut it can be fonnulated as a set of cri teria to be reached subject to some constraints.

When an optimization problem bas only one objective to be reached, the solution obtained will be an optimum point; if the model is used to generate maximum andlor minimum solutions according to multiple criteria, each one will be an optimum solution for the respective criterion and they can be grouped on an efficient set of solutions ,each one related to this individual goa1.

From an economic point of view, it is intended to synthesize a cogeneration system that reaches minimum operating and capital costs with maximum production of electrical surplus to be sent to the grid; point A can represent this speciaJ solution, that is an optimum point.

Cogeneration plant reliability can be another important criterion to be considered; these systems must be planned to operate on a continuous condition or at design point for many hours during the year (more than 6500 hours/year) and must have redundant equipment to warrant the continuity of electrical supply even in the case of a fault of some generators.

Environmental regulations are another important criterion to be considered in a cogeneration system planning. C02, S02 and NOx emission leveis are becoming more and more restrictive, conducting the optimization problem to severa! solutions, perhaps different to the previous ones.

The fact of choosing one single criterion in the design of a cogeneration system will imply in a project that reaches the maximurn (or minimum) for that goal but with unsatisfactory (or perhaps not acceptable) results to ali the others. A trade-off methodology must be defmed for helping in the decision-making process.

Multiobjective Condition of Design

The planning of a cogeneration system can be developed by means of a simulation process (varying thermodynamic conditions of a proposed plant)) with satisfactory results. However, in the case of dcsigning new plants, this too I can only be used after the choosing process of a cycle and the definition of a certain scheme.

Optimization models can be used successfully in the design and operation of cogeneration plants. OperationaJ analysis is usuaJiy done by fonnulating linear (BaJestieri and Correia, 1994).

Cogencration system operation optimization models generally consider as objective functlon an economical expression of net benefit that maximizes the revenues from electriclty and minimizes investmcnt and operational costs (Ehmke, 1990; Nath et ai .. 1989: Pungen and Macgregor, 1989).

J. A P. Ba~atieri et ai.: Cogeneration S~tem Oe~~gn Optimlzation 332

Suppose that point B of Fig. 2 represents the optimum v alue for net benefit in the feasible region S: this solution maximizes the vector b that represents this objective function .

Considering the sarne constraints but using now as objective function an expression related to the cogeneration system electrical reliability, it is expected that a new solution, represented by R, wiJJ be generated and will be certainly different to the other; vector r indicates the preferential direction of reliability function.

x2

Fig. 2 F ... lble Reglon ln the Dec .. lon Sp.~~ce

Schemes (R) and (B) of Fig. 3 are the corresponding scbemes (physical solution) of optimum v alues rea.cbed by net benefit vector b and reliability vector r.

------------------~-----~--(R) Maximum Reliability Scheme {8) Maximum Nct Benefit Scheme

Fig. 3 Optlmum Propoeed Schemn

ln the feasible region S it can be seen that only solutions placed on RB segment have a relationship of gain and loss between reliability and net benefit: coming through RB segment from B to R reliability is decreasing but net benefit is increasing and vice-versa.

Ali solutions in RB segmentare named efficient, non-dominated or Pareto's optimum and, from an efficient point, it is not possible to move feasibly to another in order to increase one of the objectives without necessarily decreasing at least one of lhe others.

ln cogeneration system planning several objectives are considered subject to the sarne constraints and it is difficult to represent such variables as a surface, so a mathematical too! must be proposed to help decision-makers define the best scheme to a case-study by adequately weighting the goals according to his preference structure.

A multiobjective interactive model that helps to define thermal cycles in this first-step and the size, capacity and type of each equipment in its second-step is presented in the next section.

333 J ofthe Braz. Soe. Mechenical Scienoes- Vol. 18, Deoember 1996

Multiobjective Approach A multiobjective interactive programming model that coosiders in each iteration the decision­

maker's preference structure for a set of objective functions was indicated to geoerate efficient solutions used in the evaluation of a specific case-study.

lnteractive procedures generally conduct an exploration over the region offeasible altematives for a satisfactory efficient solution. There are altemative steps where the computer leads to a solution in one phase, so the decision-maker inputs information and intervenes in the proposed solution in the other phase. This interaction between man and model enables the decision-maker to learn more about his problem and to have more insight to decide.

A two-step multiobjective model (TSMM) was formulated to permit a faster and more efficient search for solutions; initial phase defines the thennal cycle to be adopted, coosidering steam, gas and combined cycle as possible options. This step reduces the original problem by discarding the worst options.

Final phase uses this infonnation to define, after proposing several solutions, the cbosen scheme, type, capacity and quantity of each equipment as well as general tecbniques and economic data relevant to the cogeneration system planning - cost, reliability, pollutant emission leveis, installed power, etc. (Balestieri, 1994).

A mathematical structure of a generaliz.ed network was indicated to solve the linear programming problem because of its advantage over the conventional Simplex method - computational efforts are reduced by about 50 times (Giover et ai., 1978). Another important reason is that the network. defined in terms ofnodes andares, can be well represented in those latter mathematical entities. Equipment's energy conversion efficiencies can be associated to the gain variable of each are and this can mathematically express the physical phenomenon of converting certain flow rate of exbaust gases from a gas turbine in a fewer flow rate of steam in the heat recovery boiler.

The GeotTrion, Dyer and Feinberg procedure (Steuer, 1986) was chosen to be part ofthe search algorithm since it is compatible with network structure and don't impose new constraints to the original problem; oplimal solution for each objective is provided according to the sarne set of constraint equations and net benefit optimal solution is taken to tbe position ofthe problem's initial solution.

The model presents a new solution an four intennediate solutions are generated by means of pairwise comparison questions. Decision-maker is then invited to express his preference by choosing a solution and deciding about re-start or finish iteration (Fig. 4). ~ -- -- -- -- -- -- -- I ·I ......,._ l-' I - j

:l~--11

...,.__ .

. -.... -· -(...,......a&M)

I[F , ~-

1

L_ lnitial Phase

.,.. L----....J

L-----' r ---:- --=... _j_j _. rr......,_,.._,

I.- ....... ·-(~lot) ....... .---I . 802iNoldC02

L----

Fig. 4 TSMM Algortthm

Final Phase

J. A. P. Balestieri et ai. : Cogeneratlon System Oesign Optlmlzatlon 334

Mathematical Modellng

The following mathematical functions describe the set of objective functions and constraints for the cogeneration system planning:

ok

Min 2: (x;, ôH;) mi i • 1

2, 46

Min L EC02 ;m; ã = 1 •

2, 46

Min .L ES02 i mi i_, .

2, 46

Mrn L ENOximi i • I

ok

Max2:m; i - J

ok 5

Max L (pe;m;) - L (pu;m;) -L (FIX/F)- M i= I i= 3

subject to: mí e S

(I)

(2)

(3)

(4)

(5)

(6)

Objective function (1) aims to lead the solution to the scheme that has the machines with the least thermodynamic losses (measurcd according to the exergy loss factor 'Xt) and, therefore., it is an energy conservation criterion.

Functions (2), (3) and (4) conduct the solution to a C02, S02 and NOx less poUutant schemes. respectively, reflecting the environmental question of the 1990's. Objective function (5) is related to the plant reliability and is stated considering that a scheme will be as highly reliable as the most electric generation machines are available in the cogeneration system (Sullivan, 1977). It is still an important way to calculate the LOLP (Loss ofLoad Probability) according to a recursive equation.

Objective function (6) express the net benefit obtained with the congeneration system operation; it bas frrst term that aims ofmaximiring the revenues from electricity according to a certain buyback rate proposed to the model within the minimization of fuel consumption, capital and operational & maintenance costs on the other terms.

Problem constraints rnj e S were stated considering that mass must be conserved in ali nodes and that medium and Jow pressure process must be satisfied in its steam requirements ali the time (thermal parity). ln this study we considered that another altemative could be formulated, such as the electric parity.

A General Design Module (GDM) is represented in Fig. 5 for the final phase ofTSMM; it contains aU possible altematives to be chosen as optimum schemes to every objective function and yet to every efficient solution proposed by TSMM (number 3 to 8 are convention steam generators, 9 to 13 are gas turbine coupled to heat recovery steam generators and 14 to 22 are back-pressure or condensing steam turbines, with or without extraction).

335 J . ofthe Braz. Soe. Mechanical Sciances- Vol. 18, December 1996

Fig. 15 General Dtalgn Module

Case-Study A pulp and paper industry has an annual demand ofmedium pressure steam (1.3 MPai275°C) of

1263x l 06 kg/year and of low pressure stearn (0.3 Mpa/1700C) of 1230x 106 kg/year supplied by a generation system that has pressure stearn generators (4.0 MPa/3500C) buming fuel oil anda by a condensing turbine with two extractions. This unit is at the moment generating electricity with a little deficit that is supplied by the grid.

An analysis of process conditions showed that it is possible to condense 2500x I 06 kg/year ; considering that it is necessary to have an installed capacity of36 MW to satisfY the electrical needs of the process. TSMM is started by defming the optimal solution for each objective function individually. Table I presents information used to simulate TSMM

Table 1 General Assumptions for the Case-Study (Batestleri, 1984)

Emlsslon Leveis (kg pollutantlkg fuel) Fuel co2 S02 NOx

o ii 3.156 0.018 0.00728

gas 2.683 0.001 0.00681

biomass 1.380 0.068 0.00010

Annualized Capital Costa of Energy Generatlon Systems

Equlpment Flxed Coat Variable Cost (10' US$/vear) (US$/kg)

GT-HRSG 0.780 0.0002364

HP-boiler 0.232 0.0005857

Backpressure ST 0.185 0.0004186

Condensing ST 0.670 0.0034740

Based on these results. severa! alternative solutions can be proposed in the TSMM initial phase; in this step it was concluded that whichever altemative routes were possible to conduct mass flow except the expansion from medi um and low pressure to condensing unit.

TSMM final phase was conducted according to these conditions for two different values of buyback rate. Figure 6 presents the scheme obtained via simulation by the TSMM model for a 3 5 US$/MWh buyback rate (table 2 presents technical and economical general data related to this scheme). The analysis of such scheme reveals that it is recommended, in this case, the medium pressure stearn generation by black liquor burning and the use of steam turbines from high to medium pressure levei and from medi um pressure to condensing line.

ln the second simulation exercise, a buyback rate of 95 US$/MWh indicated a new scheme with a high incentive to the power generation ln a combined cycle, no more black liquor burning is recommended, so the disposal ofthis byproduct in the agricultura! sector must be considered, as well as its use as an altemative fucl in the high pressure steam generations or in the gas turbine cycle by gasification process (Fig. 7 and Table 3).

J . A. P. Balestíefi et ai.: Cogenetation System Oesign Optimization 336

ln both analysis the decision-maker was expected to generate a solution that could satisfy steam needs of medium and low pressure process with a Loss of Load Probability of around 30% (that is equivalent to about 6000 hours/year of continuous operation, that is a reasonable value for this technology, according to Stambler, 1989) and the lowest payback as possible.

CONOENS•NG

35 USSIM'/Vh Fig. 6 Seheme 1

Table 2 Data of Scheme 1

lnstalled capacity = 162.39 [MW)

Electric demand = 36.00 [MW)

Electric surplus = 126.39 [MW)

LOLP = 2945.27 (hlyear) = 33.62%

Power generated in each turbine:

P(18] = 5.52 (MWJ

P(32] = 5.32 (MWJ

P{33] = 3.56 (MW]

P[35] = 8.56 (MW)

P[38) = 93.59 [MWJ

P[39] = 1.82 [MWJ

P[40] =44.04 [MW)

co2 emíssion = 1195.33 [kgC021MWh]

so2 emission = 0.57 [kgSO;tMWh]

NOx emíssion = 2.29 {kgNOx/MWh]

Capital costs = 141.70 [106 US$]

Payback = 1.55 (years]

Net benefit = 91.42 (106 US$/year]

CONOENSING

MUSSIMWh Fig. 7 Scheme 2

337 J. of the Braz. Soe. Mechanical Sciences - Vol. 18, Decembef 1996

Table 3 Data of Scheme 2

lnstalled capacity c 166.65 [MW)

Electric demand = 36.00 [MWJ

Electric surplus = 130.65 (MWJ LOLP = 2761.00 (h/year) = 31.52%

Power generated in each turbine:

P[18) = 18.80 !MWJ P[32] • 10.22 [MW]

P(38] = 137.63 [MWJ

co2 emlssion = 561.05 (kgC021MWh]

502 emlssion =0.65 (kgSO:IMWhJ

NOx emission = 2.61 [kgNOx/MVVh)

Capital costs =131.79 [1o& US$]

Payback "'0.60 (years)

Net benefit = 246.73 [1o& US$/year)

Conclusions

TSMM has been special1y developed to help designers and deci.sion-makers in the pre-feasibility and planning of a cogeneration system; its most relevant contribution is the wide range of solutions generated in a short time intcrval and the possibility of interacting witb lhe model at any time to alter the preference structure designation of objective functions.

The use of generalized networks affords low computation times and a high accuracy of results. A multiobjective interactive structure teads to a conversational tool for defining severa! parameters of a cogeneration system design.

Some improvements in the auxilíary procedures must be done in lhe proposed model;; the cost functions developed to the thermal and electrical equipment are not well suiting the real values and must be rcviewed. This fact is responsible for the high specific price (868 and 790 USS/l<W, respeclively, for lhe proposed schemes in Fig. 6 and 7), that is expected to be around 650 US$/l<W for units higher than 150 MW.

References

Ahuja. R.K. et aL 1993, "Nc:twork flows", Pn:nticc HalL N. Jersey (USA). Ralestierí, J.A.P., 1994, "Pianning Cogeneration Systems: a Multiobjc:ctivc Approach" (in Portuguese), Ph.D.

Thesis, SIAte University ofCampinas (UNICAMP), Campinas, SP, Brazil. BaJestieri, J.A.P. and Correia, P.B., 1994. "Pianning Cogeneration Systems: a Multiobjective Expansion Analysis'',

Florence World Energy Research Symposium (FLOWERS), Proceedings, Flon:nce, Jtaly, pp, 1025- 1032. Boehm, R.F .• 1987, "Desígn Analysis ofTh.ennal Systems", John Wiley, N.Y, USA. Correia, P.B., 1989, "A Sectorial Model for Encrgy Supply Optimizatioo" (in Portuguese). Ph. O. Thesis, State

University ofCampinas (UNICAMP), Campinas, SP, Brazil. Ehmke,II.J., 1990, "Size Optimization for Cogeneration Units". Energy, Vol. 15, pp. 35-44. Glovcr, F. ct aJ ., 1978, "Generalized Network: a Fundamental Computer-based Planning Tool". Manag. Sei., Vol.

24, pp. 1209-20. Nath, R. ct ai., 1986, ''Joint Optimiution of Process Uoits and Utility System". Chem. Eng. Progress. Vol. 82, pp.

31-38. Puttgen. H.B., MacGregor, P.R., 1989. "Optimum Scbeduling Procc:dun: for Cogeoeration Small Power Producing

facilities". IEEE Trans. PAS, Vol. 4, pp. 957-963. Steuer, Rf., 1996, "Multi pie Cri teria Optimization", John Wiley, Siogapore. Sullivan, R.L. 1977, "Power System Planning", McGraw·Hill, N.Y., USA. Stamblcr, 1. , 1989, ''EPRI Seek:s Dramatic Gains in Maintenability and Reliability." Gas Turb. World, Vol. 19, pp.

15-20

RBCM • J. of the Bru. Soe. Mechank:al SclenCH Vo/. XVIII· No.4 • 199tl· pp. 338-354

ISSN 0100.7386 Printed in Brazíl

Recent Advances and Application of Turbulence Depth-lntegrated Two-Equation Closure for Modeling Contaminant Discharges llren Yu Antonio Marozzl Rlgheto Universidade de S4o Paulo Escola de Engenharia de S4o Carlos Departamento de Engenharia Hidráulica e Sanitána C .P. 359 13560.970 S4o Carlos, SP Bra.sil

Abstract

This paper presents a brief review and up·to-dale development of lhe depth-integrated turbulence two-equation elosure, to be speeific, lhe deplh·intcgrated K- tmodel and R- w model. A numerical simulation of lhe sidc discharge ofwaste heat into natural waters has been taken as an example oflhc practical application ofLhe K- w modcl. lo addition, lhe authors also states the existing problems and further developments oflhe presen1

deplh-integrated turbulence two-equation models. The main objective ofthe paper is Lo integrate thc researcb on lhe improvement and development of lhe present depth-integrated turbulence models witb the practice of industrial, environmental and sanitary engineering. Keywords: Turbulence, Depth-lntegrated Two-Equation Closure, Contaminant Discharge Modeling.

lntroduction

Almost ali tlows in natural watcrs are turbulent. Dealing with turbulent tlow and contaminant transport in natural waters is indispensable for scientists and engineers. For example, turbulent transport associated with the discharge of pollutant and (or) waste heat into rivers and coastal regions is of great concem to civil and coastal cngineers, and also to environmentalists because of lheir serious dam~ging efTects on the ecology and ambient. It is very important to be able to correctly simulate and predict lhe corresponding changes of currents, temperature and concentrations of certain contaminants, as a consequence of lhe anthropic activity in the area.

Although the significance oftimely simulating turbulent Oow and contaminant transport wilhin the aceuracy of engineering is clear, lhe numerica1 simulalion and prediction for natural waters is still at an unsophisticated levei. This is mainly dueto the inherent complexity ofthe problem. Any successful numerical computation and simulation for modeling tlow and transport processes in natural waters depends sensitively upon correctly understanding the mechanics of turbulence and setting up an applicable higher-order turbulence closure model. Up to now two-equation closure models wbich are characterized by comparative rigorous theoretical fundamentais with reasonable physical background remain the only practicable approach to solve engineering turbulence problems. ln arder to obtain lhe necessary information of tlow field with less labor in physica1 experiment, two-dimensional depth· integrated (or so called depth-averaged) mathematica1 models have bcen widely adopted in the last two dccades. Because lhe vertical turbulent mixing evens out the vertical non-uniformity within a limited water column rather quickly and thus contarninants are distributed fairly uniformly over a water column, either finite difference methods or finite element methods can successfully simulate and predict tlow and contaminant transport in shallow well-mixed Oow regions with different time and length scales (Heaps, 1969; Kuipers and Yreugdenhil. 1973; Yu and Zhu, 1993). However, most of these mathematical models used in practical problems merely take the turbulent viscosity and diffusivity as constants or given by simple phenomenological algebraic expressions, which are in a great degree estimated in terms of numer ical modeller's experiences. The turbulent viscosity and diffusivity are not fluid properties but really depend on lhe turbulence structure; they may vary strongly from one point in the tlow to another and from one tlow situation to another. ln well-mixed regions the horizontal distributions of waste heat and contaminants are of great interest so that lhe depth-integrated turbulcnce models, especially the advanced two-equation closure models, are found to bc accurate enough to satisfy the demands of engineering practice (Rodi et ai .. 1981 ). The results

Manuscrípt received: Aprii19G5. Technical Editor Carlos Alberto CarraS(;() Altemani

339 J . of the Braz. Soe. Mechanical Sciences- Vol. 18, December 1996

gíven by the turbulence models, together with the depth-integrated continuity equation, Reynolds cquation. and scalar transport equation can provide details for the horizontal velocity components, pressure (water-depth) field, temperature and (or) concentration field, and can also provide the distributions of the corresponding turbulent physical quantities in the domain with engineering scale.

Unfortunately, the 'standard' turbulence two-equation models cannot be applied directly to depth­integrated computation and simulation. The necessary approach is to characterize lhe local state of turbulence in the sense of depth intcgration by two turbulcncc parameters, and to relate the eddy viscosity and diffusivity at each point to these parameters. The landmark of depth-integrated turbulence two-equation closure models is depth-integrated (é - E model ( (é: depth-integrated turbulent ldnetic energy parameter, E: depth integrated dissipation rate parameter of turbulent kinetic energy) established by Rastogi and Rodi ( 1978), which was stemmed from lhe 'standard' k - E model (k: turbulent kinetic energy, E : dissipation rate of turbulent kinetic energy, Jones and Lander 1973). ln 1988, Yu and Zhang provided another depth-integrated second-order closure model K- w ( w : deplh­integrated vorticity fluctuation pararneter of turbulence) by adopting the reviscd version of k-w model (w: time-mean-square vorticity fluctuation of turbulence, see Jlegbusi and Spalding, 1982). The numerícal investigation and application of (é - w model as well as the possible improvement and development of deptb-integrated two-equation closure models are presented in this paper.

Nomenclature B = wtátn ot nver or cnannel,

jet wldth, grid point of control volume

b = souroe term ln p'- eq. C = empírica! constants.

Chezy coefficlent E = empirical constant e = empirical constant G = production of turbulent

klnelic energy g = gravitational

acceleration H = local water depth, jet

penetration or eddy height

h = localtranqull water­depth or characteristic water-depth

J = ftux k = turbulent kinetíc energy L = eddy length or

turbulenoe length scale P = production due to

bottom shear or additional source term

p = pressure a ftowrate q = discharged ftowrate S = souroe term, constant

part or coefficlent of souroe term

T = temperature U = velocity for cne-

dimensional situation u = velocity ln x- direction v = velocity in y- direction w = time-mean-square

vorticily ftuctuation of

turbulence, vorlicity ftuctuation of turbulent

x = ()(H)rdinate y = ()(H)rdinate, vertical

distance between grid polnt and wall

Z = etevation a = water-depth correction

coefficienl r = turbulent eddy diffusivlty 8 = grid span or Kronecker

delta e = dlssipation rate of

turbulent kinetic energy 1< = kárman constant ll = dynamic viscoslty v :: kinematic viscoslty p = denslty a = PrandtJ number or

Schmidt number t = shear stress + = scalar q~ = indination of rtver bank n = vorticlty of mean

movement

Subscrlpta: B = submerged outtet b = bottom or river bed c = constant e = jet lnlet, east grld point

eff = effective f = friction I = eddy, co-ordinate

posítíon or index of tenser (i=1, 2, 3)

j = ()(H)rdinate position or index oftensor 0=1, 2, 3)

k = turbulent klnetic energy kv = additional source term in

k-eq. dueto effect of vertiCal velocity gradient

m = momentum, molecular n = north grid point p = first grid node near

banks, first order term of sourceterm

s = south grid polnt, water­surface

t = turbulenoe w = wall, west grid polnt, time­

mean-square vorticily ftuctuation of turbulence

wv = additional source term in w-eq. due to effed of vertical velocity gradients

x = ()(H)rdinate y = co-ordinate E = dissipation rate of

turbulent kinetic energy E v = additional source term in

k-eq. dueto effed of vertical velocily gradients

+ = scalar ll = dynamic viscosity O = inlet of channel or rtver • = bottom frictlon

Superscrtpts: • = empírica! constant - = depth-average defined

by eq. (8) - = depth integratlon (depth

average) defined by eq. (5)

, = correction

L Yu et al.: Recent Advances ano Appticabon of 1·umu1ence ... 340

Model Formulation

The goveming equations for steady and unsteady two-dimensional depth-integrated flows in natural waters can be expressed as

(I)

oHu + oHu2

+ oHVU = -gH ollh +2~[veffH ou] +~[vcffH (ou+ õV)] 0t ôx. 0y ôx. ôx. Re ôx. ()y Re ()y ôx.

(2)

(3)

where ób stands for the variation of local tranqui l water-depth h due to the motion of water, R e and prepresents the Reynolds number and density, t

5x, ts and tbx• tby are the shear stresses on the

water-swface and the river bed respectively, usand V5 denotes the surf~ce velocities in the x- and y­directions; v eff = v+ v

1 is a depth-integrated effective viscosity, with v, being the depth-integrated

turbulent kinematic viscosity and v being the kinematic viscosity of the Ouid. Correspondingly, the temperature and concentration of pollutants are govemed by the following scalar transpor! equation:

where om stands for the m_plecular Prandtl number for temperature transpor! or Schmidt number for concentration transpor!, r+ represents depth-integrated turbulent diffusivity for temperature or concentration respectively. The symbol '-'denotes the depth integrated quantities, by which tl!e depth· integrated velocity components ü, v in x- and y- directions and the depth-integrated scalar ~ in Eqs. (2)-(4) can be defined by the following relations:

(5)

where Zb and 4+H stand for the bottom elevation and water-surface elevation respectively (see Fig. I). The depth-averaged turbulent viscosity and diffusivity ~e not lf1!e depth-integrated quantities in the sense of strict mathematical definition by Eq. (5). When v, and r+ are multiplied by the relevant gradient of tbe transported quantity, the expressions of the depth-integrated turbulcnce shear tress and mass flux will be satisfied respectively (see Eq. 8). A linearized assumption has been adopted in the numerical process to treat tbe source terms Smx, Smiand s, in control volume, produced by submerged discharge:

(6)

341 J. of lhe &az. Soe. Mechanical Sdences • Vol. 18, Deoembef 1996

where S0 stands for the constant part of S, SI!. represents the coefficient of the first-order term of S. The expressions of S mx and S my and S+ will be specified in the diffusion simulation section,

further on the t.ext.

Fig. 1 1tu1tratlon of 1 Ba1lc Control Volume

Equations (I) to (4) were obtained by simply integrating the fundamental three-dimensional continuity equation and Reynolds equation from the bottom to the water-surface. One should note that the local water-depth H cannot be treated as a constant in the process of integration, since the surface elevation Z, +H is usually function of the position, for steady flow, and of both the time and the position for unsteady flow (see Fig. I). The last two terms on the right hand sides of equations (2) and (3) appeared because of the variation of bottom topography and surface elevation, which are referred as a kind of additional source terms produced by Leibniz's integration formula ln each vertical water column. and wíll be analyzed in the Details and Skills ofNumerical Computation section, later in the text. Toe variation of the bottom topography and surface elevation needs to be included in the model and dealt efficiently. This can be achieved by introducing water-depth correction coefficients into the discretized formulas (Yu, 1991). The surface stresses tn and t 1 ycan be det.ermined by empirical expressions in the computational domain; with a spatial scale relatively small, their effect can be neglected. The terms tbx and tby stand for the bottom shear stresses i o the x- and y- directions, respectively; they are calculated by relating them to the depth-integrated velocity components via

'bx C -J-2 -2 -= ru u +v p

(7)

After determining the depth-averaged turbulent ltinematic viscosity v1

the equations (I )-(4) forma set of closed partia! differential equations that simulate and predict the depth·int.egrated turbulent flow and the corresponding transport field. ln equation (7). c, is an empirical friction factor which usually is equal to g!C1 in rivers, where C is Chezy coefficient.

Depth-integrated 1<: - t Model

The turbulence two·equation c losure model that found widest application and was tested most extensívely is the k- c model, which characterized the local state of turbulence by two parameters: the turbulent kínetic energy k and the rate of its dissipatíon c . Since 1978, the model has been developed to describe the depth-integrated open~hannel flow (Rastogí and Rodi). ln analogy to the original k· c,

L Yu et aL: Recent AdVances and Application of Turbulence ... 342

t!]e local depth-integrated state ofturbulence was assumed to be characterized by the turbulence_energy k and dissipation p_arameters e respectively, and the deplh-integrated turbulent stresses 'ti,jMd

contaminant fluxes J+i can lhen be related to lhese parameters via

:ri.i _ (õü; õüi) 2 -- =v -+- - - k~ ·. and P t ÕX · àx· 3 l,J

J I

(8)

with (9)

where Cl-l,!Uld a+ are empiriEal ~onstants. Equation (8) includes the depth-averaged turbulent visco~ity v, and diffusivity r+ . k and its ê dissipation can be calculated by two more differential equations:

( lO)

( 11 )

where G stands for lhe production of turbulcnl kinetic energy due to interactions of turbulent stresses wilh horizontal mean velocity gradients which can be expressed as follows

(12)

The coefficients C., <;. ak and a, in (to) and (li) are empirical constants. F..quations (I O) and ( 11) can be considered as deplh-integrated forms of the three-dimensional k and E equations, when a gencration term pkv has been added in lhe K-equation and p &Yin thc e -cquation to account for turbulence gencration due to the vertical vclocity gradicnt, whic~ occurs ma.Jnly near the bed. The gcneration G dueto the horizontal velocity gradient is practically v1 (ôn/ôy) in lhe deplh-averaged case. Since the additional generation term P"" and P~vdepend strongly on the bottom friction, these tenns are related to lhe bottom friction u• in thc following way (a detailed discussion is contained in Rastogi and Rodi, 1978):

(13)

where lhe local friction velocity u • '" Jc r ~2 + v

2) • The empirical constants C., and C for channel

flow determined by Rastogi and Rodi ( 1978). with the turbulent viscosity measured by L~ufcr in 1951 are:

c = -·­"JCr _ I F

and cc - 3, 6~..,c11 C r

Thc values of the empírica! constants adopted are lhe ones of the 'standard' C11 =0.09,0'.=0.9, C,• J.43, ~=1.92 , Olt -= 1.0 and Oc=l.J.

k - t model :

343 J . of lhe Braz. Soe. Mechanlcal Scienoes • Vot 18. December 1998

This model has been applied to simulate the open channel flow (Rastogi and Rodi, 1978), the symmetric expansion flow (Chapman and Chio, 1982; Chapman, 1983) and the flow and pollutant diffusion in ri ver (Rodi at ai. , 1981 ).

Depth-integrated (( - w Model

At present, the depth-integrated mathematical models are very popularly uscd in water regions with large ratio of width to depth for engineering simulation and prediction (Heaps, 1969; Kuipers and Vreugdenhil, 1973; Yu and Zhu, 1993, and so on). However, most ofthese depth-integrated models only adopted low-order 0-equation turbulence models. Before 1988, users who intended to adopt the advanced two-equation closure for modeling their goveming equations had no other choice, as there was only one depth-integrated turbulence two-equation model originated from 'standard' k- E model (the J(- t model). On the other side, two-equation closure theory and corresponding models h ave made great progress; severa! effective models and their newly developed versions have been applied in different engineering arcas. ln 1988, the depth-integrated turbulence J(- w model was established by Yu and Zhang on lhe base of the new versíon of the k-w model. reported by llegbusi ( 1983, 1984). The revised k - w model demonstrated the sarne order of accuracy as the 'standard' k. - E mode~ but bener accuracy than the k- E model for some computational examples. This revised model can be expressed as follows (llegbusi and Spalding, 1982):

{ àk ak ) à ( ôk ) 1/2 - + u . - = - J.L ...... - + pPk -C plcw a J~. ~- ~~~ . ~

J J J

(15)

and

(16)

whcre n and L stand for the vorticity of mean movement and length scale of turbulence respectively, the cffective diffusivities ~-'em• ~-'effw in the k- and w- equations are:

llefJw = J.1 + J.L, I C1 w ( 17)

and thc empirical constants for turbulence model C~, ale, aw. C1,.., C2w,C'Jw and C3w in Eqs. (15}­( 17) are equal to 0.09, 1.0, 1.0, 3.5, 0.17, 17.47 and 1.12, respectively.

The local depth-integrated state ofturbulence can be characterized by the foUowing two turbulence parameters: the turbulence kinetic energy K and the intensity of vorticity fluctuation w in depth­averaged sense; the turbulent viscosity/diffusivity can be expressed as follows (Yu, 1988):

- i:..-:-:-1/2 - i:..-:-: - 1/2 r - I (18) 1!1 = p11.w v1 = 1\.W 1 = v1 CJ+ Here I( and w are determined from following lhe two transpor! equations (Yu and Zhang, 1989):

L Yu et ai. : Recent Advances and Appllcalion of Turbulence ... 344

+C v Hl__?_(õV- õú)l2 lw t ôy ÔX ôy

(20)

The additional source terms Ph and P wv in Eqs. (19) and (20), produced mainly by the vertical velocity gradients near the bottom, can be expressed as follows (Yu, 1988}:

(21)

By using Laufer's open channel experiment result, i.e . il1 = 0.0765U•H to detennjne the empírica! constants Ck and Cw, it can be obtruned (Yu, 1988):

C = 4726C /C 312C314 w 2w 11 f (22)

Wilen Eqs,_ ( 19) and (20) are SQlved for _!llodeling the submerged discharge, corresponding source terms Sk and Sw similar to Smx , Smy and s, in Eqs. (2)-(3) should bc added to express the effects of the local diseharging outflow on the turbulence fields. Figure 2 shows a plane configuration for a waste heat díscharged from one side into a rectangular water cbannel. Tbe comgarison of computational results with experimental data in diiTerent momentum flux ratios v; s; ;u0 s 0 has been performed carefully (Yu and Zhang, I 989). The computational results by use of depth-integrated {(- & and

I< - w models configured the dimensionless eddy height 1-1;180, re-attacbment length H i/Li of recirculation eddy and dimensíonless jet width (H0-H;)/Li (see Fig. 2). 1t has bcen found that when the outlet width or flowrate of discharge is relatively small, which is a situation frequently occurred in environmental enginecring, thc rcsults computcd by thc 1< - w model are closer to experimental data than thosc simulated by the {(- & model. Because the jet width H0-H; can bc regardcd as a horizontal fine sourcc which produces the pollutant plume in tbc lower reach of thc discharge outlct, to fmely determine H0 -H; has a significant and dircct cffect on the computational accuracy for simulating and prcdicting the sprcad and distribution of thc plume. Figure 3 presents the computed vclocity vector field, stream líne, isotherm, water-depth contour, longitudinal vetocity, depth-averaged turbutcnt parameters ( {( and w) as wett as deptb-integrated turbulence viscosity ii1 for one of the computed nine cases

Fig. 2 flow ConflguratJon

B&- channel width B.- jet width U0 - channel mean veloc/ty V, - discharge mean velocity H0 - maximum jet penetration H,- eddy heighl L, - eddy length

345 J. of the Braz. Soe. Mechenlcal Sclences. Vol. 18. Oecember 1996

b) Strum lilld

~-----~-= c)bolbcnns

d} Watcf-ciqldl 00111ours

Fig. 3 Compute<! Rn utta by Uae of I( - w Model For One Olacharge Caae

Simulatíon For Flow and Thennal Diffusion in the Rhine River

General Descrlption

Recently, lhe authors further investigatcd and applied the depth-integrated k - w model to simulate lhe waste heal discharge for a practical problem in engineering: the cooling water was simultaneously discharged from lhree submerged outlets of an electric power plant into lhe Rhine River. The basic schematic diagram for lhe computation has been presented in Fig. 4. The fluvial dynamics parameters oflhe computational river which were obtained during lhe observation period are the following:

The averaged flowrate Q0 = 1550 m3/s;

The average width 8 0 = 2~0 m;

The average water-depth h = 4.85 m;

The friction velocity u. = 0.068 m/s;

The water temperature T 0 = 7Cf C, for this computation • stands for the temperature T;

The averaged si de slop: about I: 4.7:

The flowrates and temperatures of thermal discharges from three outlets q., T .. q2, T2 and cu. T3 are 4.0 m3/s, 20°C, 1.75 m1/s, J8.4°C and 6.0 m3/s, 17.9°C respectivcly .

• .. lt,l4t • • 'lt. ' )1 •• 4)0,)1 )'

·~

' Fig. 4 Olacharge Schematlc Olag,..m

L Yu et ai.: Recent Advances anc:t Applicatlon of Turtlulence ... 346

The number of nodes in lhe computational grid on a 1,217 m Jong segment is 256 along lhe ri ver direction and 66 across lhe river direction, respectively. ln lhe longitudinal direct.ion (flow direction), lhe grid adopted piecewise unifonn distribution and the grid !ines near the three discharging outlets were denser than others; in the cross river direction, the grid spans near lhe two ban.ks were finer than the ones in the center part of the river. The 'standard' or 'original' empírica! constants of both ~-E and f(- w models were adopted witbout any alteration, except lhe coefficients Cw and Ct , whidl had to be calibrated in terms ofthe practical problem by the authors.

Boundary Condttlons

The boundary conditions include those at the ri ver inlet and outlet, those on the river banks and at the submerged outlets. At the river inlet, the velocity profile was taken to be the observational data, the temperature was set as the natural water temperature (70°C). the turbulent parameters for f(- w model were expressed by the empírica! expressions k = Jo.6 U 0 U ~ I C~ and w • u 0 u • 1 o.6 c • H , for the K- E model, & = u 0 u ~ I H .

At the downstrearn section of the ri ver, the zero gradient condition of ali variables except v was · ed · aü a'f ak' aw a& d · d th h n · th tmpos , 1.e. _ .. _ = _ ., _ ., _ , an tt was assume at t ere was no ow m e

õx õx õx õx õx

cross-river direction, i. e. 9 . On !_he river banks, the impervious condition wt\re imposed on the normal velocity, and the temperature T was assumed to satisfy the adiabatic condition. ln addition, tbe longitudinal velocity o and turbulence pararneters f( , w or E at lhe first grid node P near the banks satisfied the wall functions of the corresponding turbulence models (llegbusi, 1984, Yu and Zhang, 1989):

(23)

(24)

where 1C = 0.4, E= 0.9 for smooth wall and Tw"' Tp in the simple flow situation. For the computation ofpractical engineering problems, E needs to be calibrated from site data.

Apart from the additional source terms P""' P wv or P tvresulting from the vef1ical integr~tion, there are some extra source t~rms resulting from submerged oujlets, su~h as t~e Smx and Smy in the momentum equations, S+ in the transport equation and S~c and Sw or s, in the corresponding turbulence equations. These extra source terms can be dealt with as if they were part of the boundaries of a computational domain. ln our computation. the submerged discharging outlet (see Fig. 5) was approximated by a vertical line source at the nodal point in the momentum equations, transport equations and the pressure-correction equation (p ·- eq.) derive<! from the continuity equation and tini te volume approach as shown in the literature (Patankar, 1980).

N

E

Fig. 5 Conflguratlon of • Submtrved Outltt 1nd Nod•l Gr1da Amngtment

J . d tlle Braz. Soe. Mechanlcal Scienon • Vol. 18. Dec:embef 1996

Since lhese outlets were treated as part of the boundary for the computational domain when lhe discretization was carried out, lhe boundary conditions for submerged discharging, as the ' inner' boundary oonditions, could be derived according to the conservation principie of various dependent variables at the control volume centered ata grid point B (Yu, 1988). For exarnple, the pressure correction equation now conta.ins a source term b, whlch is called 'mass source' in the literature (Patankar, 1980 ) . For ali ceJls, where no pollutant outlets are located., this source tenn is set to zero. On the other hand, if lhere is an outlet with flow rate <Is in the control volume B, this source tenn is lhen set to be equal to q8 . ln the discretized form of tho depth-integrated mathematical model, lhe additional source term appeared as ScóxóyH + SpÃX<lYH+ (Yu, 1988). With óxand óy being x­and y- direction widths ofthe control volume respectively, in the pressure correction equation, the steady discharge of the vertical linear source at nodal point B should not produce the pressure correction p·. Therefore, we set lhe pressure correction + = p' = O, lhe 'source tenn' at nodal point B is only ~ = ScâxâyH and the linearized coefficient s. of constant term should be equal to q..a / 6xóy H (see Fig. 5). The derivations for lhe ' source terms' in other equations are similar (see Yu, 1988 for deta.iled derivation) and are listed below without further derivation:

u - eq . S =-~i!!_ I PIJ 2 H ây(óx)w

S = - .!_ .9.!!_ I PI+IJ 2 H 6y(ôx)e

(25)

v - eq . (26)

s = -.!. .9..!!. ---Pi.J+I 2 H àx(óy)n

(27)

p'-eq. (28)

~-eq. (29)

k -eq. (30)

E - eq. S = qoEo s = Go (3 1) CIJ ~y ci,J Hàxtly

w - eq. S = qBwB s = Qs (32)

ÇIJ Hàxày ci,J Hàxày

where <Is stands for lhe discharge strength ofthe verticallíne source at point B; (ax) e' (Sx) w express the x-direction spans betwcen B and the east and west control volume faces; (Sy) n• (Sy) are lhe y­direction spans between B and the north and south control volume faces respectively (see ~ig. 5). The equations (25)-(29) have been successfully taken as the outflow boundary conditions for simulating

L Yu et ai.: Recent Advanoes and Application of Turbulenoa. .. 348

severa! immersed outlets discharging waste heat and contaminants into the south estuary of the Yangtze River (Yu and Zhu, 1993), while an simple algebraic formula was used to express the turbulent viscosity and diffusivity oftidaJ currents. The Eqs. (30)-(32) were then firstly published, wbere K8, &8 and w8 could be detennined by the empirical expressions used in tbe river inlet.

Details and Skllls of Numerical Computation

Cholce of turbulent dlffualon coefflclent

The given constant vaJues ofthe coefficients C and Cw in Eqs. ( 14) and (22) were determined by using Laufer's experimental result in small wat~r channel. rn order to obtain tbe usable values adaptable to model natural ri ver, the dimensionless diff\.isivity e• = v l ( U. ho •) in the center part of the ri ver reacb was adopted, where tbe flow was assumed to be uniform. Fisber et aJ., ( 1979) reported that the vaJue of the coefficient e• was in the range of 0.4 - 0.9 for a number of different rivers. Rodi at ai ( 1981) adopted e• = 0.6 to simulate the Rhine Ri ver and obtained a fair good agreement with the field measurements. ln our computation, e• = 0.6 was also adopted.

Computatlon of Bed..Shear Stl'8$a

The two components of bed-shear stresses in momentum Eqs. (2) and (3) should be computed by use of the quadratic friction law considering the etTect of tbe bank slope (Rodi et ai. , 1981)

(33)

wbere ~ stands for tbe inclination defined in Figure 6 and the friction coefficient Cr is not equal to a constant vaJue similar to the situation in the rectangular channel. The Cr increases towards the edge of the shallow water regions. The variation of Cr can be evaluated by using Rodi 's proposition presented in 1981.

1----y .,_. ___ B

012----.j

H

·~~~~~., "---. Pw/2

Fig. 6 The cro .. -sectloll of a Typlc:1l Rlver B1nk

I Treatment of Addltlonal Source Tenns Produced by Verticallntegratlon

Since the variation of bottom topography and (or) water-surface are unavoidable, the terms , produced by using Leibniz formula to integrate the three-dimensionaJ goveming equations should be

analyzed carefully, as the effects of these terms could not be neglected at ali. Our numerical investigation has shown that the produced terms whicb have a direct effect on computational results are (Yu, 1988):

in u- equation (34)

J. of lhe Braz. Soe. Mechanat Sdences • Vot 18. Deoember 1996

(35)

where ~ and v, stand for lhe surface velocity components in x- and y- directions, respectively. ln lhe long and straight ri ver reach witb trapezoidal cross-section, while tbe variation of water-deptb mainly occurred along the transversal dlrection, if the variation of water-surface elevation relatively to tbe water-dcpth is small enough, Eqs. (34) and (35) can then be expressed approximately in terms of following suggestions (Yu, 1988):

The gradients of depth-integrated velocity components wilh respect to the x- and y- directions can be substituted for surface velocity components and

The tranquil water-depth h can substitute lhe dynamic water-depth H.

ln lhis case, Eqs. (34) and (35) can be simplified as follows:

in u- equation S x = -i!, ( ~ + :) : (36)

in v- equation Sy = -2;:l, : : (37)

Af\er we approximated lhe depth-integrated vel<><:_ities ü and v to t}le surface velocities in Eqs. (2) and (3) and calculated the turbulence viscosity v, and díffusívity r. by Eq. (18), tbe equation system ( 1 )-( 4 ), (I 9), (20) or ( II) became closed. These equations formed a set of two-dimcnsional non-linear partial dífferential equations with ellipticaJ sítuation in !!Je spatial scope. However, as t}le unknown variables O, 9, H. ~. w or & and the cffective viscosity v cff as well as the diffusivity r+ are ali coupled in a highly non-Iinear fashion, an adequate numerical solution must be sought for these equations.

ln our computation, the SIMPLE-c (Semi-lmplicit Method for Pressure-Linked Equation­consistent, Doormaal and Raithby, 1982) algorilhm for iteratively solving ali momentum and transpor! equations as well as pressure-correction equation has been adopted. Figures 7 and 8 clearly show tbe comparisons of variation of 'mass source', in which the 'mass source' excluding the two terms (36), (37) and including these terms on the right sides of the momentum Eqs. (2) and (3) has been recorded. respcctively. lflhese two terms were neglected ín the goveming equations, the iteration process would show that larger fluctuation of 'mass source' appeared in the initial iteration stage and conversely, the iteration rapidly tendcd to be convergent only in limited 45 iteration steps

1 o·• '--..-..-...-...--r--r--r--r--r-., o 20 40 60 80 100

Fig. 7 Maxlmum Mau Soure• ln Control Volume

L Yu et ai · Recent Advances and Application of Turbutenee .• 350

10 -S~.Sy~O

- ·-Sx.Sy=O

....., rj\ \ \ J' 10'1 \ \,I \

....

10'' o 20 ~o 60 ao 100

Fig. 8 Ma .. Souret ln AllttMt Computatlonal Domaln

Water-Depth Correctlon Coefftcients

Since the variations of bortom topography and ( or) water-surface are generally quite large, the variations of battom tapography and (ar) surface elcvation should be included into lhe mathematical model. This can be achieved by introducing the water-depth correction coefficients, which were defined by Yu in 1991 as:

a .. = H · /h= (b .. +.M .. ) / h IJ \) IJ lj > (38)

where H;J, h,J and 6h; i stand fo! the local dynamic water-depth, the local tranquil water-depth and the variatton of H;j respectively; h represents a characteristic water-depth, usually equal to the average tranquil water-depth ofthe main cbannel for a river reach. The water-depth c.orrection coefficients are functions of space, and also of time for unsteady flow. The variation of a 1 j provides a consistent means to deal with the variation of water-depth and boundary geometry', for example, «i j = O corresponds to land or islands and non-zero ai . represents a region occupicd by water. The water­depth correction coefficient technique has becriJadopted to treat thc variations of bottom topography and waler-surface for modeling the tidal currents in the estuary ofthe Yangtze River (Yu and Zhu, 1993). By using finite volume approach, the water-dcpth correction coefficients can be tlassified according to twa possible locatians af a grid point within a control volume. Those 8.SS()tiated with the center point of a control volume. the values of which approxiroately represent the fraction of volume occupicd by fluid in the contrai volume. and those associated wilh the ccnter points of contrai faces of a contrai volume, the v alues of which approximately represent lhe void fractions of correspanding contrai faces.

Computatlonal Results

Ali thc unknown variables were discretized with spatially staggered grids. The depth-intcgrated velocities o and 9 are computed at the faces of each main control volume, whereas the water-depth H. the temperature T and the turbulence paramcters 1<, w, or tare computed at thc center of each cell. The SIMPLE-c algorithm (Doormaal and ~aithb~, 1 982) was adopted to couple the vectorial quantities o , 9 and the scalar quantities H, T, 1<, w , or & • With the power-law scheme (Patankar, 1980) having lower numericaJ dissipation and highcr numcricaJ accuracy being used for lhe convective and ditTusion terrns, Eqs. (I )-(4), (19) and (20) or (li) are discretiz.ed with the finile volume approach. Iterations were perforrned at each solving step with an initial guess solution of the pressure (water­depth) field. The velocity componeots o and 9 corresponding to the gucssed pressure field are solved from the momentum Eqs. (2) and (3) by usíng the so-called line·by-line method (Patankar, 1980)

351 J. of lhe &az. Soe. Mechanical Sdence$-V~ 18. Decembef 1996

associated with the wall functions. The pressure (water-depth) values are corrected with the line-by­line method again, from the pressure correction equation. ln our model the used pressure correction equation and two velocity-correction equations, which calculated the corrected velocities, are similar to those in the literature (Patankar, 1980) and, therefore, have not been included. By using the corrected velocities and wall functions, the temperature transport equation and turbulence parameter cquations could be solved, and the turbulent viscosity and diffusivity needed to be renewed at the end of each iteration. After that, the corrected pressure (water depth) is taken as a new guessed pressure (water-depth) and the iteration repeats itself until a pre·specified convergence is achieved. Due to the tedious process ofthe discretization ofthe governing equations and their low relevance to the theme of this paper, the detailed numerical procedures have not been included.

The computational results from IC - w model were compared with the observed depth-integrated longitudinal velocity and temperature distribution, and also compared with the corresponding numerical results from constant coefficient viscosity and diffusivity model and IC - & model, respectively. These computations were performed by adopting a general-purpose computational prograrn for elliptic situation ínvolving the water-deplh correction coefficients (Yu, 1994) under the sarne convergente indicators, which are the sum ofthe absolute values ofthe 'mass source' in the whole computational domain and the maximum mass source of one cell within the domain. The computations were finished when these indicators became smaller than 2.5% for solving the pressure correction equation. The max.imum residue for solvíng temperature transport equation between two adjacent iteration cycles was less than lx 106• The single block correction technique (Prakash and Patankar, 1981; Patankar, 1981) was adopted for accelerating the convergence of these computations. ln order to investigate the newly developed turbuJence k- w model, different turbulence models such as constant coefficient viscosity and diffusivity model and IC -E model were adopted in different side wall situations (vertical side and sloping bank). Figure 9 presents the computed local now field near the thrce submerged outlets, in which reasonable flowing situatíons and details in the ri ver have been shown clearly. The distributions of computed longitudinal velocities and temperature compared with the field data have becn shown in the Figs. lO and 11 , respectively. For convenience, the x- and y- coordinates in these two figures represent the lower reach direction and cross-river direction of the ri ver. ln terms of the results from the temperature equation, the isoth.erms have been presented in Fig. 12. which clearly illustrates the development of poUutant plume in the arcas both near the three submerged outlets and of lower reaches of these outlets. From the objective analysis and comparison, the following conclusions have been obtained (Yu, 1990; Yu, 1993)

Without considering the etTect ofbottom topography on velocity and temperature distributions, for example, only adopting the simple and rough flat-bottomed approximation, the thermal plume will concentrate on the ri ver bank, the computational result does not coincide with the measurements:

By using depth-integrated IC - w and ~- Ê models, the comyuted temperature distribution is different; the results computed by the newly developed lê - w model are better than the ones computed by the lé - Ê model and

The temperature distribution using the model with constant coefficient diffusivity has a relatively larger error in comparison with field data; this is very evident even at the lower reaches r ar from the discharging outlet.

Fig. 9 Loeal Flow Fleld Neer Th,.e Subme!11ed Outlets

L Yu et ai .: Reoem Advai10N and Application of Turtlulence ... 352

!11

/0

o

Fig. 1 O Dtatrtbutlon of l ongitudinal Veloc;ltln

Fig. 11 Dlatrtbutlon ofTemperatunt

o site data

-- r.; modtt

-·- k-&mOtkl ----· k-&model

(vertical wall)

• • • • • COIU .. diffu.. coeff. lfiOdel

5'C

--- ·- · --·-·~· -- ·- ·-· ----- · -

---------7.0-------7.0 ----7.0 7J 7J 1.J ----

F lg. 12 lsottlerma

Because of the poor infonnation about constructions near the banks and on the bottom topograpby, such as the locarions and sizes of docks, obstacles and pipeline systems in the lower reach of outlets, and also Jack of detailed data about the side slopes, which in our computations had to be simplified as an averaged single vaJue, there ex.ist rather large differenoes between the computed velocity profiles and site data in the near bank regions at the sections located at 430.621c:m and 430.88 km, respectively.

I Discussion

lt is recognized that up to the prcsent. the advanced and practicable method for illustrating turbulent mean behavior in two-dimensional depth-integrated simulation and prediction is to adopt two-equation turbulence closure models. The newly developed depth-integrated turbulence two­equation closure IC - w model is available to be used to compute turbuleot mcao bchavior i o engineering. The authors' computational results have shown that thc temperature fields using depth-

353 J. ot the Braz. Soe. Medlanlcal Sciences- VoA. 18. Oecembel' 1998

integrated R- w and K- Ê models have some important ditTerences, that is, in the case ofnarrower widths of jet outlet or smaller jet flowrate, the results computed by the fé- w model are better than the ones computed by the I<- t model. ln the presented computatiooal example of practical engineering significance, the widths and flowrates of the three outlets of cooling water discharge are smaller relatively to the width and global flowrate of the river. Hence, it can be concluded that the jet widths simulated by the K - w model should be more correct than the widths computed by the tC-e model and the pollutant plume in the lower reach simulated by the K- w model, including the position and temperaturc distribution, would have a reasonably higher accuracy.

Possible lmprovement and Development

ln the pas1 18 years while modcllers tried to adopt turbulence two-cquation model to simulate and predict flow and pollutant transport in engineering, the development and application of depth­intcgrated two-equation closure models was not satisfactory indeed. Recently, Booij ( 1989) pointed out that in Rodi 's depth-integrated K-Ê model the turbulence production due to horizontal velocity gradients and the production dueto the presence ofthe bottom were not weighed correctly; this had lcd to a much exaggerated bottom influence. A modifted model suggested by Booij (1989) was shown to yield better results for modeling mixing layer between a ri ver anda harbor in a se ale model, in which a man-made weighing to regulate these two parts of turbulence production has been used. lo fact., it is still doubtful that the artificial weighing of the turbulence productions dueto horizontal and vertical velocity gradients is in lack of sufficient theoretical basis. For the K- w model developed by the authors, further numerical and experimental investigations are also needed. Moreover, the applicability of a few well-behaved. valuable versions of turbulence two-equation models in the depth-integrated sense has not been invcstigated sufficieotly, even after the utilization of some of these models has become quite popular in various industrial departments. On the otber hand, englneers still tend to adopt the traditional constant viscosity and diffusivity with those coefficients being roughly estimated empirically. Such questionable estimation combined with obsolete models doubtless led to large differences between numerically simulated results and observed ones. The severe gap between the theoretical research in turbulence modeliog and its practical developments has further retarded improvement of the highly-accurate depth-integrated turbulence two-equation models and their application to natural waters. The failure ofthe theory to meet practical demands should be rectified as soon as possible. lt is expected that further research, investigation and application of advanced turbulen~ two-equation closure models will become very popular among researchers and engineers in the arca of depth-integrated modeling. At present. it is possible to develop some new depth-integrated turbuleoce two-equation closure models based on some oewly developed versions, which have corrected some shortages of the ·standard' k - e model. The authors expect thatthe newly developing featurcs for the existent K- Ê and K- w models can set up a higher standard for numcrically modeling turbulent flows and pollutant transport phenomena in tenns of computational efliciency, a1gorithm extensibility and universality, as well as model robustness.

Acknowledgment

The support ofCNPq (Bra.zilian National Council for Scientific and Technological Development) through Process No. 30 J089/94-9(NV) is gratefully acknowledged. The authors also thank Prof. Marcius F. Giorgetti of the School of Engineering at São Carlos- USP for h is useful assistance in the preparation of thls paper.

References

Soo, R., 1989, Dcpt.h-Averaged k- a Modellng, ln: Procecdíngs of IAHR XXUI Congress, Tcchnical Section A, Turbulence in Hydraulics A, 199-206, Ottawa, Canada, 21-2S, August.

Chapman, R.S. and Kuo, Y. C., 1982, "A Numcrical Simulation ofTwo-Dimensional Separated Flow in a Symmetric Open-Channcl Expansion Uslng thc Dcpth-Jntegrated Two-Equation ( k- ~) Turbulence Closure Modcl», Dcpl ofCivil Engineering, Virgínia Polyteebnic lnstitute and State University.

L Yu etal.: Recent Advances and Applicalion of Tulbulence 354

Chapman, R.S., 1983, "Two-Equation Oepth-Averaged Turbulence Closure for Modeling Geometry-Dominated Flowsft, Environmental Laboratory U.S. Army Engineer Waterways Expcriment Station.

Fisher, H.B. c! ai ., 1979, "Mixing in lnland and Coastal", Acadcmic Press. Heaps, N.S., 1969, "A Two-dimcnsional Numcrical Sea Modcl", Phil. Trans., Roy. Soe., A26S, pp.93-137. llcgbusi J.O. and Spalding. D.B., 1982, "Appllcation of a New Version ofthc k - w Model ofTurbulence to a

Boundary Laycr with Mass Transfcr", CFD/112/1 5. flegbusi, J.O., 1983, "Revised Two-Equation Model ofTurbulence", Imperial College ofScience and Technology,

Rep. CFD/83/S. Tlegbusi, J.O., 1984, "Proposal for Wall Function for Friction and Mass Transfer", CFD/84/12. Joncs, W.P. and Lander, B.E., 1973, "Tbe CalcuiBtion ofLow Reynolds Number Phenomena With a Two-Equation

Model ofTurbulcnce", Int. J. Hcat Mass Transfer. Vol. 16, pp. 1114-1130. Kuipers, J. and Vrcugdenhll, C.B .• 1973, "Calculations ofTwo-Dimensional Horizontal Flow", DeU\ Hydraulics

Laboratory Rep. Sl63, PartI. Laufcr, J., 1951 , "lnvestigation ofTurbulcnt Flow in a Two-Dimcnsional Channelft, NACA Repon 1053. Patankar, S. V., 1980, "Numerical Heat Trans fer and Fluid Flow", Hemisphere Publishing Corporation and

McGraw Hill Book Company. Patankar, S.V., 1981 , "A Calculation Procedure for Two-Dimensional Ellíptic Situations", Numcrical Hcat

Transfe.r, Vol. 4, pp. 409-425. Prakash, C. and Patankar, S.V. , 1981, "Combined Frec and Forced Convection in Vertical Tubes with Radial

Internal Fines", ASME J. Heat Transfer, Vol. 103, pp. 566-572. Rastogi A.K. and Rodi, W.A., 1978, ''Prediction ofHeat and Mass Transfer in Open Channels", 1 ofthe Hydraulics

Divislon, ASCE, no. Hy3. Rodi, W.A. et ai ., 1981, "Prediction offlow and Pollutant Spreading in Rivers". Transport Models for lnland and

Coastal Watcrs, Proceedings of a Symposium on Predictive Ability, Academic Press, pp. 63-11 1. Van Doormaal, J.P. and Raithby, GD., 1982, "Enhancemcnt oftbc SIMPLE Method for Predicting lncompressiblc

Fluid Flowsft, Numerical Heat Transfcr, Vol. 7, no. 2, pp. 147-152. Yu, L. and Zhang, S., 1988, "A New Depth-Averaged Two-Equation (i - w) Turbulent Closurc Model",

Proceedings ofthc Tbird lntcr. Symp. on Refined Flow Modcling and Turbulcnce Measurements, Tokyo, Japan. July 26-28 or J. ofHydrodynamics, Vol. I , no. 3, pp. 47-54 (1989).

Yu, L., 1988, "The Turbulent Model and Nurnerical Researcb for The Turbulent Transpor! ofThe Contaminants in Watcr Environment" (in Chincsc), Ph.D. Dissertation, Hohai.

Yu, L , 1990, "Applicar.Jon ofNew Depth-Averaged Two-Equation ( ~ -w) Turbulent C1osure Modcl to Numcrical Simulation for a Rivcr", Proceedings ofthc 7th Congress of APD-IAHR, Vol. 3, Beijing, Cbína, Nov. pp. 157-162.

Yu, L., 1991, "Two-Dimensional Variable Water-Depth Mathematical Model for Tida! Flows and Numerical Computation for Discharges of Waste Heat and Contamlnate ln South Estuary of the Yangtzc River", Proccedings of lntemational Symposium on Envirorunental Hydraulics, Oct. Hong Kong.

Yu, L. , 1993, " Establishment. Numerical lnvestigation and Applícation ofDepth-Averaged Two-Equation Turbulcnt Closure Model (i - w )", Procecdlngs of Annual Confcrence of Australian Mathematical Society, Wollongong, Austral ia, July.

Yu, L. and Zhu, S., 1993, "Numerical Simulation of Discharged Waste Heat and Contaminants into thc South Estuary ofthc Yangtze Rivcr", Mathl . Comput. Modeling Vol. 18, pp. 107-123.

Yu, L, 1994, " Eicmcnts of Numerical Computation and Simulation of Flow and Contaminants Transport", A Lecture Course, Oept. of Hydraulics and Sanitary Engincering, Sio Carlos School of Engincering, University of Silo Paulo, Brazil, March-July.

RBCM - J. of the Braz. Soe. Mec:hanlcal SCieneea Vol. XVIII- No.4- 1996- pp. 355-373

Added Mass and Damping on Bottom­Mounted Rectangular Cylinders Paulo de Tarso T. Esperança Sergio Hamilton Sphaler Umversidade Federal do Rio de Janeiro Departamento de Engenharia Naval e COPPE - Programa de Engenharia Oodnlca 21945 Rio de Janeiro, RJ Brasll

Abstract

ISSN 0100-7386 Printed in Braz.il

Looking for a major comprehension from the reason ofthe stroog variations in the hydrodynamic ooefficients of added mass and damping when a two-dimensional body oscillates vertical close to the free surfàce in finite depth we analyz.e lhis problem using eigeofunction expansions. Tt can make evident lhe parameters involved oo resonances and the beb.avior for low and high frequencies. lntroducing scparation ofvariables the bouodary valuc problem is transfonncd in ao eigenvalue problcm ofSturm-Liouville kind. This procedure was used by Yeung (1981) to detennine lhe hydrodynamic ooefficients for a truncated vertical cylinder with circular section piercing thc free surface oscillating in water of finite depth, Yeung and Sphaier ( 1989) to consider tbe interfercnce of vertical walls of a canal in the hydrodynamic coefficients, Mclver and Evans ( 1984) in lhe problem of a submerged vertical cylinder oscillating near the free surfacc and Esperança (1993) for the two­dimensional case. ln the present paper we confum the occurrence of negative added mass when the body oscillates close to the free surface. and obtain tini te values of the added mass and damping coefficient in dimensional fonn for lhe zero frequency limit as expectcd. The negative values of the added mass are related to two kind of resonant

1wave

modc.s occurring close to values ofthe frequency parameter, wave number times halfrec~angle width, m a, is equalto n1t and n1tl2. For tllese frequencies lhe darnping coefficie.nt is equal to zero or acbieves a ~ocal maximum value, and statio.nary waves cbaracleriz.e lhe internal flow. Negative addcd mass for the frcquency parameter going lO zero is obtaincd for shallow water, wben the cmergence oflhe rectangle is Jarger than 113 of lhe water depth. A shallow water approacb is also presented, allowing us to dcscribe the wave motion for small water deplh and clearance. The resu lts derived from the shallow-water solution is in good agreement with the one of the complete solution. Further we present lhe solution for the case of rectangular cylinder beaving on the frce surface and campare the results with results presented by Bai and Yeung (1974), Bai (1971) and by Drimcr, Agnon and Stiassnie (1992). Keytrords: Motion of Floating Bodies, Waves, Ship Design, Added Mass and Oampiog, Bouon • Mountcd Rccungular Cylinder.

lntroduction

ln the design of modem ships their motion behavior in waves plays an important role. ln the case of ship, the viscous effcct can be neglccted, the n uid considered incompressible, the body forces be derived from a potential function and the tlow assumed irrotational. Further, we can assume that tlle wave length are ofthe sarne order as the ship width (high frequencies), the wave amplitude small and the body considered slender. Under these condition the problem can be linearized and a strip theory developed: the body is subdivided in many scctions along the length; two-dimensional problems are solved for the motion ofthese sections and the hydrodynamic reaction forces on the body under action of incoming waves can be calculated. lncorporatin g the buoyancy and gravlty forces, the motion of the floating body can be detennined according to the strip theory.

The hydrodynamic forces are subdivided in two groups: excitation forces, resulting from the action ofthe incoming waves and their diffraction, considere<! the body fixed; and radiation forces, resulting from tlle oscillatory motion ofthe body in calm waters.

l n the linearized theory of oscillatory motion of a body close or o n the free sur face the hydrodynamic reaction is split up in two different components: one in phase with the acceleration and tlle other in phase with the velocity. The coefficient ofproportionality force in phase with acceleration is called added mass and the relation ofthe force in phase with velocity and vdocity is called darnping coefficient.

Manuacript 1'110111Ved: JanUI/y 1895 Technical Editor. Cal1ot Alberto carrasco Altemanl

J . ofthe 8raL Soe. Mechamcal Sciences- Vol. 18. Deoember 199$ 356

Tbe added mass for a body in finite domain has a coostant positive value. When the body oscillates on the f.ree surface in depends on lhe frcquency, on lhe local depth and whelher we bave a two- or three- dimensional problem. ln lhe case of a two-dimensionaJ body oscillating on the free surface in finite depth in the low frcquency limit the added mass goes to infmity, while for the case of finite deptb it goes to a finite value. The damping coefficient is related to lhe energy tlux in the outgoing waves and necessarily is always positive. For high frequency it must goes to zero.

The determination ofthe added mass and the damping coefficient for a two-dimensional cylinder has motivated research over lhe last decades. Ursell in 1949 sotved the problem of a heaving circular cylinder on the surface in deep water. He used a multipole series and a wave source poteotial. Havelock source, the rcpresent the fluid tlow motion. Later on be developed solutions for lhe sway and roll motions. ln 1953 Orim presented a solution for lhe beaving problem using also a multipole expansion. After lhat Tasai ( 1959) used conformai mapping to extend the solution to non-circular sections. Tbe use of a therc parameter Lewis transform cause some rcstrictions in the method. To avoid the conformai mapping Frank (1967) discretized lhe body contour and used sources distribution with constant intensity to rcpresent lhe velocity potential.

Extending lhe approach used by Ursell (I 949), Yu and UrseU ( 1961) devetoped a solution for the case of finite deplh. ln evaluating the added mass for low frequency their results showed a tendency of growth to infinity.

Bai and Yeung (1974) using two numerical procedures, one based on a variational method (finite­element and eigenfunction expansion) and another using lhe Oreen's second identity with a logarithmic fundamental solution, obtained finite limit for added mass in low frequency too. Bai ( 1997) confirmed these rcsults using also a numerical procedure. Yeung ( 1975) developed a hybrid model combining integral equation and eigenfunction expansion and determined the hydrodynamic forces for two-dimensional bodies near or in lhe free surface. This behavior was obtained analytically by Ursell (1976) rcviewing the earlier paper with Yu. For the low frequency problem lhe boundary conditions are of lhe Neumann kind, that is, the boundary conditions contain only derivativo of lhe potential and, consequently, lhe potential has a unkoown constant. wbich was the center of a long polemic about the problem.

Continuing the investigation Sayer and Ursell ( 1976) determined the behavior of the derivative of lhe added mass in relation to the wave number times the radius ofthe cylinder. ln the discussion ofthe paper Yeung and Newman (1976) prcsented a procedurc based on the method ofmatched asymptotic expansions to determine the constant and then the value of tbe added mass for low frequency. They proved that the eonstanl does not depend on the section form and confirm that is value is the sarne obtained by Ursell (1976) for a circle. Similar results were obtained by Keil (1974), extending the multipole cxpansion used by Orim ( 1953 ).

Mclver an Linton (1991) following a procedure similar to tbe one used by Yeung ( 1976) ana1yzed lhe differences in the added masses for a bonom mounted cylinder and a tloating cylinder for low frequencies. Sioce lhe boundary value problems are the sarne, the difference in added mass corresponds to ditferences in the constants.

For a horizontal circular submerged cylinder close to the free surface Ogilvie ( 1963) concluded that the added mass assumes negative values for some frequencies. This was supported by experimental results obtained by Chung (1997). Newman, Sortland and Vinje (1984) used lhe matched asymptotic expansion method to calculate the added mass and damping coefficient for a submerged two-dimensional vertically infinite rectangle executing beave motion.

Yeung ( 1981) solved the problem of a vertical truncated circular cylinder in finite depth using an eigenfunction series solution. This procedure allows us to take into account more precisely the behaviour ofthe added mass for low frequencies. Yeung and Sphaier (1989) showed lhat lhis limit is finite when the cylinder oscillates in a cbannel. Mclver and Evans ( 1984) used the eigenfunction method for the case of a bottom mounted vertical circular cylinder and determine negative values for lhe added mass in some frequency range.

ln tbe present work we apply lhe sarne procedurc used by Yeung (1981) and Yeung and Sphaier (1989) to s tudy the hydrodynarnic properties for a bottom-mounted two-dimensional rectangular cylinder using series of eigenfunction to reprcsent the velocity potential. A complete linear solution is

357 P. T. T. Eaperança et el.: Added Mess and Oemplng on Bottom-Mounted ...

obtained and lhe behavior of added mass and damping is presented. Further an asymptotic solution for high frequencies and a shallow-water solution for low frequency are developed. For small clearance, between the of the cyl inder and the calm water levei, we determined the range of frequency for negative added mass. We also h ave applied lhis melhodology for lhe case of a Ooating rectangle and evaluated the added mass and damping. We obtained good agrcement with lhe results presented by Bai (1977) and Bai and Yeung (1974), and presented by Drimer, Agnon and Stiassnie (1992), which used a simplified model based on the eigenfunction expansion to analyze lhe behavior of a rectangular breakwater.

The Boundary Value Problem

The water deplh is detined as h •• the width of lhe bottom-mounted rectangle is equal to 2a an the clearance between lhe top oflhe rectangle and lhe free surface at lhe rest position is hj. The rectangle height is d = he- hj. To describe lhe fluid flow we use a cartesian system of coordinates Oxz. The origin is located on the mid point of lhe top of rectangle and the axis Oz points vertically upwards.

We assume that the viscous effect can be neglected, the fluid considered incornpressible, lhe body forces derived from a potential function and the flow is irrotational. Under lhese assumptions, lhe velocity can be derived from a potential v = '17'1', and the potential function satisfies Laplace equation, '172 '1' =O. Function 'i' is eOxpressed in the form:

\11 (x, z, t) = 9l [ - iat; C1> (x, z) exp ( - iat)] ( I)

where a rep!esents the Jrequency of the oscillation, t; the complex amplitude of lhe vertical motion t; = t; exp ( iõ) , t; the module of the amplitude, õ the please angle and i = H · 9l means lhe real part.

Function C1> must satisfy the Laplace equation over lhe fluid domain

and lhe following boundary conditions:

2 ôCl> (x, z) = a C1> (x, z)

õz g z =h.

I

ôCl> (x, z) az

ôCl> (x, z) = I ôz

ôCl> (x, z) = 0 ax

z = O;lxl S a

z = -d;lxl ~a

-d S z S O;lxl = a

ôct>(x,z) 'f imeCJ>(x z) -+0 õx o ,

Internal and Externa! Boundary Value Problems

(2)

(3)

(4)

(5)

(6)

(7)

The fluid domain is divided in two parts, an internal and an externa! rcgion as showed in Fig. I. The internal domain is a fl uid rectangle located over the top ofthe cylinder wilh lhe sarne width as the body, 2a, and height equal to hj. {-as x s a; O~ z ~ hj} . The externa! domain is located at lhe side o.f the body and the internal domain so lhat {x ~a; -d $ z s hj}. We introduce lhe potential functions Cl>1

and Cl>e , so that C1> = eJ>i in lhe internal region and Cll = CJ)t in lhe externa! region.

J. ofthe Braz. Soe. Me<:hanieal Scíen08$- Vol. 18. December 1996 358

free surface Z-fl

t z h; -- ------ ---- ---i

top z= o

h. 2a

d

bottom z- -d

Fig. 1 Bottom-fflounted Rectangular Sectlon

cJ>i satisfies equation (2}, the boundary conditions (3) and (4). Further, the fluid motion is symmetric about the Oz axis, then we can study the fluid Oow in a half part of the domain and impose õcJ>i("'z) / ôx = O forO :S z s hp=O.

cJ>C satisfies equation (2}, the boundary conditions (3), (5), (6) and (7).

Both boundary value problem can not uniquely determined, since we did not defined boundary conditions for OS z S h;, Jxl = a We introduce complementary boundary condition representing lhe compatibility between the two potential given by the continuity on lhe pressure:

and continuity of lhe horizontal velocities:

o for -<1 S zs O

õcJ>i ("' z) l ôx for Os z S h;

Eigenfunction Expansions

x= a

x= a

(8)

(9)

The internal solution, satisfying (2), (3), (4) and the symmetry condition, can be written as lhe sum of a particular C!>~ and a bomogeneous solution C!>~ :

..... ; = ..... ; + .... ; 2 h;b0 cos(mbx)Zb(z) ~ h;b0 cosh(m~x)Z~(z) ~ ~P ,.,h = z- h;+ g / a + . + L. .

cos (m0a) n . 1 cosh (m~a)

(lO)

where bn are unknown complex coefticients, b = bR + i b I , z i the eigenfunctions, m0 and mi are th e eigenvalues. mb is the first eigenv~lue, "the s'blu,ion of the dispers ion equa~ion a2f g- mtanh (mh.) = Q. mi are the others eigenvalues, solutions of a2/ g + mtanh (mh.) = o. The eigenfunctions áre: "

1

Zb (z) cosh (mbz)

Po and z~ (z)

cos (m~z) for n ;:::: I (I I)

l59 P. T. T Esperança et ai.; Added Mass and Damping on Bottom-Mounted ...

with

and 1[ sin(2m~hi)J Ni = - I + for n 2:: 1 " 2 2m i h.

n •

(12)

Similarly, the extemal solution is given by:

d0exp(im8x)Z8(z) ~ d1 exp(-mfx)~c(z) <l>e (X. z) = + L. me2 exp (imea) me2exp ( - mea)

O O I . I I I

(13)

. R . I e . . e e where d, are unlcnown fOIDplex coeflic1ents, d0 =do + ld0 , Z0 the elgenfunCtiOnS, IDO and m0 are lhe eigen values . m0 is tpe first eigen va lue, t he solution of tb.e d ispersion equation a 2 / g - m tanh (mhe)= O. ID0 are the other eigenvalues, solutions of <f / g + mtan(mh.) = O. Tbe

eigenfunctions are:

e _ cosh [ m~ (z+ d) J Zo (z) - X and ~ (z) =

cos [ m~ (z + d) J R for/ 2:: I (14)

with

e 1 [ sinh (2m~hJ ] • ~ _ 1 [ 1 sin (2m~hc) ] N z - 1 + and N 1 - - + -----'-~ O 2 2 eh 2 2meh

mo c I e

for I~ I (IS)

The eigenvalues m~and m~, related to the dispersion equation, indicate the presence of progressive waves in the internal and externa! domains.

Determination of the Coefficients

The deterrnination ofthe coefficients ofthe internal and the externa! solutions can be accomplished by, applying the compatibility conditions (8) and (9) and utilizing lhe orthogonality property ofthe z; and Z~ functions. Let us define the operators:

c.t = h ~h t1g (x,z)~ (z)dz (16) 2 I I

where the superscript j is equal to i o r e and the superior limit the integral h2 is equal to h; and the inferior limit is O or -d depending on whether we are dealing with the internal or the externa.! eigenfunction z; or z: respectively.

and

Applying C.~ and Ç~ respectively to equations (8) and (9) and using (lO) and (13) we obtain: .. . fk + bk = L c~, d1 ,_o

00

qjdj = L cj: bm m- o

(17)

(18)

where f.k ' c~,'. qj and cj~ are the results of the application of ç~ and ç~ to the potential functions <t>\ <I>~ and <I>~ . Substituting ( 17) in (18) we can solve the system for b and the recover d.

J. of the er.z. Soe. Mechanical ScienCoe$- Vol. 18, Oecember 19G& 360

Hydrodynamic Coefficients and Wave Profile

The hydrodynamic force acting oo the body is given by the integral ofthe hydrodynamic pressure p = -qà'f'1 1 at obtained from lhe linearized Euler's integral on the top ofthe rectangle:

F(t) = fp0:1

ds = -11Ú(t) -Ã.U(t) a -(11+i~)Ú(t) {19)

where the body vertical velocity is given by U (t) = 9t [ -ia~exp (iô- iat) J "" a~sin (ô- at) , Ú is the accoleration and IA and I. are the added mass and the damping coefficient respectively. lntroducing equation (tO) in ( 19) a straightforward integration leads to the following expressions for the hydrodynamic complex coefficient:

11 + •- .. 2ph.a - L . + I - -"'- ( .À, [b0tan(m~a)Z~(O) ~b0tanh(m!a)Z!(O) ( o)] 20) (J I I I 2

m0 a n - 1 m0

a a h.

The free surface elevation n (x, t) is also obtaincd ftom lhe linearized Euler's integral:

{ 1 a~(x h. t)} ~ . .

'lj(x,t) = 9l -g a~ 1' =

0

8 [<~>~ (x,h1}cos(S-at)-~ (x,h;)sin(ô - at)J

(21)

where <%>~ and <1>'1 the real and imaginary parts of the complex potential ~ = <%>~ + i ~1 and the superscript j depends on whether we are considering the internal domain, j =i and lxl ~ aor the ext.emal domain, j = e and lxl ~ a. The real and lhe imaginary potential in these domains are givcn by:

,.,i 2 h;b~ cos(m~ x)Z~ (z) .ç. h;b~ cosh(m! x)Z! (z) (22) ....-R • z - h. +g/a + . + .t- .

1 cos (m~ a) o • 1 cosh (m~ a)

I i i I i i <ti = h1h0 cos(m0 x)Z0 (z) + f. h;b0 cosh(m0 x)Z0 (z) (23)

1 cos (m~ a) n & 1 cosh (m~ a)

<O' { d: "'' [ m; (x-a) ] - d~ ,;" [ m; (> - a) J) Z,' (') + f. d~ oxp (-m~ >) z; (') R • 2 L.. e 2 e

<I> e I

(m~ ) I ~ 1 (m1 ) exp ( - m1 a)

(24)

{ d~ cos( [ m~ (x-a)]+ d: sin [ m~ (x-a)])} Z0e (z) "" d: exp (-m~ x)Z~ (z) ------------:2:---------- + L e 2 e

(m~) 1• 1 (m1 ) exp(- m1 a)

(25)

Behavior ln Shallow Waters

To analyze the behavior fQ_rjsmall clearances and water depths we use the shallow-water approach. The mean potential function <%> (x) is defined as:

381 P. T. T. Esperança et 81.: Added Mass and Damplng on Bottom-Mounted ...

q,i (x) = h ~h fl4} (x, z) dz 2 I ~~ (26)

where the index j means internal or external domain, the inferior limit in the integral h1 is O or -d, the superior limit is h; e h. respectively and <i>J (x) is the potential in shallow waters.

lntegration of Laplace equ~tion in the ~Srtical coordinate, application of the boundary conditions and use ofthe approximation CllJ (x, h i) ,_ «ll (x) yields to

82cr,i (x) + 0 2cr,i (x) 2 h ~L~ ~n

ôx g j

where r.s. v. means the right si de value depending on index j. For j = i the equation refers to the internal domain and r.s.v. is equal to llh;. while for j = e the domain is the externa! and r.s.v. is equal to O.

From the dispersion equation we obtain for shallow water

( . c)2h ("i)2 21 mo e = mo = 0 g (28)

Equation (27) conducts to the solutions:

-. -i 2 «ll' (x) == b0 cos (m0 x) + g/ a (29)

(30)

ln the internal region we have a stationary wave b0 cos (m~ x) , the homogeneous solution. The particular solution is a constant representing the motion ofthe fluidas a rigid body. For the externa! domain we have a progressive wave. The complex coefficients b0 and do are to be determioed using equations (8) and (9)

by

-i . -i .·e 2 m0 hism (m0 a) exp ( - 1m0 a) g/a

do = im~ becos (m~ a) + m~ hisin (m~ a)

(31)

Accordiog to lhe shallow-water approximatioo the vertical force F (t) acting on the body is given

and

2 jll - i F(t) • - pa \Jt{2Çexp(-iat)j «ll (x)dx}

o (32)

lntroducing the solution for the internal potential in shallow waters, equation (29), we obtain for the hydrodynamic coefficients:

(33)

J . ofthe 8raL Soe. Mechanical Sclenoea · VoA. 18. Deoember 1996 362

2pg (m0 ehe) (m~ h e) sin 2 (m~ a) À "" ------=---=--___,::....:: __ __:: ___ _

. . . 2 . 2 m~ o { [ Õl~ hisin (m~ a)] + [ m~ becos (m~ a) ] }

(34)

with the limit values for m~ a ~ o

• 2 a (he ) a "' -- - - 3 3he h;

(35)

where ll = m~ /m~ = {h/ he) 112

in (33) and (34) is lhe "clearanceness" parameter.

lntroduc ing (29) in (21 ) we obtain the expression of the free surface elevation in the internal dornain:

2 i -i o Çcos (m0 x) (36) 11 (x, t) = [bRcos(/5 - ot) - b1sin(li - ot)) + ~cos(S - ot)

g

We observe that for m~ a = nn , with n ~ I the imaginary part of the coefficients b0 and d0 vanish,

(37)

and the wave eontribuition of the added mass and the damping coeflicient are ~· The added mass is negative and given just by lhe motion oflhe tluid as a rigid body 1..1 = - 2pga/ ~ . Ali the properties depending on lhe est.emal potential vanish, which means that the fluid in lhe extemal domain is at rest On lhe olher hand a stationary oscillation is present in lhe internal domain wilh potential, horizontal velocity and wave pro file given by:

- i 2[ n + 1 i J «<> (x) - g / o I+ (-1) cos(m0 .x) (38)

- i ael> (x) - i b . (-i ) ax = - m0 R SIO m0 X

(39)

(40)

The wave e levation has a behavior similar to a resonant mode: there is a standing wave in the internal domain, while outside there is no mot.ion. For x .. a the horizontal flux and the elevation ofthe free surface are always zero. At this point the standing wave has a vertical motion (crest/through) but lhe surface elevation remains equal zero since lhe fluid motion as a rigid body compensate the wave CO)llponent and lhe wave lengh L; • 7nlm~ • 2a/n, is a fraction of lhe rectangle width, since .... . . m0 a - n7t where n IS an mteger :2:: 1.

Now we observe the behavior when m~ a • (2n - I) n/2 with n ~ I . We can identify another resonant mode for this case. The internal and lhe externa! potential functions and the horizontal fluid velocities are:

363 P. T. T. Esperança et ai.: Adc:led Mass anel Damping on Bottom-Mounted ...

and

<i>c(x) = ~{cos[m~ (x - a)]+isin[m~ (x - a)]l (J

-c -e ô<ll (x) mo g ( . [.c --- = -- -sm m ax (J2 o

(x-a)]+icos[m: (x-a)]l (42)

For x =a

-i - i -e <ll (a) = if>e (a) = g / o2 and ô<ll (a) = igm i I (to2) .. ô<ll (a) = . -e I 2

ôx o ôx •&mo o

Now the wave profile in the internal domain is given by

lii(x,t) = (- l)n + ll_cos(m~x)sin(ô-ot)+~cos(õ - ot) E

(43)

(44)

and the wave length L1 = 2a [ (2n- I) 12) with n ~ I , that rectangle width is equal an integer

number of wave lengths plus half wave length with wave nodes located at lxl = a. The nodes move

vertically because of the tluid motion as a rigid body. There is a time oscillatory tlux under the nodes.

Reduction o f lhe clearance, for a constant water depth, implies in an increase of lhe wave amplitude.

for lhese frequencies lhe wave component of the added mass also vanishes, but the darnping is

different from zero. Writing lhe darnping coefficient, Eq. (34), as a funcrion ofthe product m~ a is:

't • i ~o.(m0 a)=

• i . 2 • I 2pgathihcmo asm (m0 a)

. . 2 . 2 a { [ rr.~ ah1sin <m~ a> J + [ rr.~ aMcos <m~ a> J }

(45)

The analysis of its derivative wíth respect to m~ a índicates lhat 1.. assume local maximum values at this poínts.

The High Frequency Problem

lntroduction of the high frequency assumption in boundary condition on the free surface (3) yields to

i <ll(x,z) =O for z • h

I (46)

and consequently a modification in the radiation condition, not prescribing lhe presencc of djverging waves.

The procedure to obtain the asymptotic solution of the problem is the sarne used to reach the complete solution

Thc potcntial in thc internal domain is given by:

i ~ h.bncosh (m~ x) z~ (z) <ll..,(x,z) = ~ 1

1 + (z-h;) n •O eosh(mna)

where lhe eigenfunctions. lhe nonnalizing factors and lhe cigenvalues are given by:

i cos (m

0 z) .

z~ <~ = Ni R n

(2n + 1)1t

2h. I

with

(47)

(48)

J. of lhe Braz. Soe. Med!ank:al Sciences - Vol. 18, Deoembef 1996

ln the externa! domain the eigenfunction expansion is given by:

e L .. dlcZ~(z) exp(-rn: x) <I> (x,z)=-

ao e 2 e 1c • o mie exp ( - mk a)

(49)

whcre the eigenfunctions, the normalizing factors and lhe eigenvalues are given by:

z: (z) = cos(m~z+d). e . e_ {2k+l)1t . X • Nk = 1/2• mk - 2 h0

wtth (SO)

Two differences frorn lhe complete solutions (I O) and (13) are to be noted. First., the particular solution is modified, and lhe progressive mode does not appear. Consequently, lhe damping coefficient is zero.

Solving the problem and evaluating the hydrodynamic force we obtain the expression of the nondimensional added mass:

R . ~ bntanh (m: a) Ji

11.., = -2phi L i + 2pahi n - o m

D

(S I)

Rectangular Cylinder Piercing the Free Surface

We consider now lhe case oflhe section heaving on the free surface. The x-axis coincides wilh lhe free surface at rest position. h; = O and d is lhe rectangle draft. The solution is giveo by

2 2 i <I>i = <I>i + <I>i = z + 2hez-x + f b cosh (p0 x)Z0 (z)

P h 2 (he -d) 11

_ 0

° cosh (p0a)

(52)

where b. are unknown coefficients of the homogenous solution, Z~ thc eigenfunctions and Pn are the eigenvalues. Applying a similar procedure as presented above we obtain the added mass and damping coefficien t:

ll +i!:: = 2p(boa+ f A(- l)"tanh(pka)) +2p(d2 -2heda - aJ ) C1 k - 1 Pk 2 (h

0- d) 6 (h

0- d)

(53)

Results and Conclusion

Using the methodology developed here we determined the added mass and damping coefficient for a bottom mounted rectangular cylinder. For ali calculations we used 20 internal tc rms of the eigenfunction expansion and 60 externa I tenns. First we observe the effect of the clearance on the hydrodynamic coefficients for fixed rectangle width. Fi~ures 2 to 5 present nondimensional added mass ( 11/ pal) and damping coefficient (1..1 pa2 {g/a) 12 )as functions of the nondimensional parameter m~ athe internal wave number times halfwidth. The relation halfwidth to water deptb is given by a/he = 0.10 and 0.20. Results are presented for values ofthc n:lation clearance to water depth h/h. = O.OS, 0. 10, and 0.50, using the complete solution and the solution for high frequency. From these results we can say that: I) the added mass and the dimensional damping coefficient assume tini te values for low frequency ; 2) for large clearance, small va.lues of alh;, the added mass is positive for ali frequencies. Negative values appear for small clearance, large values of alh;. The frequency range of negative added mass increases by increasing alh;; 3) results for the added mass from the complete solution converge to the asymptotic solution for high frequency (see also Table I).

365

-1.5

P. T. T. Eape111nça et ai.: Added Man and Camping on Bottom-Mounted ...

- N/he • 0.05 a~ • 2.0 ~N'/he • Q.IO o~ • 1,0 <JBBBE) N';t. - 0.50 0/N - 0.2 -~ ...... ~~-.. OBI!BO~ . ....._ a;t. • 0 .10

-~5~~~~nm~~~~~~~~~~Trl 0.0 3.0 ' · 5.0

mo'a Fig~ 2 Nondlmenalonal Added M.ua Ve111ua Internal Wave Numbe

Lo~~~~~~~~~~nT~~~~n,

5.5

l..O

-2.0

- hilhe - 0.05 M6M hfl.he - o. 10 Q99BD híjh. - 0..50 -...,.np ......... 6MM...,.np. ........ CJIIB9E) ...,.np . .......

e,,h. - 0.20

o/,!'j - 4.0 0/.fM - 2.0 oi'hl - 0.4

-7.0 -h.,,....,., ....... "'"""""" ............. r?Y'I,...,.,...,........,......, ...... ,...,....-1 o. l

m..'a Fig. 3 Nondlmenalonal Added ..... Vel"llua lntamal Wave Number

[l.O Cl

';;a ;;; s2.o ~

1.0

- 111/he • 0.05 o/.hl • 2.0 666Mhl'/he • 0.10 o/hl • 1.0 CJIIB9E) h(;t. - 0.50 o/hl - o .2

o/he- 0 .10

Fig. 4 Nondlmenalonal Camping Ve111ua lntamal Wave Number

J. oflhe Braz. Soe. Medlank:al Sdences- Vot. 18. Oecember 1996

~ - 0 .10 ct/hl- 1.0 - ~~ he - o.~ 11/.h! - 4.0

OIIJEI) - 0.50 a/hl - 0.2

olhe - 0.20

Fig. 5 Nondlmenalonal Oamplflill Versua lnt~ Wavw Number

Figures 6 and 7 show results of added mass and damping coefficient for a case with a much smaller clearance, h/he = 0.01 and alh. = 0.10. We also present results obtained with Newman, Sortland and Vinje (1984) solution for a rectangular cylinder with the sarne widtb and infinite deptb, and calculation using Frank.'s close-ftt method for a submerged square cylinder in finite water depth. The frequency pararneter and the non-dimensional forms are the sarne used by Newman, Sortland and Vinje. Using the eigenfunction expansion and Newman, Sortland and Vinje solutions the hydrodynarnic coefficients are determined by integrating the pressure only on the top of the cylinder while in Frank's method bottom and top side are considered. Discrepancies appear for low frequencies due to the different approaches. lt is to note that an infmitc limit for the addcd mass is expectcd in the solution presented by Newman et ai.

1Ul

10.0 - [Joenfunctlon

a. o -~ OIIIH)f ......

1.0

•• 4.0 '$:- 0.01 he - 0.10

~ z.o a/111 - 10.0

::t o.o

- 2.0

-4.0

- 1.0

- 1.0

· 10.0 0.0 0.4 O ..I 1.2 1.1 2.0

a (a/a)'11

Fig. I Nondlmenalonel Addld Maaa Veraua Fraquancy

367

•

P T. T. Esperança et ai. : Added Mass and Camping on Bottom-Mounted ...

20.0 ....,...,...,...,...,,..,....,,..,....,...,...,...,...,...,...,...,...,..,...,..,...,..,...,..,.....,..,.....,..,.....,...,...,....,...,....,...,....,...,...,...,.,

18.0

16.0

14.0

--- Eigenfunction ~Newman GEa:3a Frank

~ 12.0

Q.. 10.0 hi/he • 0.01 o/he • 0.10 o/hi • 10.0 :i 8.0

8.0

4.0

2.0

o.o -&r""T"T"'f""f;......T"T"1""T"l,..,..;:;::;:=-;;::,.,.,..I"T""nEil~llr"r'.,.,..,.,...-..n-t 0.0 0.4 0.8

cr(a/g)1/2 1.2

Fig. 7 Nondlmenelonal Camping Versus Frequency

1.6

ln Figures 8 and 9 we present the influence of lhe rectangle width in lhe bydrodynamic coefficients for a fixed clearance. The clearance is 5% of tbe water depth h;lh,. .. 0.05, and the values for the relation width to water depth are a!h,. .. 0.01, 0.5, 1.0 and 2.0 corresponding to width to clearance a!h; = 0.20, 10.0, 20.0 and 40.0. These leads to the sarne conclusions presented above.

25.0 ..,...,.....,..,,.....,..,...,...T""r.....,....,..,....,...,...,...,.....,..,,.....,..r-r,...,.."Tõ....,..,...,r-r-,..,

15.0

.. IS

5.0

~ ::t -5.0

-15.0

-25.0

..... a/he ~a/he a:Hfla/he ......-a/h•

0.05 a/hi • - 0.25 o/hl • • 0.50 o/hi

1.00 a/hi •

1.00 5.00 10.0 20.0

hf/he • 0.05

-J5.0 -h-....,...,r-rT...+rrr-rrrr-r-rr...,...,"TTTT'rTT'TT"rrT"TTT"T"1r-r! 0.0 4-.0 8.0 8.0

"In4•a

Fig. 8 Nondlmenslonal Added Masa Versus Internal Wave Number

J. of lhe Braz.. Soe. Mechanical Sciences- Voi. 18, December 1998

.-..15.0 ~ --., ............ «30 -• ., 10.0

s ............ ;<

5.0

- o/he • 0.05 a/hl • bbl:lbiJ. oZ'n• - 0.25 a"/hi • a:a:E> oZ'n• - 0.50 a"/hl )8888( a/he - 1.00 afhl •

'ni/h• - 0.05

2.0 mo1a

4.0

1.00 5.00 10.0 20.0

Fig. 9 Nondlmenslonal Oamplng Veraus Internal Wave Number

388

e.o

ln Figures lO and II we compare the complete and the shallow-water solutions. ln this case h/h.= 0.5 and alh. = 5.0 corresponding to a/hi = 100.0, a very small clearance case. The agreement is very good for ali frequencies.

40.0

-10.0

-eo.o

-110.0

hijh. - 0.05 al)t• • 5.0 a"fhl • 100.0

-lfiO.O +r"'I"TTT'I'"T"''T-r-TT'...,....'T"T'TT'IrTTTT"rT"'..,...,.TT"rT"IrrTTT"rl

0.0 2.0 l 4.0 6.0 8.0 moa

Fig. 10 Nondlmenalonal Hldrodynamlc Coefflclenta

3&9

-1~.0

P. T. T. El~nça et ai.: Added Mass and Damping on Bottom-Mounted ...

hl./he - 0 .05 .,__~- 1.0 fi/N - 20..0

Fig. 11 Nondlmenalonal Hldrodynamlc Coefflclents

140.0

40.0

hl/loe • O.DI!I ,.,.. - ~.00 e;N - 1Cl0.0

m.'a a.o

Fig. 12 Nondlmenalonal Hldrodynamlc Coefftclenta ln Shallow Water

Two contributions detennine lhe added mass (20). One part dueto the tluid motion as a rigid body in the upper part of lhe rectangle. Thc other part h as an undulatory character due to lhe waves. ln Fig. 12 lhe total added mass, lhe contribution of lhe undulatory tenn to the added mass and tbe da_mping coefficient are plotted in nondimensional fonn as function of the frequency parameter m~ a for shallow water. The contribution due to lhe wave motion is zero when lhe fn:quency parameter is equal to n / 2, and simultaneously the damping coefficient is maximum. This corresponds to the second kind of resonance presented in the behavior in shallow water section. A standing wave appears o ver the rectangle wilh length equal to two times the rectangle width. The modes are located on x = ±a , originat ing an !>Scillatory tlux from tbe internal to the externa! regions. A second mode can be observed for m~ a - 3~t/2 . ln this case lhe wave length is equal to 4a/3.

J . ol lhe &az.. Soe. Mechanical Sdeoces · Vol. 18, Decembef 11196 370

V alues of added mass for different parameters h;fhe and alhe obtaíned from the complete solution (EIGEN) were compared with values for lhe high frequency approacb (HASYMP), and values published by Mclver and Linton (1991) for low frequency limit (LINTON) . Table I contains these results. By comparison the agreement can be considered very good.

Table 1 Added Maas for Hlgh and Low Frequencies

ald h/h. alh. IJr IJ IJoo IJ Linton Eigen Hasymp Eigen

0,5 0,333 0,333 1,398 1,403 1,398 1,389

0,5 0 ,500 0,250 1,102 1,109 2,049 2,045

0,5 0,666 0,167 1,29-4 1,296 2,732 2,727

1,0 0,500 0,500 0,626 0,629 1,406 1,397

1,0 0,666 0,333 0,704 0,710 2,083 2,078

1,0 0,875 0,125 1,722 1,739 3,452 3,448

2,0 0,333 1,333 1,080 1,085 0,460 0,457

2,0 0,500 1,000 0.084 0,085 0,844 0,830

2,0 0,875 0,250 1,030 1,035 2,750 . 2 .743

For a flxed relation ald we observe from Table I that the added mass for mi a-+ o decreases and then increases for increasing h;~- Further the limit value of the added mass for the shallow-water solution assumes negative values when h; > hJ3 (Eq. 35). This limit does no! means necessarily shallow waters in the internal domain since we are dealing with three parameters: m~ a, m~ 111 and m~ h e .Since the externa! do~ain is infinite shallowness is defmed by small m~ "e . ln the domaín it does not depend only on m~ h i but also on how large alb; is. Figure 13 prese.nts results of the nondimensional added mass and damping coefficient for a very small value of m~ a as function of h/he and alh;, using the complete solution. Since lhe values of alh; are small, lhe results for added mass do not agree with the limit value obtained ln section on behavior in shallow waters. But we have a exact agreement for the damping coefficient. Table 2 presents results for added mass obta.ined from (20) and (33). For large values ofalh; • 10, 100 and 1000, the differences are about 50, 5 and .5% respectively, confuming what was expected.

12.0

10.0 •ct

~ 11.0 rn.' a • 0.001 -~ ... 11.0 -Cl

.......... 4.0 III -• Cl .s 2.0

.......... .<

Fig. 13 Nondlmenalonal added mau aa functlon of a/1\1

371 P. T. T. Esperança et ai.: Added Mass anel Damplng on Bottom-Mounted ...

Table 2 Added Mass for Shallow Waters for ~a =0.001

h/h. .3 .3 .4 .4 alh, Compl. Shallow Compl. Shallow

0.001 8.74 0 .00022 10.75 -0.00033 0.01 3.43 0.0022 3.49 -00033 0.1 2.05 0.022 3.49 -0.0033 1.0 1.45 0.22 0.61 -0.33 10 3.32 2.22 -2.57 -3.33

100 23.30 22.22 -32.59 -33.33 1000 223.3 222.22 -332.59 -333.33

ln Figures 14 and 15 we present resulls for added mass p • .. /(p2ad)and damping coetlicient :t • 'J... I ( pa2ad) for a rectangular cylinder, with half width to draft = 1.0 and water depth I draft = 1.5, 2.0, 3.0, 4.0 and 5.0, heaving on the free surface. Results for these cases are presented by Bai and Yeung ( 1974) and Bai ( 1977). Although we do not plot their results here, we can say the agreement is very good except for some differences for low frequencies.

2.50 ~~ ............ ~ ....... ....,..,..;..... ............... ....,..,.,.,.,,..,.,..,..,,.,.n-TI

2.25

2.00

"01.75 o Q.. ~1.50 ::t

1..2.5

1.00

0.75

<B38E) h /. d - 1.5 IBB:lh z d - 2.0 ~h Z d • J.O ~h z d- 4.0 - h"/ d- 5.0

Fig. 14 Added Mus for a Rectang,. on th• Fr .. Surfaee

3.0

2.5

1.0

0.5

<388E!Oh /. d - 1.5 IBB:lh z d - 2.0 ~h Z d • J .O ~h z d- 4.0 - h"/ d- s.o

Fig. 15 Camping Coemclent for 1 Rectangle on the Free Surfaçe

373 P. T. T. Espeornça et al: Added Mass anel Damplng on Bottom-Mounted ...

Yeung, R. W., 1981, "Added Mass and Damping of a Vertical Cytinder in Finite-depth Waters". Applied Ocean Rcsearch, 3:119-133.

Yeung, R. W. and Newman, J . N ., 1976, wln the Discusion of Ursell ( 1976)". ln li th Symp. on Naval Hydrodynamics.

Yeung, R. W. and Sphaier, S. H. , 1989, "Wave-interference Effects on a Truncated Cylinder in a Channel". Joumal ofEngineering Mathematics, 23:95-1 17.

Yu, Y. S. and Ursell, F., 1961, wSurface Waves Generated by an Oscillating Circular Cylinder on Water offinite Oepth: Theory and Experiment". J. Fluid Mech., 11:529-551.

J . oftlle Elfaz. Soe. Mechanicar Sciences. vor. 18, Deoembef" 1996 372

Drimer, Agnon and Stiassnie (1992) used a simplified model based on the eigenfunction expansion to anal(ze the behavior of a rectangular breakwater. ln Fig. 16 we present the added mass ~ - Ilhe / p and the damping coefficient >. - >.h; MI p for halfwidth lo draft = 1.0 and water deplh

I draft = O. 7, in the sarne dimensional form of Drimer, Agnon and Stiassnie to be easily comparable. The results are in good agreement; in much better agreement with the numerical results showed in their pape r.

-- odcled mon --- domping

4.0

J.O

2.0

Fig. 11 Added Ma11 and Oamplng Coefftc:lent for 1 BnNkwater

References

Bai, K. J., I 977, "The Added Mass ofTwo-dimensional Cylinders Heaving in Water of Finite Depth". J. Fluid Mech., 81 :85-105.

Bai, K. J. and YeWJg, R. W., 1974, "Numerical Solutions to Free Surface ftow Problerns". ln Tenth Symposium on Naval Hydrodynamics, Cambridge, Mass ..

Chung, J. S., 1977, "Forces on Submc.rged Cyllnders Osciltating Near a Frec Surface". J. Hydronautics, 2:100-105. Drimer, N., Agnon, Y. and Stiassnie, M., 1992, "A Simplified Analytical Modcl for a floating Brealcwater in Watcr

ofFinite Depth". Applied Ocean Research, 14:33-41. Esperança, P. T. T., 1993, "Sobre o Comportamento Hidrodin&mico de Corpos Proximos a Superllcie Livre". PhD

lhes is, Programa de Eng. Ocdnica/Coppe, Universidade Federal do Rio de Janeiro. Frank, W., 1967, "On the Oscitlation of Cylinders in or Below thc Free Surface of Deep fluids". Report 2375,

Naval Ship Res. & Dev. Ccntcr, Betbesda, MD. Grim, 0., 1953, "Berechnung Der Durcb Schwingungen Eines SchiffslcOrpers Erzeugten Hydrodynamischcn

Krlfte~. Jb. Schiffbauteeh. Ges., 47:277·296. Keil H., 1974, ''On the Oscillltion of C)'linders in or Below the Free Surface of Decp Fluids". Report 305, lnstitut

fiir Schiffibau Hamburg, Hamburg, Oermany. Mclver, P. and Evans, D. V .• 1984, "The Ocurrence ofNegativc Added Mass ln Frce Surfacc Problems lnvolving

Submerged Oscillating Bodies". J. Eng. Math., 18:7-22. Mel ver, P. and Linton, C. M., 1991, "The Added Mass ofBodies Heaving at Low Frequcncy in Watcr of Finite

Depth". Applied Ocean Research, 13:12-17. Newman, J. N., Sortland, B. and Vinje, T., 1984, •Added Mass and Dampíng of Rectangular Bodics Close to the

Free Surfacc". J. Ship Res., 28:219·225. Ogilvie, T. F., 1963, "Firs1- and Second-order forces on a Cylinder Submerged Under a Free Surfacc". J. Fluid

Meeh., 16:451-472. Sayer, P. and Ursell, F., 1976, ·0n the Virtual Mass, at Long Wavelength.s, of a Half-immersed Circular Cylinder

Heaving on Wateroffinite Deptbs".ln IIth Symposium on Naval Hydrodynamics, London, England. Tasai, F., 1959, "On tbe Damping Force and Added Mass ofShips Heaving and Pitching" . J. Zosen Kiolcal , 105:47-

56. Ursell, F~ 1949, "On the Heaving Motion of a Circular Cylinder on thc Surface of a Fluid". Quart. Joum. Mech.

and Applied Math., 2:218-231. Ursell, F., 1976, "'n lhe Virtual-mass and Damping Coefficicnts for Long-waves in Watcr offinite Depth". Fluid

Mech., 76:17-28. Yeung, R. W., 1975, "A Hibrid lntegral-equation Method for Time-harmonic Frce Surface Flows". ln 1st

lntcmational Conferencc on Nwnerical Ship Hydrodynamics.

RBCM • J. of tha Bru. Soe:. Meehanlcal Sc:lanc:.s Vol. XVIII - No.4 • 111!14 - pp. 374-382

ISSN 0100-7386 Prlnted in Brazil

Vibration Analysis of Structures Subjected to Boundary Condition Modifications Using Experimental Data Domingos Alves Rade Universidade Federal da Ubelllndia Depto. de Ciência$ Flslcas • C.P. 593 38<400-902 Ubelllndla, MG Brasil

Gérard LaJiement Univerlité de Franch&-Comté laboratolre de Mécanlque Appllquée R. Chaléat 25030 Besançon, France

Abstract This paper addresses the Analysis of Modified Stroctures by using experimental data. ln particular, modifications ofthe boundary conditions ofthe s~ture by grounding ofooe or severa! ofits degrees of freedom are considered. Three methods are proposed, wbicb are conceived to obtain the eigenvalues and eigenvectors of more constrained configurations, giveo tbe experimental Frequency Rcsponse Fuoctions measured on a less conslrained configuration. The formulations oftbe thrce methods are first presented and tbeir performances are then evaluated through applications to an automotive structure tested io laboratory. Keyworcls: Structural Modifications, Modal Analysis, Antiresonances

lntroduction The fundamental problem of the Analysis of Modified Structures (AMS) consists in, given the

dynamic characteristics of a structure ln a reference configuration, named herein the initial configuration (IC), to determine the dynamic behavior ofthe so called modified configurations (MCs), which are obtained by introducing in the IC, known modifications of its physical properties, i. e., mass and/or stiffness and/or darnping.

Severa.! techniques of AMS, based either on the use of theoretical models, like finite element models, or on the use of experimental data, were developed and became wide spread the last two decades (Soyder, 1986). The latter c lass oftechniques presents, with respect to lhe fonner, the advantage of being applicable to complex structures for wbich ao accurate fioite e lement modeJ is difficult to be obtained.

ln the present paper, some techniques of AMS, based oo the exploitation of experimentally measured frequency respoose functions (FRFs) are examined. The structural modificatioos considered are those representing modifications ofthe boundary conditlons ofthe IC by grounding of one o r several of its degrees offreedom (DOFs). ·

The study reported herein was mainJy motivated by the fact that physicaJ grounding of DOFs is, as a rule, ooe ofthe major drawbacks arising when perfonning vibration tests, sioce the mechanical devices designed to inhibit the motion are oever ideally rigid for large frequency bands. Tbis f&ct can entail serious difficuJties. As an example, if the experimental data are to be used for the adjustment of theoretical models, the true flexibilities ofthe grounding devices, which are in general dlfficult to estimate, may become additional unknowns of the identification problem.

ln this paper, three methods are presented, which enable the real eigensolutions, namely, natural frequeocies and oonnaJ modes of more coostrained, undamped, MCs to be obtained indirectly from the FRFs measured on a less constrained, damped, IC. Since the structural modifications are oumerically introduced, ideal groundings of DOFs can be achieved. This way, the difficulties mentioned above can be overcome.

The methods focused in this paper have proved to be very useful in the context of a strategy of enrichment of experimental data, based oo the simultaneous exploitation of different boundary condition configurations. This strategy bas been associated with some techniques for the adjustment of finite elemeot models and for the ideotification of structural faults using the dynamic responses (Rade et ai. , l994a and Rade, l994b). Presented the 13th ABCM Mechanical Englnering Conference • COBEM 95, Belo Horizonte, December 12-16, 1995. Tectmlcal Edltorshlp: COBEM Editorial Comrntttee

375 J . ol lhe Braz. Soe. Mechanical Sciences • Vot 18, Oecembef 1996

Fonnulations of the Methods

First Method

The frrst method examined, which is applicable to the case of grounding of a single DOF, exploits the equivalence between the eigensolut ions of the MCs and the zeros of the FRFs of the I C. This equivalence is stated by the following assertion (Miu., 199 I):

"For undamped structures. the eigensolutions ofthe MC obtained by grounding the i-th DOF ofthe IC are equivalent to the antiresonance eigensolutions associated with the zeros of the FRF H ii ( ro) of the IC."

The demonstration that the zeros of point and transfer FRFs are, respectively, the solutions of symmetric and non-symmetric eigenvalue problems. is given by Flannely ( 1971 ).

According to the statement above, t~.e natural frequencies of the MC, obtained by grounding the i-th DOF ofthe IC, denoted herein as "ro ... , v = .. 1, 2, .... ,are given by the zeros ofthe FRF A;; (ro) and the corresponding eigenvectors ofthe MC, {

11Yv} • v = .J, 2, .. . ,are collinear to the i-th column of

the dynamic flexibility matrix ofthe IC, evaluated for each "rov. Symbolically:

ii Hii( ro) =O;

where { I;} denotes the i·th column ofthe identity matrix.

Considering the case of lightly damped ICs, it was demonstrated, by using a first order perturbation-type technique (Rade, 1994b), that any FRF of the IC, H;~d) (ro), is related to the corresponding FRF of the undamped IC through the expression (in the reminder, superscripts (d)

indicatc lhe quantities pertaining to the damped structures)

( I)

where lm [ Hód) ( ro)] is of the sarne order as the arbitrarily small perturbation parameter E , defined in the sense of the perturbation method. This result allows to obtain a first order approximation for the real eigensolutions ofthe undamped MC, through the application ofthe previous statement, by simply operating on the real part of the FRFs of the JC. So, lhe eigenfrequencies of the undarnped MC, detined by grounding of the i-th DOF, are given by lhe frequency values between two successive eigenfrequencies ofthe IC which satisfy:

v= 1,2, .. . (2)

Alternatively, the eigenfrequencíes of lhe undamped MC can be obtained from the phase of the FRF H;~d)( <o) according to:

(3)

The corresponding real, non normalized eigenvcctors are calculated according to:

(4)

D. A. Rade et -'·: Vlbration Analyals ot Strudures Subfedd to ... 378

where c is the number of instrumented OOFs for which the experimental FRFs ofthe IC are available.

The practical procedure for calculating the eigensolutions of MCs using this metbod is given in detail by Rade (1994b).

Second Method

This method can be used for the simuJtaneous grounding of any nwnber of DOFs. To simplify the presentation, its formulation is presented bere forthe case ofgrounding oftwo DOFs, its extension to the case of grounding of any other number of DOFs being straightforward. The case of undamped structures wiiJ be first examined.

Assuming that a submatrix of dynamic flexibilities, measured on a number c of instrumented DOFs is at disposal, two control forces ~ and ~. are applied on DOFs i and j, respectively, leading the harmonic responses on both these DOFs to vanish. The forcing frequencies and the cQrresponding displacement shapes which satisfy these conditions are defined respectively as the eigenfrequencies and the eigenvectors of the MC obtained by simultaneously grounding the DOFs i and j.

The pertaini:ng dynarnic flexibility relations are written:

[Hi 1 (ro) Hji (ro) Hjj(ro) v1=o

[H;, (ro)] H11 (ro) H1J (ro) {O}

Yi=O ..

[H 11 (ro) { Hli (ro) { H 1J(ro)} {y,(ro)} f;

The partitioning indicated in the equation above leads to the relations:

t--Hi_i <_ro)-+-_HiJ_· C_ro)--; D = ~o H 11 (ro) H;j(ro) D ~ (5)

(6)

Equations (5) and (6) fonn ao eigenvalue problem. Developing (5}, lhe following frequency equation is obtained;

(7)

The real solutions of (7), noted ijmv, v = I, 2, are the eigenfrequencies of the MC. The introduction ofthe 'Jrov into (5) leads to the following expressions for the ratios between the control forces :

( f,) H .. ( 1iro ) H .. ( ijro ) .J =~" = 11 v f. ij l ij )

I v H .. (I) H.j ()) JJ V I V

(8)

J . of lhe Braz. Soe. Med1anical Sclences- Vol. 18. December 1998

lntroducing (8) into (6). the following expression for non-nonnalized eigenvectors of the MC is obtained:

(i\,.} ... {Hii(ijro)} +(f) {H;/jro)} I V

(9)

Consider now the case of Jightly damped structures. Taking into account relation (I), the perturbation method, in first order approximation, Jeads to:

( l O)

where the imaginary part is agajn of the sarne order as the arbitrarily small perturbation parameter E •

defmed in the sense of the perturbation method. The eigenfrequencies of the undamped MC satisfy the following relations, which can be considered as homologous to (2) and (3):

v = 1, 2, ... (11)

ij (d)(ii ) + [ iio(d) (iiro ) ] = tan_1 lm[ O rov J .. ±!! .

v ReCjn(d>CiroJJ 2 ' v= 1, 2, .. . (12)

The ratios between the control forces and lhe corresponding real eigenvectors ofthe undamped MC are given by the expressions:

{13)

{ ijY) = Re( t H li (ijro)}) + ( t) Re( { Hij (1jroy)}) I V

( 14)

The practical procedure for calculating the eigensolutions of the MC by using this method is presented in detail by Rade ( 1994b).

Third Method

The method comprises two steps:

1st) calculation ofthe FRFs ofthe MC, from the FRFs ofthe IC;

2nd) application of a modal identification algoritbm to the calculated FRFs of the MC to obtain its eigensolutions.

The fonnulation of the method is presented bere for the case of grounding of a single DOF. lts extension to case of grounding of a larger number of DOFs is immediate. It should be mentioned that analogous formulations have been used by Crowley et ai. (1984) and Ewins (1989).

Let:

i be the DOF for which the harmonic response is to vanish and wbere the control ~ force is applied;

k be the DOF where the excitation force fk is applied.

The harmonic responses at these two DOFs are given by the following flexibility relations:

O. A. Rede et ai .. Vlbfation Anatysis of Structures Subjectd to ... 378

( 15)

(16)

From equation (15), the ratio between the forces f; and fk is written:

f. = - H.(d> "' (m) H.<d>(m) [k I 11 1k (17)

lntroducing (17) into (16), the following relation is obtained:

+ f~

This last equation is re-written under the fonn :

(18)

Vector { if<d>(m) }contains the FRFs ofthe MC. The a.pplication ofa modal identification algorithm to these FRFs le ads to tbe complex eigensolutions of the damped M C. 'JmSd)' {'Jy sd)} ;v .... I' 2, .. . When damping is small and the eigenfrequencies are not too close one

to the others, fairly accurate approximations to the real eigensolutions of the MC can be obtained by simply taki.ng tbe real parts of the corresponding complex eigensolutions. ln case these requirements are not fulfilled, the real eigensolutions can be either directly identified from the FRFs of the MCs by using special modal identification a.lgorithrns (Otte et al., 1993), or be calculated a posteriori from the identified complex eigensolutions by using techniques such those suggested by Zhang and Lallement (1985).

Numerical Applications

ln the sequence, the results obtained by applying the tbree methods to a full-scale muffier of an automotive exbaust line (see Fig. I), tested in laboratory, are presented.

® 'outooe.piAno­_ :la-p...,.a.-

Fig. 1 ScheiiMI of lhe Teat Stnlcture.

r,

379 J . ofthe Braz. Soe. Med\anica.l Sciences • Vol. 18, December 1~

The main characteristics ofthe test procedure are:

free-free boundary condilions provided by a flexible suspension;

lnstrumentation with 36 piezoelectric accelerometers;

3 single point random excitations, f1, f2 and f3 apptied successively and independently in the locations indicated ln Fig. 1;

acquisition of 507 frequency points in the frequency range (320 Hz; 520 Hz};

the final values of amplitudes and phases ofthe FRFs were obtained by calculating the mean values over 150 samples, for each value ofthe frequency.

Groundlng of 1 DOF

The results presented in the following were obtained by using the first and the third methods to calculate the eigensolutions of the MC defined by grounding the DOF where the excitation force f2 is applied.

Figure 2 shows the FRF H22 ( ro) of the I C, superimposed to some FRFs of the MC, which were c.alculated in the frrst step ofthe third method. lt can be clearly seen that the antiresonance frequencies of H 22 ( ro) coincide with the eigenfrequencies of the MC.

-Cl)

8 tb G)

"O '-" G)

] -100 p,.

360 380 400 420 440 460 480 500

frequency (Hz) 10'5

~ 10-e

'-" 10.7

.g .a 10-e --õ.. 10-9

~ 10·10 340 360 380 400 420 440 500

F.lg. 2 FRF H22 (ro) ofttle IC and some FRFa oftha MC

Figure 3 allows the comparison of the eigensolutions corresponding to the antiresonance frequency designated in Fig. 2 as 21ol5, obtained by using the first and the third methods. A good agreement can be observed between the results provided by both methods. lt should be noted that, since the structure is lightly damped and the natural frequencies ofthe rc are well separated, the real eigenvectors ofthe assoe iate conservative MC, resultíng from the application ofthe third method, were obtained by taking the real part ofthe oomplex eigenvectors identified in the experimental modal analysis procedure. The eigenvectors were normalized soas to render the magnitude oftheir largest componcnts equa.l to one.

O. A. Rade et ai.: Vlbration Analysis of structures Subjectd to ...

0.8

~ 0.6

~ 0.4 &. ~ ã 0.2. 8 .... o

i .0.2

5 -~ .0.4 u

.0.6

.().8

+ First Method

0 Third Method

fi.

' ' /': ..

• •

22 (1)5 = 4662 Hz 22 IDs = 466.2 Hz

• '· .·•. . .. ··•··. ••

-1oL-----~5----~1o~_.~1~5----~2~o-----25~----~37o----~35

instrumented DOFs

Fig. 3 22

Elgenaolutlona Aaaoclllled wtth (J)s

Simultaneous groundlng of two DOFs

380

The second and the third methods were used to calculate the eigensolutions ofthe MC obtained by simultaneously grounding the DOFs where the excitation forces f2 and f3 are applied.

Figure 4 shows the Bode diagram for the frequency equation, defmed in Equation (7).

ln Figure 5 the eigensolutions corresponding to the frequency 23(J)4 , indicaled in Figure 4, can be compared. Once more, a good agreement between the results provided by both methods is observed.

200~r---~----~----r---~-----r----~---,-----.~

380 <400 420 440 460 480 500

10'12 frequency (Hz)

-~ 10·13 -.g .... ~ 10

Q. 10''8

~ 10·18

340 360 380 400 420 440 460 480 500

~ 4

Flg. 4 23

Bode Dlagram for D ( d) ( (J))

O. A. Rade et ai. : Vibralion Anary.is of Structures Subjectd to ... 382

Acknowledgments

The first author wishes to express his gratitude to the agency CAPES ofthe Brazilian Ministry of Education and Sports, for the financial support which has enabled him to develop his Ph.D. prograrn in France.

References

Aitrimouch, H., 1993, • Analyse de Structures Mécaniques Modifiées", Th~e de Doctoral, Université de Franche­Comté, UFRST, Besançon, Francc.

Cogan, S., 1990, "Rtanalyse de Structurcs à Partir de Donnêes Expérimentales: Rigidiflcation et Substitution", Th~e de Doctorat. Université de Franche·Comté, UFRST, Besançon, France.

Crowley, J.R. et al ., 1994, "Direct Structural Modification Using Frequency Response Functions", Proceedings of the 2nd lntemational Modal Analysis Conference, pp. S8-6S.

Ewins, D.J., 1987, "Vibration Analysis of Modified or Coupled Structurcs Using F.R.f. Data", Internal Repport, Imperial College of Sciences and Technology, London.

Flannelly, W. G., 1971, "Natural Antiresonanccs in Structural Dynamics", Internal Report, Kaman Aerospace Corporation, Bloomfield, Connecticul

Miu, D. K., 1991, "Pbysical lnterpretation ofTransfer Funct.ion Zeros for Simple Control Systems with Mechanical Flcxibilities•, Joumal of Dynamic Systems, Measurement and Control, Vol. 113, pp. 419-424.

Otte, D. et al ., 1993, "Enhaced Force Vector Appropriation Methods for Normal Modc Testing•, Proccedings of thc IIth lntemational Modal Analysis Conference, pp. 1310-1316.

Rade, O .A., Lallement, G., Cogan, S., 1994a, "lm.provement of Structural Change Assessment by Using Multiple Boundary Condition Configurations•, Proceedings ofthe 19th lntemational Seminar on Modal Analysis, Leuven, Belgium.

Rade, D.A., 1994b, "Correction Paramétriquc de Mod~les Élements Finis - Élargissement de !'Espace de Connaissance", TMse de Doctorat, Universitê de Franche-Comté, UFRST, Besançon, France.

Snyder, 'V. W., 1986, "Structural Modification and Modal Analysis- a Survey", lntematlonal Journal of Analytical and Experimental Modal Analysis, Vol. I, n° I, pp. 45-52.

Zhang, Q., Lallement, G., 1985, '"New Method of Determining the Eigensolutions of lhe Associate Conservative Structures from the ldentified Eigensolutions• , Proceedings of the 3rd lnternational Modal Analysis Conference, pp. 322-328.

381 J . ot lhe Braz. Soe. Mechanical Sciences • Vol. 18, Decembef 1996

•• •

+ Second Method

O Third Method

• • •

. :

• •

• :· : ~

~ : .

•

•

23 (1)4 "' 452.6 Hz 23 (1)4 = 452.6 Hz

. .. : . ••

•

• -tL---~----~--_.~----~----2~----30~----35~

5 10 15 20 5

instrumented DOFs

23 Fig. 5 Elganaolutlona Anocl1ted wlth cs>4

Discussion and Conclusions

lbe various nwnerical applications perfonned have shown that the thrce methods studied are well adapted to the Analysis of Modified Structures based on the use of experimentally measured FRFs. The accuracy ofthe results obviously depends upon tbe quality ofthose FRFs. Applications performed to numerically simulated structures - for which tbe exact eigensolutions of tbe MC are known - have shown that an acceptable accuracy can be achieved even for relatively high leveis of simulated experimental noise contarninating the FRFs ofthe IC (Rade, \994b).

The good correlation obseived between the eigensolutions provided by the different methods does not assure the accuracy of these eigensolutions. Nevertheless, it reveals that, as expected, the different methods perform in an equivalent way, despite the fact that each of them follows clearly distinct calculation procedures. Furthermore, the coherence between the results can be seen. at least. as a necessary condition for their accuracy.

Severa! authors (Ewins, 1989, and Aitrimouch, 1993) have perceived, when using AMS techniques similar to the third method presented in this paper, the appearance of spurious poles in FRFs of the MCs. which nearly coincide with the eigenfrequencies of the IC. These authors ascribe this phenomenon to: i) the incoherences affecting the measured FRFs, caused by experimental noise and measurement erros, and ii) the numerical iii conditioning of the submatrix of dynamic flexibilities, which _has to be inverted in Equation ( 17), in the case of simultaneous grounding of severa! DOFs. Contrarily to the reports of the mentioned authors, such spurious poles have not been observed for the experimental test structure examined in this paper. lt is believed that this is dueto the excellent quality of the FRfs of the IC, for which the v alues of amplitude and phase were obtained by calculating the mean v alues over a large number of samples, for each frequeocy line.

Some aspects conceming the practical application of the methods focused in this paper deserve further investigation a.nd are currently being examined, such as : i) the performance of the second and third methods in the case of simultaneous grounding a larger number ofDOFs, and ii) the evaluation of the accuracy of the eigensolulions ofthe MC calculated from experimentally measured FRFs.

RBCM - J . of the BI"'IZ. Soe. Mechanlc•l Sclencea Vol. XVIII- No. 4 • 1996 • pp. 383-394

Gerando o Caminho de Corte de Reentrâncias Utilizando Diagramas de Voronoi

ISSN 0100-7386 Printed ln Brazll

Use of Voronoi Diagram for Defining the Cutting Path Marcos S. G. Tsuzuki Lucas Antonio Moscato Universidade de Sto Paulo ESCOla Potitéalica Departamento de Engenharia MecânicaiMecatrOnlca Av. Prof. Mello Moraes. 2231 05508· 900 Sto Paulo, SP Brasil

Abstract

ln tbis work the use ofVoronoi Diagram for defining tbe cuning path is discussed. ln tbis fashion tbe cutting path is deterrnined in two steps· 6rst, tbe Voronoi Diagram is created, and secood, the cuttiog patb is calculated based on lhe Voronoi Diagram. A possible implementatioo for both steps is presented. The main purpose ofthis algorithm is to facilitate the comprehensioo of the evolved concepts. Resuhs obtained from a prototype are presented. KeywoniJ: Pocket Milling. Cutting Path. Voronoi Diagram

Resumo Neste trabalho realizamos uma rápida introduçilo, onde destacamos a ímportância das pesquisas na área de usin.agem 2Yill Em seguida apresentamos uma rápida revisilo bibliogrâfica do assunto. Detalhamos a proposta de definir o caminho de corte com base nos diagramas de Voronoi. Assim. a definiçilo do caminho de corte passa a ser feita em duas etapas: primeira, definir o diagrama de Voronoi e, a segunda. definir o caminho de cone propriamente dito. Apresentamos um algoritmo para cada uma das etapas. Estes algoritmos foram desenvolvidos com o principal objetivo de facilitar a compreendo dos conceitos envolvidos. Finalmente, apresentamos alguns resultados obtidos com o protótipo desenvolvido e tiramos algumas conclusões. Palavras chave: Caminho de Corte, Diagrama de Voroooi, Usinagem de Reentrâncias

Introdução

Os principais tipos de usinagem por comando numérico geralmente são classificados de acordo ~om o número de eixos cont.rolâveis simultaneamente e independentemente um do outro:

Usinagem 20: permite o controle de dois eixos rranslacionais. Logo, a usinagem 20 suporta apenas aplicações planares,

Usinagem 3D: permite a interpolação linear utilizando todos os três eixos translacionais simultaneamente. Interpolações circulares podem ocorrer apenas em um plano de coordenadas (planos -xy, -yz ou -xz).

Usinagem 50: adicionalmente às caracterlsticas da us inagem 3D, dois eixos adicionais de rotação são controláveis.

A maior dificuldade em usinagens 50 e 3D (este último em menor escala) é a necessidade de dados de controle precisos para o controle simultâneo dos vários eixos de movimento. Estes mesmos problemas estão presentes no campo da robótica. Detetar e evitar colisl\o entre a ferramenta e o "mounting", de um lado, e a peça e o fixador, de outro lado. são problemas matematicamente e algoritmicamente complexos. Devido a estes problemas, a usinagem 5D é utilizada apenas quando não existe nenhuma outra solução posslvel. Entretanto, usinagem 20 sofre uma série de restrições drásticas que a eliminam de ser uma alternativa prática para usinagem 50. Este vazio entre os requisitos práticos, de um lado, e as deficiências de programação, de outro, é sanado por uma forma hlbrida entre usinagem 20 e usinagem 30: usinagem 2~0. Em principio é passivei realizar usinagem 3D, utilizando usinagem 2 '1,0. A Figura I ilustra os vários tipos de usinagem apresentados.

Manuscript received: October 1995. Technlcal Editor: Agenor de Toledo Fleury

J. of lhe Btaz. Soe. Mechanical Sclences- Vol. 18, December 1996

Fig. 1 01 Prtncl~ Tipo. de Uelnagem: 20, zy.o. 3D • 50.

A grande maioria dos processos de usinagem industrial pode ser realizada utilizando usinagem 2~D. segundo Held Keppler ( 1989) mais de 80% de todas as peças mecânicas podem ser usinadas segundo o conceito de usinagem 2Y.O, ele explica esta afirmação pelos dois aspectos seguintes:

Surpreendentemente, um grande número de peças mec4nicas possui a forma de curvas de nlvel em terraços, i.e. suas faces de contorno silo simultaneamenle paralelas ao plano-xy ou constantemente normais ao plano-xy. Por convenção, estes objetos são denominados como 211, dimensionais.

Objetos mais complexos- como cavidades com contornos definidos por superflcies esculpidas -são produzidos usualmente por um desbaste 2Y.O e um acabamento 3D ou 50.

Desta maneira, conclulmos que o domlnio de técnicas de usinagem 2~0 é de fundamental importância

Revisão Bibliográfica

As técnicas conhecidas para usinagem de reentrâncias podem ser classificadas de acordo com a estratégia de usinagem e pela maneira que o caminho de corte é definido. Apesar do caminho de corte variar de algoritmo para algoritmo existem apenas dois tipos de estratég.ias de uslnagem -vide Fig. 2:

Usinagem paralela ao contorno. e

Usinagem paralela a uma direção (zigue z.ague1)

Fig. 2 A Uslnagtm Paralela ao Contorno • a Uslnagem Paralela 1 uma Dlreçlo.

1 Existe uma variaçlo da usioagem zigue-zague chamada de usinagcm :r:igue-zague, onde a usinagcm ocorre cm apenas um sentido, de modo que a usiugem nJo OCOITI na direçlo oposta à da ítvore da màquina de comando numérico ("spiodle").

385 Tsuzukl M. S G. et ar.: GerandO Caminho de Corte de Reentrâncias ...

A estratégia de usinagem paralela ao contorno utiliza o conceito de deslocar elementos do contorno paralelamente para definir o caminho de corte. lsto significa que a área é usinada segundo a forma de uma espiral. A idéia bâsica da estratégia de usinagem paralela a uma direçllo é aparentemente simples: após selecionarmos uma linha de referência inicial a usinagem ocorre utilizando segmentos paralelos a esta linha de referência. A principal desvantagem desta estratégia é a necessidade de realizarmos um ciclo de usinagem.

Podemos dividir as estratégias de usinagem paralela ao contorno em duas técnicas principais:

Definindo o "offset" do contorno a partir do contorno das reentrâncias,

Definindo o ''offset" do contorno a partir de diagramas de voronoi.

Os algoritmos que implementam a primeira técnica estão concentrados na definição de "offsets" sucessivos a partir do contorno original. Este processo é executlldo segundo os três passos seguintes -vide Fig. 3:

Para cada elemento do contorno um elemento de "offset" elementar é construido;

Espaços entre elementos de "offset" consecutivos são preenchidos por arcos de conexão, isto pode resultar em curvas fechadas e, em alguns casos, auto-intersectantes;

As auto-intersecções das curvas são eliminadas e partes da curva são eliminadas determinando a curva de "offset" final.

O problema principal nesta técnica é a necessidade de determinarmos todas as auto-intersecções. Até o momento, as vllrias propostas conhecidas n!o eliminaram a necessidade de interseccionar todos os pares de elememos de "offset" entre si, além de apresentarem complexos algoritmos para eliminar as autointersecções (Suh e Lee, 1990). Vários artigos sobre o assuntojâ foram publicados, inclusive no Bras i I. Ferreira ( 1993)propõe uma técnica que denominou por "shrinking" que é baseada na teoria de grafos. Cota et ai. { 1993) propõe uma técnica baseada em subcontornos que são deslocados paralelamente ao contorno original.

Fig. 3 Definindo o "Of'lllet• do Contorno Dlretllmenta a Partir do Contorno da Raentrtncla, Cuo com VAn .. Auto-lntllrsac:çOea.

Seguindo uma linha complementar de pesquisadores, Tsuzuki e Miyagi [ 1991 ] apresentaram uma técnica para definir o caminho de co rte baseando-se em técnicas de reconhecimento de "Form Features" diretarnente a partir do Modelo Sólido. As técnicas de reconhecimento de "Form Features" baseiam-se no conceito de que é posslvel associar um conjunto de comandos de usinagem a uma forma geométrica simples, e é possivel interpretar uma forma geométrica complexa como sendo a composição de formas geométricas simples. Este último é um problema que ainda permanece aberto. Toledo eL ai . ( 1993) propuseram uma técnica que fornece um conJunto de interpretações a uma forma geométrica complexa e ficaria a cargo do usuário decidir qual a interpretaçllo mais adequada. Entretanto, o problema de múltipla interpretação de formas geométricas complexas ainda está longe de ser resolvido.

Diagramas de Voronoi

Person, ( 1978) propôs uma nova e eficiente técnica baseada em duas etapas. Inicialmente, toda a reentrância é dividida em sub-áreas independentes c, em seguida, o caminho de corte é criado diretamente a partir das sub-áreas. Person apresentou em seu artigo a seguinte idéia: selecionc um ponto final p entre dois segmentos de caminho de corte consecutivos, e uma observação óbvia, mas importante, é que este ponto possui a mesma distância miníma de dois elementos de contorno e está a uma distância maior de qualquer outro elemento de contorno

J . ofthe &az. Soe. Medlanlcal Sclenoes • Vol. 18, Deoember 1996

Fig. 4 Propriedade element.r doa porrtoe flnala entre doia aegmentoe de caminho de corte coneecuthlo. e o reapectlvo diagn~ma de Voronoi.

386

Considerando o conjunto de pontos que possuem esta propriedade obteremos um grafo muito conhecido na geometria computacional como diagramas de Voronoi. Apenas para auxiliar a intuição, consideremos um valor de "offset" inicial que aumentará continuamente até que a curva de "offset" degenere tomando-se um ponto. O conjunto dos pontos p definidos segundo as curvas definidas com os vârios valores de "offset" d definem o diagrama de Voronoi. A Figura 4 ilustra este conceito.

Qual será a vantagem em considerar dia.gramas de Voronoi? A disponibilidade dos diagramas de Voronoi facilita consideravelmente a criação das curvas de "offset". Bastaria definir um algoritmo que defina os pontos extremos de uma curva de "otfset" elementar:

Construa um segmento de "offset" elementar a partir do elemento de contorno, e

lnterseccione o segmento de ·•offset" e lementar com os bissetores que delimitam a área de Voronoi associada ao segmento de "offset".

Se o diagrama de Voronoi for representado de maneira conveniente então este algoritmo poderá ser executado de maneira muito eficiente. Em particular, a idéia originalmente proposta por Person foi que os bissetores devem ser expressados por funções cujos parâmetros representam a distância mínima do bissetor ao contorno da reentrância, ou simplesmente, distância de "offset". Desta maneira, a determinação das intersecções do passo 2 do algoritmo será reduzida a simples avaliação das fórmulas parametriz.adas dos bissetores. Segundo Held [ 199 1 ], e le foi o primeiro a publicar e implementar esta idéia, originalmente proposta por Person, segundo um conceito mais omplo para reentrâncias com contornos complexos.

Determinando Diagrama• de Voronoi

As seguintes suposiçOes foram realizadas:

É comum que os contornos das reentrâncias possuam casos matematicamente problemáticos, como segmentos de linha paralelos, arcos concêntricos, entre outros. Desta maneira, todos estes casos devem ser previstos.

O contorno da reentrância é definido pela composição de segmentos de linha e arcos de circunferência.

A maioria das máquinas de comando numérico permitem que realizemos apenas interpolações lineares e interpolações circulares. Desta maneira, a segunda suposição está dentro de nossas necessidades do momento.

Parametrizaçlo dos Blssetores

Uma vez que a distância de "offset" é muito utilizada durante todo o sistema, tanto no que diz respeito à construção do diagrama de Voronoi como na determinação do caminho de corte, decidimos expressar os bissetores como funções de seu "offset" em relaçito aos elementos de contorno. Um problema porém surge nesta representação: um bissetor b pode possuir dois pontos PI e P2 com o mesmo "offset" em relação aos elementos de contorno- vide Fig. 5.

Pi"---JP2

Fig. 5 Pontoa P1 e P2 Pouuem o Meamo "Offaet~.

387 Tauzukí M. S. G. et ai.: Gerando CaminhO de Corte de Reentrâncias ...

A solução adotada em nosso sistema foi proposta por Held [ 1991 ], ela consiste em dividir os bisseto~s em dois tipos: bissetores geométricos (as arestas de um diagrama de Voronoi) e bissetores anallticos (bissetores geométricos se dividem nos pontos de extrema distAncia de "offset"- mini ma ou mâxima distância).

Como é posslvel parametrizar um bissetor analltico? Inicialmente, nllo distinguimos entre um segmento de linha limitado e sua reta suporte. ou entre um arco circular e o circulo que contém o arco. As linhas e clrculos são representados pelas sua equações impllcitas, e desta maneira é muito fãcil representar as curvas de "offset" . Para um circulo temos:

(x - xc)2 + (y- yc)2 =r onde o "offsef' do círculo- com "offset" t., é fornecido por:

(x(t)- xc)2 + (y(t) - yc)2 =(r+ k.t)2

em que (xc, yc) representam o centro do circulo e r representa o raio do circulo. Por analogia, para a reta:

a.x + b.x + c = O

onde n2 + bl = O. O"otfset" da reta é fornecido por:

a.x(t) + b.y(t) +c+ k.t = O

onde a , b e c silo os coeficientes normalizados da reta. Em ambas as fórmulas a direção de escorregamento é dada por k. O caso em que k é positivo pode ser interpretado como um engrandecimento do círculo, e opostamente, o caso em que k é negativo significa uma diminuição do círculo. Para a reta, o parâmetro k indica se o "offset" da reta está situado à direita ou à esquerda da reta.

As fórmulas de parametrização para os bissetores (utilizando o "offset' ' t como parâmetro) pode ser obtida pela solução das equações de intersecção dos elementos de "offset" apresentados. As fórmulas finais já foram apresentadas por Person [19781 e constam do Apêndice I. Observe que estas fórmulas são vàlidas apenas para retas que nllo são paralelas e para círculos que não são condntricos.

Representando o Diagrama de Voronol

É necessário definir uma representação adequada para os diagramas de Voronoi. Como o diagrama de Voronoi é um grafo planar conexo, é natural que analisemos representações de grafos conhecidas. As mais apropriadas são: estrutura "winged-edge", proposta por Baumgart [1975] e descrita em detalhes por Mãntylll. ( 1988]; e a estrutura "doubly-connected edge list" (DCEL) proposta por Preparata e Shamos (1988).

d [) i

·<J· Aresta ant_2 ant_1 prox_2 prox_1 Faoe_esq Faoe_dir

3 2 d a

2 1 3 c d

3 2 5 4 c a 4 3 5 b c

5 4 3 a b

Flg 6 Um Dlagl"'ma de Voronol e aua RMpec:Uva ~rtSent.çio. O. ponttlroa ant_1, ant_2, prox_1 1 prox_2 Repreaentam Al'lat .. que Incidam na Preaentt Al'lata 1 Face_Mq 1 Face_dlr alo Pontelroa Pll"' ArMtas do

Contorno que Dlflnam 1 P,...entt Aresta

J . of lhe 8raz.. Soe. Mechanlcal Sdences. Vol. 18, December 1996 388

O principal componente de ambas as estruturas de dados é a aresta. Em nossa aplicação uma aresta possui associada a si as seguintes informações: quatro ponteiros para as arestas que incidem na presente aresta e dois ponteiros para os objetos que definem a presente aresta. Desta maneira, cada bissetor está associado ao contorno da reentrância. A Figura 6 ilustra um diagrama de Voronoi e sua respectiva representação.

Cone de Influência

O cone de influência é definido para os elementos de contorno da reentrância. Para um arco de circunferência menor que 1800 , o cone de influência é a região delimitada pelo par de raios originários do centro do arco e que passam pelas extremidades do arco sendo a parte hachurada da Fig. 7.a. Para um segmento de reta, o cone de influência é a região delimitada pelas retas normais ao segmento que passam pelas suas extremidades, como ilustrado na Fig. 7.b.

Detenninando as Intersecções dos Bissetores

Determinar a intersecção de dois bissetores anallticos b1 e b2, com parametrização de coordenadas (x1• y1) e (x2, yJ, significa solucionar um conjunto de equações lineares para t1 e~:

x1(t1) = x2(~

y,(t,) = Y2(~)

Como em nosso algoritmo, t1 e t2 devem ser iguais, porque ambos os parâmetros constituem o "offset" do mesmo ponto, determinar a intersecção do bissetor equivale a determinar o ponto equidistante de três objetos (ponto, reta ou circulo). Por exemplo, no caso de dois bissetores defmidos por três arcos circulares devemos solucionar o seguinte sistema de equações:

(x(t)- xc1)2 + (y(t)- yc1)2 = (q + k1.t)Z

(x(t) • XCV2 + (y(t) • y~ =(r, + k,.I)Z

(x(t) - XC))2 + (y(t) - YC-))2 = (r3 + k3.t)2

No Apêndice 2 apresentamos os passos para determinar a solução dos vários casos. Como possu[mos dois objetos posslveis, foram deduzidas quatro intersecções posslveis: três clrculos, dois clrculos e uma reta, um circulo e duas retas e, finalmente, três retas. Também foram estudados os casos em que dois clrculos são concêntricos, ou, duas retas são paralelas. Geralmente, ao se resolver o sistema não linear, resultamos numa equação quadrática. Desta maneira, há a necessidade de verificarmos qual dos dois pontos realmente pertence ao diagrama de Voronoi. O ponto que pertence ao diagrama de Voronoi deve pertencer ao cone de influência dos três elementos de contorno que definem estes bissetores.

Algoritmo Proposto

Vários algorimos para a defmiç!o do diagrama de Voronoi já foram propostos [Held, 1991, Srinivasan e Nackman, 1987, e Preparate e Shamos, 1988). Entretanto, para este projeto preocupamo­nos em definir um algoritmo de fácil compreensão e que apresentasse os conceitos propostos inicialmente por Person [1978]. A Figura 8 ilustra os vários passos do algoritmo de definição do diagrama de Vorono~ proposto neste trabalho.

O algoritmo é definido abaixo:

i. Todos os ângulos internos maiores que silo transformados em arcos de circunferência de raio nulo.

ii. Os bissetores de todos os pares de arestas adjacentes são determinados

389 Tsuz:ukl M. S G. et ai. Gerando Caminho de Corte de Reentrâncias ...

Fig .. 7 Cone de lnflutncla de um Arco de Cln:unfertncla e de um Segmento de Reta

Fig. 8 Oa V6rfoa Pauoa na Deflnlçlo doa Blueto'" que CompOem o Diagrama de Voronol.

iii. Determinamos a intersecção entre todos os bissetores que possuem uma aresta em comum, e selecionamos a aresta cujo ponto de intersecção possui o menor valor de "offset" t. Cada bissetor possui duas intersecções, uma com cada bissctor adjacente, assim devemos verificar qual o ponto de intersecção está mais próximo do ponto inicial do bissetor. Após esta verificaçllo, o ponto mais próximo ~ selecionado. Esta verificação ~ necessária pois em certos casos (bissetor de segunda ordem), o bissetor está direcionado no sentido do maior para o menor "offset". Neste caso, se outro bissetor interseccionar o bissetor em análise entre seus pontos final e inicial, então devemos selecionar o ponto de intersecção mais próximo ao ponto inicial do bissetor. A Figura 9 ilustra um exemplo onde esta verificaçllo se faz necessária

Fig. 9 Ellemplo de Crftçlo do Dltgrama de Voronol em que 1 Verffiellçlo Definida no Ten:eiro Puao 3 do Algoritmo se Faz Neceadrlt. Poaaulmoa taac> t, mu o Ponto

Definido por Tese • mala Próximo de T que o Ponto Definido por t.

iv. Este lado~ removido da lista auxiliar de arestas. Desta maneifa, as arestas adjacentes a esta aresta removida, passam a se comportar como se fossem adjacentes.

v. Determinamos o bissetor associado ao novo par de arestas adjacentes que surgiu com a remoção da aresta do te rceiro passo. E realizamos novamente este passo, att que possuamos apenas uma única aresta na lista auxiliar de arestas.

O conceito básico neste algorimo é definir o diagrama de Voronoi através da busca de pontos de intersecção entre bissetores com menor valor de "ofTset" para pontos de intersecção entre bissetores com maior valor de "offset''.

J o4 lhe Braz. SOe. Mechanical Sdeoces- Vol . 18, Deoember 1996 390

Lista Auxiliar de Arestas do Contorno da Reentrlnda: As arestas do contorno da reentrância estão annazenadas em uma lista auxiliar. Esta lista permite determinannos o estado atual do algoritmo. Após deteri'JÜnar que dois bissetores devem se interseccionar, removemos a aresta comum aos dois bissetores. Com a remoção da aresta, um novo par de arestaS adjacentes surge e utilizamos esta informação para criar um novo bissetor associado a este novo par de arestas adjacentes. E , assim, procedemos continuamente.

Verificaçlo do Funcionamento do Algoritmo: Também implementamos uma função que permite verificar se o algoritmo está caminhando corretamente. Esta função verifica se o ponto de intersecção entre dois bissetores obedece à propriedade abaixo:

d(P,L) = d(P,prox_J) = d(P,prox_2) < d(P,L)

onde L é a aresta comwn aos dois bissetores, as arestas Lccw e Lcw sll.o adjacentes à aresta L. O ponto P é o ponto de intersecção entre os dois bissetores e L' são todas as demais arestas do contorno da reentrância. Esta propriedade significa que o ponto de intersecção entre dois bissetores deve estar à miníma equidistância entre as arestas associadas a estes dois bissetores.

Resultados

Conforme apresentamos anteriormente, o diagrama de Voronoi, representado por funções cujos parâmetros representam a distância mínima do bissctor ao contorno da reentrância, permite definir um algoritmo de geração de caminho de corte simples, conforme ilustrado abaixo:

O usuário especifica um raio de ferramenta qualquer, este raio é adotado como o "offset" inicial. Determinamos o ponto mais interno com o seu respectivo "offset", desta maneira é posslvel calcular o n(Jmcro de passadas da ferramenta. O contador de passadas é inicializado com l.

Através do raio da ferramenta e do número da passada aLuai, determinamos quais bjssetores possuem "offset" inicial menor que o nlvel do "offset" atual e quais bissetores possuem "offset" final maior que o nlvel de "offset" atual.

Monta-se uma lista contendo a aresta, os pontos final e inicial e o "offset" correspondente. Esta lista é ordenada segundo a ordem crescente de identificação da aresta. Desta maneira, o caminho de corte é traçado sempre no sentido anti-horário.

Fig. 10 Olagl"'maa de Voronole Caminho de CortAI oetennlnados por Nosso Protótipo.

391 Tsuzukl M. S. G. et ai.: Gerando Caminho de Corte de Reentréncias ...

v,.... ""'<::ll

~

~

~

Flg.11 Exemplo do Processamento de Reentrtncla com Ilha.

O sistema criado a partir deste algoritmo foi testado em um grande número de reentrâncias dos mais diversos tipos que se enquadrassem nas limitações do projeto. Diversos resultados estão ilustrados na Fig. 10. No caso de ilhas, podemos unir o contorno externo da reentrância à ilha. A Fig. li ilustra um exemplo de reentrância com ilha.

Conclusões

Os maiores problemas encontrados para o cAlculo do diagrama de Voronoi e do caminho de corte estavam relacionados com erros de arredondamento. Para as pessoas que estejam interessadas em desenvolver um sistema semelhante ao que apresentamos neste artigo, podemos fornecer alguns conselhos para solucionar este problema. Primeiro, defina uma tolerância para comparar números reais. E. segundo, defina tolerâncias geométricas para que retas sejam consideradas paraJelas e para que circunferências sejam consideradas con~ntricas.

Um fator muito positivo dos diagramas de Voronol, é a definição do caminho de corte em dois passos. A definiç!lo do caminho de corte pode ser aprimorada de forma independente do aJgoritmo que determina o diagrama de Voronoi. Desta maneira. é posslvel acrescentar um algoritmo que defina a profundidade radial de avanço da ferramenta para que a vibração seja minimizada. conforme apresentado por Tsai et ai. ( 1991 ) .

Agradecimentos

Este projeto foi parcialmente suportado pelo CNPq. Agradecemos aos Engenheiros Cláudio J. F. Alves e Roberto K. Sato pela implementação do protótipo do sistema que gera os diagramas de Voronoi.

Referências

Baumgan, B.G., 1975, "A Polihedral Representation for Computer Vis ion", ln Proc. Natiooal Computer Sciencc Conference, pp. 44:589-596, AFTPS Press, Arlington. VA, USA.

Cota, F.E.; Queiroz, A .A.: Filho, E.V.G., 1993, "Usinagens de Cavidades em Formas Arbitrárias em MAquinas Ferramentas de Comando Numérico", Anais do XII Congresso Brasileiro de Engenharia MecAnica, vol. III, pp. 1599·1602, Brasllia

Ferreira, J.C.E.,I993, "An A1gorithm for Generating the Tool Paths for Machining Complex 2V.O Components", Anais do Xl1 Congresso Brasileiro de Engenharia Mec!nica, vol. UI, pp. 1587-1 S90, Brasil ia.

HeM, M., Kepplc:r. J., 1989, "GEOPOCKET. A sopbisticated oomputational geomctry soluction of gcometrical and techno1ogical problems arising from pocket machlnlng", Computer Applications ln Production and Engineering · CAPE89, Kimura. F .. Rolstadls, A. (editors), Elsevier Sciencc Publisbers B.V. (North Holland), IFIP

Held. M., 1991, "On the Computational Geometry ofPocket Machining", Springer-Verlag.

J. ol lhe Braz.. Soe. Mect\anlcal Sciences - Vol 18. Qecember 1~ 392

Mantyla, M., 1988, "An lntroduction to Solid Modeling". Principies ofComputer Science Series, Computer Science Press, Rockville, ML, USA.

Person, H., I 978, "NC machioing of arbitrarity shaped pockets~, Computer Aided Design, 10(3): 169-174, May. Preparata, F.P., Shamos, M.J., 1988, "CompuWional Geometry- An Jntroduction", Texts and Monographs in

Computx:r Scienoe, Springer- Vertag, second edition, October. Srinivasan, V., Nack:man, R., 1987, "Voronoi Diagram for Multiply Connected Polygonal Domains, 1: Algorithm",

IBM J. ofResearcb and Developmenl, 31(3):36l-3n, May. Suh, Y .S., Lee, K.,June 1990, ~NC milling tool path generatíon for arbitrary pockets defined by scu1ptured

surfaces", Computer Aided Oesign, 22(5):273·284. Toledo, C.F.M.. Garcia, M.B., TsiU.Uici, M.S.G., Miyagi, P.E., 1993, "CAPP AutomAtioo", Anais do Xll Congresso

Brasileiro de Engenharia Mecânica, vol. UI, pp. 1767-1770, Brasília. Tsuzuki, M.S.G., Miyagi, P.E., 1991, "Reconhecimento AutomAtioo de "Features" Baseado na RepresentaçAo B·

REP", X1 Congresso Brasileiro de Engenharia MecAnica, vol. 8, pp. 611-614, Slo Paulo. Tsai, M.D., Takata. S.; lnui, M., Kimura, F.; Sata, T., 1991 , "Operation Plannlng Based on Cuttlng Process

Models", Annals ofthe ClRP, Vo140/1, pp. 95-98.

Apêndice 1

As fórmulas de parametrizaç!o para os bissetores (utilizando o "offset" t como parâmetro) podem ser obtidas pela soluçAo das equações de intersecçao dos elementos de "otTset" reta e circulo. As fórmulas finais já foram apresentadas por Persoo (I 978).

Reta-Reta

a 1.x(t) = b1.y(t) + c1 + k1.t =O

a1.x(t) = b2.y(t} +c,+ k1.t =O

onde,

Bissetor Linear

onde

Círculo-Reta

(x (t)- xc1f+ (y (t}- yc1f = (r1 + k1.1)2

~ . x. (t) + ~ . y (t) + c2 + k2.t = O

onde

393

Bissetor Parabólico

onde

r1(t):= r1 + k1.t

h:=~.xc1 + b2.yc1 + C2

h(t):= h + kl.l

Círculo-Círculo

Tsuzukl M. S. G e1 ai .: Gerando Caminho de Corte de Reenlrtndas ...

{x{t) • xc1) 2 + (y(t)- yc1)l = (r1 + k1.tY

(x(t)- xc2)2 + (y{t)- yc2)2 = (r2 + k2.tY

Bissetor Hiperbólico/Eifptico

r1(t):= r 1 + k1.t

r2(t):"' r2 + k2.t

onde dx:= (xc2 - xc1) I d

dy:= (yc2 - yc1) I d

6 :={k2.r2 - k 1.r1) I d

h:= (ri- rr- d2) I (2d)

h(t):= (r2(t)2 • r 1(t)2 • d2) /2d

Apêndice 2

Como possulmos dois objetos posslveis, foram deduzidas quatro intersecções possfveis: três clrculos, dois circulos e uma reta, um circulo e duas retas e, fmalmente, trts retas. Considere o caso em que os dois bissetores silo defnidos por tres arcos de circunferência:

(x(t)- xc1)2 + (y(t)- yc1) 2 "'(r1 + k1.t)2

(x(t) • xc2) 2 + (y{t)- yc2)2 = (r2 + k2.1)2

(x(t)- xc3) 2 + (y(t) • yc3)2 = (r3 + k3.tf

J . of lhe Braz. Soe- Mec:hanlçal Sclences- Vol 18, December 1~ 394

Através de um método de substituição de variáveis, x(t) e y(t) podem ser expressos em termos lineares de t, ou seja x(t) = a · t + 13 e y(t) y · t + õ. Depois de substituirmos essas expressões em uma das fórmulas dadas, um polinOmio de segundo grau é obtido. A solução deste polinOmio de segundo grau em t nos fornecerá os valores que definem a intersecção dos dois bissetores O caso em que possulmos dois arcos de circunfer~ncia e uma reta pode ser solucionado de modo semelhante ao método proposto para o primeiro caso.

O caso em que os dois bissetores são defmidos por um arco de circunferência e duas retas:

(x(t)- xc1f + (y(t)- yc1)2 = (r1 + k1.t)l

~ . x(t) + ~ . y(t) + ~ + k2 . t = O

a3 . x(t) + b3 • y(t) + c3 + k3 . t =O

O sistema formado pelas duas equações lineares é resolvido pela regra de Cramer. Determinando os valores de a, 13 , y e ô recai mos novamente no método proposto no primeiro caso. E para o último caso em que os dois bissetores são deímidos por tr!s retas, utilizamos diretamente a regra de Crarner para encontrar a solução do sistema.

RBCM • J . of lhe Braz. Soe. Mechanlcal Sclencea Vol XVIII · No.4 • 1P96 • pp. 3115-402

ISSN 0100.7388 Printed ín BrazU

Projeto de Cascatas de Pás para Turbina Baseado num Método de Progressão no Tempo Design of Turbine Blade Rows Using a Time-Marching Method João Eduardo de Barros Teixeira Borges Luis Manuel de Carvalho Gato Rui Manoel Roma de Jesus Pereira UniversJdade Técnica de Usboa Instituto Superior Tácnloo Departamento de Engenharia Mecânica 1096 Uaboa Codex Portugal

Abstract

An invcrse method based on the itera.tive use of a direct time-marcbing code is described. lbe design input specification used consists of lhe imposition of the mean tangencial velocity distribution and the desired distribution of blade lhickness normal to the blade ca.mber given as a funclion of the axial distance. The described inverse method was applied to the design of a turbine cascade typical of a rotor blade row representing a high aerodynamic load. This exa.mple brings out clearly the relationship between the blade pressure loading and thc dcrivation o f thc meand velocity distribulion used as design input. Keywords· CFD, Turbomachine, lnverse Methods, Turbine Cascade

Resumo Descreve-se um método inverso baseado no uso iterativo de um programa de análise do tipo Progressão no Tempo. A e.specificação de desenho para o escoamento consiste na imposição ao longo da distAncia axial da distribuiçao da velocidade tangencial média, e na distribuição desejada da espessura normal à linha de esqueleto. O método de projeto desenvolvido foi aplicado ao desenho de um rotor de uma turbina sujeito a uma carga aerodinArnica elevada. O exemplo mostra claramente a relação existente entre diferença de pressões entre as duas faces da pá e a derivada da velocidade tangencial médla utilizada como uma especificação de projeto. Palavru chaves: CFP, Turbomáquinas, Métodos Inversos, Cascatas de Turbinas.

Introdução

O projetista confrontado com a tarefa de desenhar turbomáquinas eficientes pode recorrer à métodos inversos em vez de utilizar o processo mais simples e usual de adivinhar uma geometria e analisar o escoamento produzido por esta geometria utiJ izando um código de análise, a fim de verificar se as caracterlsticas do escoamento assim produzido são aceitáveis. As metodologias inversas tratam o problema de fonna oposta, na medida em que primeiramente se especificam as caracte.risticas que se deseja para o escoamento, e se calcula em seguida a geometria da cascata de pás que pennite obter essas caracterlsticas. Nos últimos anos tem-se evidenciado um interesse renovado no estudo dos métodos inversos, o qual se manifesta no elevado número de artigos publicados e conferências cientificas organizadas sobre o assunto- veja-se a titulo de exemplo Dulikravich ( 1991, 1992), Van den Braembussche ( 1990), Slooff ( 1989). Como res ultado de todo este esforço, pode-se dizer-se que surgiram três grandes tendências para a solução do problema inverso: técnicas de otim izaç4o, métodos de fonnulação completamente inversa e uso iterativo de códigos de análise.

Nas técnicas de otimização descreve-se a geometria da coroa de pás em função de alguns parâmetros e procura-se entllo. no espaço de desenho definido por esses parâmetros geométricos, a combinação de parâmetros que dê as caracterlsticas do escoamento mais próximas posslvel das desejadas ou que conduza a um ótimo no que diz respeito a uma função-objetivo dependente do escoamento. Uma fonna usual de por em prática os métodos de otimizaçâo consiste em minimizar as diferenças entre a distribuição desejada de pressllo nas superficies da pá e a calculada para urna dada

Presented at lhe 13th ABCM Mecha nicai EngiMering Conference. COBEM !15, Belo Horizonte, Oecember 12-15, 1995. Tedlnical Editorship: COBEM Editorial Commlltee

J ot the Braz. Soe. Mechanical Sc:íences ·VoA. 18, Oeoember 1996 396

combínaç!o de parâmetros. Em contraste, nos m~todos de formulação completamente ínversa calcula­se d iretamente a geometria da coroa de pás que nos dá o escoamento pretendido, utilizando uma formulação do problema que faz uso imediato da ínfonnaçâo correspondente às condições desejadas para o escoamento (por exemplo, através da imposiç~o destas condições como condições de fronteira para as equaçOes diferenciais a resolver). A estratégia de usar iterativamente códigos de análise consiste na junçllo de um simples programa de cálculo do escoamento com um procedimento automático de a.lteração da geometria da cascata de pás, por forma a atingir os objetivos desejados. Nesta estratégia começa-se por adivinhar uma geometria inicial e calcular o escoamento em tomo dela. As diferenças entre o escoamento calculado e o desejado à partida são usadas para modificar a geometria inicial. Em seguida, procede-se a novo cálculo do escoamento e entra-se num ciclo de iterações, até se obter convergência da geometria.

Estas três metodologias têm sido empregues em conjunção com os chamados M~todos de Progressão no Tempo. Estes últimos procuram obter a solução do cálculo do escoamento como o estado final, esta­cionário, a que se chega por intermédio de uma evolução das variáveis que definem o escoamento no tempo, evoluç~o esta regida pelas equações de Euler para regime transiente.

No presente artigo descreve-se um método inverso baseado no uso iterativo de um programa de análise do tipo Progres~o no Tempo. A especificação de desenho para o escoan1entQ.consiste na imposição ao longo da distância axial da distribuiçllo da velocidade tangencial média, V y , (a qual é definida como uma média ponderada pela massa, ao longo de linhas x =constante) e na distribuição desejada da espessura normal à linha de esqueleto.

A motivaçllo para o uso deste tipo de especificaçllo reside no fato de o valor da velocidade média tangencial depender fortemente do valor do ângulo local da linha de esqueleto, de taJ forma que será fácil desenvolver um algoritmo automático de alteração da geometria da linha de esqueleto. É de esperar que este algoritmo dê lugar a uma convergência rápida, devido à forte interdependência existente entre o valor da velocidade tangencial média e o ângulo local da linha esqueleto.

A cspecificaçllo da distribuiçllo da velocidade tangencial média pode soar algo estranho a principio, mas na realidade já foi utilizada anteriormente, como por exemplo em Novak e Haymann­Haber, 1983; Hawthome et ai., 1984; Pereira, 1993. Como é claramente mostrado em Novak e Baymann-Haber, 1983, a derivada da média ponderada pela massa da velocidade tangencial em ordem à distância axial é proporcional à carga aerodinâmica das pãs, definida como sendo a diferença entre os valores da pres~o entre o intradorso da pá. Atendendo a este fato, pode afirmar-se que sob o ponto de vista das especificações de projeto, o presente procedimento é equivalente aos métodos que especificam a carga aerodinâmica e a espessura.

Método Inverso

O prese nte método inverso aproveita a relaç!o existente, por um lado, entre o valor médio da velocidade tangencial ao longo de um passo, ponderado pela massa, Fig. I ,

JY,, pVxVydy

Y., ( I)

e a distribuição da carga aerodinâmica sobre as pàs e, por outro (.ado, entre os valores do ângulo de inclinaçllo da linha de esqueleto c a velocidade tangencial média, V y. A primeira relaçllo é utilizada para especificar a distribuição, ao longo da corda, da velocidade tangencial média que se pretende seja introduzida pelas pás, enquanto que a segunda é usada no algoritmo de alteração da geometria para atualizar a geometria da pá.

Efetuando um balanço de quantidade de movimento segundo a direção tangencial y, a um volume de controle do tipo [ 1234] indicado na Fig. I , facilmente se conclui (ver Novak e Haymann-Haber, 1983, ou Pereira, 1993) que a derivada da velocidade tangencial média em relação à coordenada axial x, é proporcional à carga aerodinâmica sobre as pás. Pode-se dizer que a diferença de pressão entre o extradorso e o intradorso das pás dá origem a uma força tangencial sobre o fluido que é responsável

397 J . e . B. T. Borges et at: Projeto de C.acatas de Pu para Turbina Baseado num Mêtodo ...

pela variação da velocidade tangencial mtdia Vy. Isto mostra que os valores da derivada de Vy devem ser nulos quer a montante, que a jusante das pás, por não serem exercidas forças segundo a direçllo tangencial nesses locais, devido à condição de periodicidade do escoamento.

"

Fig. 1 Casçata de Pb Tlplca e Linha de Esqueleto do Perfil.

Na posse destes resul!ados, nllo é dificil escolher uma evolução apropriada para a evolução de Vy ao longo da cascata. De fato, se o objetivo for evitar variações indesejáveis da distribuição de pres.sAo na superflcie das pás, será conveniente distribuir uniformemente ao longo da corda a força efetuada pela pá, o que se consegue teQ_do ~a carga aerodinâmica constante, a qual, no nosso caso, se traduz numa derivada constante de V y. E óbvio que a carga aerodinâmica e, consequentemente, a derivada da velocidade tangencial média devem apresentar valores nulos nos bordos de ataque e de fuga das pás, tal como acontece a montante e a jusante da cascata. Assim, devem existir regiQes de transiçllo na vizinhança dos bordos de ataque e de fuga das pás, onde os valores da derivada de V y, variam de zero até ao valor constante desejado.

A distribuiçllo da velocidade tangencial média utilizada no presente trabaJho foi escolhida de acordo com o exposto anteriormente, tomando em consideração uma dificuldade adicional originada por, na presente implementação, se ter definido a espessura do perfil como sendo medida na direçAo normal à sua linha de esqueleto. Consequentemente, exceto se a pá tiver uma línha de esqueleto axial no bordo de ataque, a geometria que resulta da adição da espessura na direçlio normal à linha de esqueleto apresenta, em geral, uma pequena regillo que se estende para montante do próprio bordo de ataque. Assim, esta regillo está, normalmente, sujeita a uma pequena carga aerodinâmica que dá origem a valores nllo nulos da derivada da velocidade tangencial média junto do bordo de ataque. Outra conseqüência deste fato é que o valor da própria velocidade tangenciaJ média no bordo de ataque é, em geral, um pouco diferente do valor constante verificado a montante da cascata. O mesmo se aplica junto do bordo de fuga, com as necessárias alterações.

Tomando em consideração o exposto anteriormente, a especificação da velocidade tangencial média foi efetuada usando as seguintes expressões:

~ ~ K{r1 +(1 - r1>t(2-t)}. o~çsç, dV ~ = K, Ç1 sçsç2 (2)

~ ~ K{l + (r2- l) (Ç - Ç2)

2

}. !;2

sÇ S I dÇ (I- Ç2) 2

J. of lhe IQ.z. Soe. Me<:Nnical Sc:ieflces- Vol. 18. Deoember 1996 398

onde Ç é a coordenada axial adimensionalisada pela corda das pás medida na direçilo axial, e rl'r2,Ç1 e 1;.2 são parâmetros cujos valores são escolhidos pelo p[.S)jetista {com O~ r1 < 1 e O~ r2 < I ). A constante K pode ser calculada a partir dos valores Vy nos bordos de

ataque e de fuga_

A transição entre os valores especificados para a derivada de V y nos bordos de ataque e de fuga, r1K e r2K e o valor constante na parte central da pá, é efetuada ajustando duas parâbolas com derivada nula. respectivamente, em Ç = Ç 1 e Ç - 1;2 .

A escolha da velocidade tangencial média como parâmetro de projeto justifica-se, não só por permitir controlar de uma forma direta a carga aerodinâmica dos perfis, mas também por facilitar a alteração da geometria do perfil das pás, entre duas iterações. Esta atualização é realizada com base na diferença entre os valores especificados para a velocidade tangencial média e os valores calculados com a geometria da pá na iteração anterior. A expressllo utilizada é a seguinte:

tg9 n+ 1 = RF{(\\{x)) - [(Vt{x))n - tg9 "l} + ( 1- RF)tg9 ~ (V,.) 1 sp (Vx) 1 J (3)

onde RF é um fator de relaxação (com valores entre 0,3 e 0,6), 9 é o ângulo da linha de esqueleto (Fig. I}, o lndice sp indica os valores especificados, o lndice I refere os valores muito a montante e os

lndices n e n+ I dizem (~)o ao núm(e; (~:)i~eração.

Atendendo a que <v.> 1 sp e ~ determinam os ângulos do escoamento médio, respecti-

vamente, nas geometrias analisadas e naquela que se pretende desenhar, facilmente se conclui que a Eq. (3) se baseia na hipótese de que a diferença entre os ângulos locais do escoamento médio e da linha de esqueleto é a mesma, que para as geometrias analisadas, quer para as geometrias a serem desenhadas. Em geral, esta é uma boa aproximação mas, não sendo exata, são necessárias algumas iterações para obter a geometria final da linha de esqueleto. Depois de conhecer os valores do ângulo da linha de esqueleto, a coordenada y da linha de esqueleto das pás pode ser finalmente calculada integrando numeric~ente a equação dy c/ dx = tg ( 9 { x)) .

Resultados O método de cálculo apresentado na seção anterior foi testado tomando como exemplo o projeto

duma coroa de pás de um rotor, com ftngulos de entrada e saída de valor elevado, para um escoamento a baixo número de Mach. A relaçAo passo-corda da cascata tem um valor de O, 7 e o número de Mach do escoamento à salda é 0,4, cm condições isentrópicas. As condições de projeto sAo p01 = 2, O MPa, T01 = 340K, p

2 = I, 791 MPa, 13

1 .. 65°e o lngulo pretendido para o escoamento à salda é

132 = - 70°. A espessura do perfil foi definida usando uma distribuição de espessura convencional em projetas de cascatas deste tipo (veja-se Pereira 1993).

2.5 2

1.5 1

0.5 1t .... o I> -o.5

-1 -1.5

-2 -2_5

-3

o &,e clncM.t

--IDldal

--F1•••

o 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9

Fig. 2 Velocidade Tangencial Mtdla Eapeelfteada, lnlei.I t Final, em Funçlo da Coordenada Axial.

399 J . E. B. T Borges et ai.: Projeto de cascatas de Pés pa.ra Turbina Baseado num Método ...

O programa de análise incorporado no método de desenho das pás utiliza uma formulaç!lo explicita. baseada no método da progressão no tempo, e uma discretizaç4o em volumes finitos, com uma malha estruturada gerada algebricamente e linhas quase-ortogonais curvas. Os cálculos foram efetuados usando uma malha com 25xl41 pontos e quatro nlveis de malha múltipla ( lx l , 2x2, 4x4, 8x8). O método convergiu ao fim de 15 iterações impondo um valor máximo, inferior a 0,025%, para a variação relativa da coordenada y, entre duas iterações. Como resultado dos cálculos obteve-se um ângulo do escoamento à salda p2 = - 70, 2° e um valor médio da queda da pressão de estagnação através da cascata que representa cerca de 0,2% do valor da pressão de estagnação à entrada. Estas diferenças são originadas pelos erros numéricos inerentes ao método da progressão no tempo, e slo indicativas da precisão com que é efetuado o cálculo do escoamento. Como resultado, após a convergência do método de projeto, o valor obtido para o número de Mach à saída, MBc2 = 0,398, , é ligeiramente menor do que o obtido considerando o escoamento como sendo isentrópico.

o -1

Fspeciflcada o -2 --Inicial -3 --Flnal

u..r -4 "O ·- -5 1>>. "O -6

-7

-8

-9 -lO

o 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9

Fig. l Oer1vada da Velocidade Tangencial JMdla Ea.,.eme.cta, lnlelala Final, em Funçao da Coordenada Axial.

Na Figura 2 mostram-se os valores especificados, iniciais e finais da distribuiç!lo da velocidade t~encial média em funç!lo da coordenada adimensional Ç. Aqui até o fim desta seç!lo, Vy representa a velocidade tangencial média adimensionada pela componente axial da velocidade_d.o

escoamento muito a montante das pás. Estes resultados indicam que a distribuição desejada para Vy foi obtida no final do processo iterativo com uma boa precisão. Os resultados de um teste mais exigente, sob este ponto de vista. estão apresentados na Fig. 3, na qual são apresentadas as derivadas das distribuições de velocidade média especificada, inicial e final, em ordem à coordenada axial, em -· função de Ç. Esta figura revela algumas discrepâncias entre a distribuição imposta para dV /dÇ e a obtida, após a convergência do método, na proximidade dos bordos de ataque e de fuga dos perfis. Tal como já foi referido. isto deve-se às pequenas cargas aerodinâmicas existentes nestas regiões, em resultado da adoçâo da detiniç!lo da espessura das pás na direçAo normal à linha de esqueleto dos perfis. Esta figura mostra igual~ente que entre 7 e 93 por cento da corda da pá especificou-se um valor constante para a derivada de Vy , enquanto que no bordo de fuga se impôs um valor nulo para a carga aerodinâmica. Uma vez que p1 tem um valor muito elevado, verificou-se ser necessário especificar um valor não nulo para a carga aerodinâmica no bordo de ataque. Nas regiões de transiçllo adotou-se uma evolução parabólica tomando para a Eq. (2), r1 = 0,5, r2 =O. 1; 1 = O, 07 e Ç2 = O, 93 .

J . ol lhe Bf82. Soe. Mechanical Scienoes - Vd.. 18, December 1996 ~00

- - Inicial

Fig . • G.ometrfa Inicial a Final da c .. cata

0.5

0.4

0.3

0.2 (.) - 0.1 -a ......

o .0.1 -- Inicial

.0.2 --Anal

.0.3 .Q.l o 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1.1

ç Flg.S Linha de Esqueleto Inicial e Final

A Figura 4 ilustra a alteração de geometria sofrida pelos perfis durante o processo iterativo. O perfil inicial foi gerado supondo urna variaçilo linear para a derivada da coordenada y da linha de esqueleto entre tg ( f1 ~ 1) e tg (f2~2 ) , com f 1 = I, 06 e f2 = I, O. Como se pode ver, a geometria obtida após a convergência do método é bastante diferente da geometria de partida. Este exemplo ilustra a capacidade do método para modificar substanciahUJWte a geometria da linha de esqueleto do perfil de fonna a aproximar a distribuição calculada de Vy da especificada. No entanto, a Fig. 5 mostra que o método de projeto nl\o modificou apreciavelmente a linha de esqueleto da pá nos primeiros I 00~ da corda. Este resultado indica que, neste caso, a distribuição de pressão junto ao bordo de ataque é uma função acentuada da distribuição de pressão e, somente em muito menor grau, do ângulo da linha de esqueleto.

401 J . E. B. T. Borges et ai.: Projeto de Cascatas de Pés para Turbina Baseado num Método .. .

2

o -2

-4 Q.

-6 u -8

-10 - -Inicial

-12 --Fmal

-14 ..Q. 1 o 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1.1

. Fig. 6 Dlatrtbulçlo de Preulo em Tomo daa P ..

A Figura 6 mostra as distribuições de pressão inicial e final, sobre a superflcie das pás. O coeficiente de pressão representado na figura foi normalizado pela energia cinética média na cascata, calculada com base nos valores médios das componentes da velocidade dos escoamentos nas regiões muito afetadas, a montante e a jusante, da cascata. Comparando esta figura com as as Figs. 2 e 3, toma­se evidente que se obteve uma car&!! aerodinllmica constante para 1;.2 >I;> Çl' região onde se especificou um valor constante para d v• I dÇ . Confirma-se assim que a carga aerodinâmica sobre os perfis apresenta uma evolução semelhante à de dV• / dÇ, tal como foi referido na seção anterior.

0 .9

0.8

0.7

--Inicial

--Final

0 .6 ~

0.5 ~ 0.4

0.3

0 .2

0.1

o -0.1 o 0.1 0.2 0.3 0.4 0.5 0.6 o. 7 0.8 0.9 1 1.1

Fig. 7 Dletribulçlo do NOmero de Mach na Superfle .. da P6

A Figura 1 compara as distribuições inicial e final do número de Mach na superficie da pá. Com exceçào das regiões muito próximas dos bordos de ataque e de fuga. pode ver-se na figura que a distribuição do número de Mach apresenta uma evo!llsão muito suave e quase-monotônica em resultado da especificação de um valor constante para dV /dÇ na região central da pá.

J. ot lhe Braz. Soe:. Mechanical Scieocea- Vol 18, December 19116 402

Sumário e Conclusões

O presente trabalho descreve uw método inverso que usa como dados de entrada uma distribuição da velocidade tangencial média, Vy, juntamente com a imposição da espessura normal da pá. Este método, que foi desenvolvido para o projeto de cascata de pâs bi-dimensionais, pode tratar pás com cargas aerodinâmicas elevadas e é baseado no uso interativo de um programa de análise.

Uma vez que a especificação da velocidade tangencial média ( Vy) não tem sido utilizada com freqUência, indicou-se a equivalência desta especificação com a imposição da carga aerodinâmica e mostrou-se que o seu uso é vantajoso no processo de alterar a geometria das pâs. Esta alteraçl!.o foi levada a cabo modificando a linha de esqueleto da pá, seguida da adiçiio da espessura na direçiio normal à linha de esqueleto.

A metodologia inversa proposta foi aplicada ao desenho de urna cascata de pás para turbina com elevada carga aerodinâmica, com o intuito de mostrar as suas capacidades. Através deste exemruo, mostrou-se que a evolução da carga aerodinâmica é semelhante à variaçiio da derivada de Vy. confirmando assim a afirmaçao feita acerca da equivalência entre estas _duas quantidades. Esta equvalência pode ser utilizada pelo projetista para a escolha da evolução de V y, permitindo um certo controle sobre a distribuição de pressiio nas superflcies da pá. Verificou-se também que esta distribuição de pressões é afetada apreciavelmente pela especificação de espessura normal, exigindo assim que a escolha desta espessura seja feita com cuidado a fim de evitar flutuações indesejáveis na pressão. A escolha apropriada da espessura requer alguma experiência por parte do projetista, o que constitui uma desvantagem deste método. embora esta fato seja comum a todos os métodos que especificam a carga aerodinâmica e a espessura da pá.

Referências

Dulikravich, G. S., 1992, "Aerodynamic Shape Design and Optimization: Status and Trends", JoumaJ of Aircraft. Vol. 29, n• 6, pp. 1020-1027.

Dulikravich, G. S. (Ed.), 1991, "Procoedings oflhe 3rd lntemationaJ Confercnce on lnverse Design Concepts and Optimization in Engineering Sciences (ICIOES-111)", Washington O. C.

Hawthome, W. R. et 11 .• 1984, '"Theory ofBiade Oesign for Largc Oeflections: Part 1- Two-dimensiooal Cascade", ASME Joumal of Englneering for Gas Turbines and Powcr, Vol. I 06, pp. 346-353.

Novak, R A. e Haymmann-Haber, G., 1983, "A Mixed-flow Cascade Passagc Oesign Procedure Based on a Power Serics Expansion", ASME Joumal of Engineering for Power, Vol. I 05, pp. 231-242.

Pereira, R. M. R. J., 1993, "Cálculo do Escoamento Compreenslvcl em Cascata de Pás", Tese de mestrado, Instituto Superior Ttcnico, Universidade T6cnica de Lisboa, Lisboa, Portugal.

Slooff, J. W. (Ed.), 1989, " Proceedings of AGARD Spccialists' Meeting on Computacional Methods for Aerodynamic Oesign (lnvcrse) and Optimization", AGARD CP-463 U\en, Norway.

Vand deo Braembusschc, R. (Ed.), 1990, "Procecdings ofa Special Course on lnvcrse Mctbods for Airfoil Oesign for AcronauticaJ and Turbomachinery Applications", AGARD-R-780, Rhode-St.-Gcnese, Bclgium.

RBCM • J. of the Bru. Soe. Mechank:al Sclencee Vol. XVIII- No.4 • 1g96

Abstracts

ISSN 0100.7386 Prinled ln Brazíl

da Cruz, J . J . and Matuoka, C. A. , "Posltlon Control of Mechanical Manipulators Uslng the Optimal Root Locus", RBCM- J. of the Braz. Soe. Mechanical Sclences - Vol. 18, No. 4, pp 308-317.

The design of a robust joint independent position cootrol scheme for mechanical maoipulators using a poJe placement approach based on lhe Optimal Root Locus is discussed in lhis paper. l ts performance is oompared to that oflhe traditional PD stnJcture, widely used in industrial applications, ín situations where lhe last one is not recommended. The influence ofvaríous feedforward schemes is also analysed. Keywords: Position Control, Mechanical Maoipulators, PoJe Placement, Optimal Root Locus, Linear Quadratic Regulator.

Barbosa, R. S. and Neto, A. C. "Rallway Wheelset Oynamics", RBCM - J. of the Braz. Soe. Mechanical Sclences - Vol.18, No. 4, pp 318-329 (ln Portuguese). A velocity paramelrized railway wheelset model was produced for dynamic behavior analysis. Lateral wheelset excursion related to the track is the first vibration mode. ll was observed in lhe eigen-properties that the wave lenglh of lhis movement is approximately const.ant and dependent on lhe wheel conicity. Modal damping of tbe first mode is inversely proponional to velocity and becames negative for high speed, resulting on system instability. Atlow speed, real, distinct and overdamped second mode e igenvalues are invcrsely proportionalto velocity and strongly coupled with contact stiflhess. Dynamic wheelset behavior through a variable track trajetory may be observed during model simulation. Keywonts: R.ailway Wbeelset, Dynamic Behavior Analysis, Símulation

Balestleri, J . A. P. and Correia, P. 8 ., "Cogeneratlon System Deaign Optlmizatlon", RBCM - J . of the Braz. Soe. Mechanical Sciencea- Vol. 18, No. 4, pp 330-337. Cogencration systcm design deals wilh several para.meters in lhe synthesis phase, whcre not only a lhermal cyclc must be indicated but lhe general arrangement, type, capacity and number ofmachincs need to be defincd. This problem is not trivial because many parameters are considered as goals in lhe project. An optimizalion technique that considers costs and rcvenues, reliability, pollutant emissions and exergetic efficiency as goals to be reached in lhe synthesis phase of a cogeneration system dcsign process is presented. A discussion of appropriated values and lhe results for a pulp and paper plaot integration to a cogeneration system are sho"''ll in order to illustrate lhe proposed melhodology. Ktywonts: Cogeneration, Optimization, Multiple Prognuruning, Design, Efficient Solutions.

Yu, L. and Righeto A., M., "Recent Advances and Appllcation of Turbulence Depth­lntegrated Two-Equatlon Closure for Modeling Contamlnant Dlscharges", RBCM • J. of the Braz. Soe. Mechanlcal Sciencea- Vol. 18, No. 4, pp 338-354.

This paper presents a briefreview and up-to-date development oflhe deplh·integrated turbulence two-equation closure. to be specilic. lhe depth-integrated I( - emodel and E- w model. A numerical simulation ofthe side discharge of waste heat into natural waters h as been taken as ao example of the practical application of lhe K- w model. ln addition. lhe aulhors also state lhe existing problems and further development of lhe present

deplh-integrated turbulence two-equation models. The maio objective oflhe paper is to integrate the research on the improvement and development of the present depth·integrated turbuJence models wilh lhe practice of industrial, environmentallllld sanitary cngineering. Ktywonts: Turbulence, Deplh-lntegrated Two-Equation Closure, Conwninant Discharge Modeling.

Esperança, P. T. T. and Sphaler, S. H., "Added Maas and Damping on Bottom-Mounted Rectangular Cylindera", RBCM- J . of the Braz. Soe. Mechanical Sciencea-Vol. 18, No. 4, PP 355-373. Look.ing for a major comprehension from the reason oflhe strong variations in lhe hydrodynamic coefficients of added mass and damping when a twQ·dimensional body oscillates venical close to lhe free surfaoe in finite depth we analyze this problem using eigcnfunction expansions. lt can make evident the parameters involved on resonances and lhe behavior for low and high frequencies. Jntroducing separation ofvariables lhe boundary value problem is transformed in an eigenvalue problem ofSturm·Liouville k.ind. This procedure was used by Yeung (1981) to determine the hydrodynamic coefficients for a truncated vcr1ical cylinder wilh circular section piercing lhe free surface oscillating in water of finite depth, Yeung and Sphaier ( 1989) to consider lhe interference of vertical walls of a canal in the hydrodynamic coefficients, Mel ver and Evans (1984) in lhe

J . of lhe Braz. Soe. Mechanical Sdencas - Vol. 18. December 1996

problem of a submerged vertical cylinder oscillating near the free surface and Esperança (1993) for the two­dimensional case. ln thc present paper we co.nf111Tl the occurrence of negative added mass when lhe body oscillates close to lhe surface, and obtain finite values ofthe added mass and damping coefficient in dimensional form for the zero frequency li mil as expected. The negative values of the added mass are related to two kind of resonant .wave modes occurrlng cJose to values oflhe frequency parameter, wave number times hal f rectangle widlh, m~a, is equal to n1t and nx/2. For these frequencies lhe damping coefficient is equal to zero or achieves a local maximum value, and stationary waves characterize the internal flow. Negative added mass for lhe frequency parametcr going to zero is obtained for shallow water, wbcn lhe emergence ofthe rectangle is larger lhan 113 of lhe water deplh. A shallow water approach is also presented, allowing us to describe the wave motion for small water deplh and clearance. The results derived from the shallow-water solution is in good agreement with the one of the complete solution. furtber we presenl lhe solution for lhe case of rectangular cylinder beaving on lhe free surfàce and compare lhe results wilh results presented by Bai and Yeung (I 974), Bai (1977) and by Drimer, Agnon and Stiassnie (1992). Keywords: Motion of Floating Bodies, Waves, Ship Design, Added Mass and Damping, Botton - Mounted Rectangular Cylioder.

Rade, 0 ., A. and Lallement. G., "Vibration Analysls of Structures Subjected to Boundary Condltion Modiflcatlons Using Experimental Data", RBCM - J . of the Braz. Soe. Mechanlcal Sciences - Vol. 18, No. 4, pp 374-382.

This paper addresses tbe Analysis of Modified Structures by using experimental data. ln particula r, modifications ofthe boundary conditions ofthc structure by grounding of ooc or severa! of its degrees of freedom are considered. Three mctbods are proposed, which are conceived to obtain the eigenvalues and eigeovectors of more constrained configurations, given lhe experimental Frequcncy Response Functions measured on a less constrained configuration. The formulations of lhe three melhods are first presented and their performances are lhen evaluated through applicatioos to an automotive structure tested in laboratory. Keywords: Structural Modifications, Modal Analysis, Anliresonances

Tsuzuki, M. S. G. and Moscato, L. A., "Use of Voronoi Diagram for Oeflning the Cutting Path", RBCM - J . of the Bra.z. Soe. Mecha nicai Sclences - Vol. 18, No. 4, pp 383-394 (ln Portuguese).

ln this woric the use ofVoronoi Diagram for dctining the cutting path is discussed. ln tbis fashion the cutting path is dctermioed in two steps: first, lhe Voronoi Diagram is created, and second, the cutting path is calculated based on lhe Voronoi Diagram. A possible implementat.ion for both steps is presented. The main purpose ofthis algorithm is to facilitate lhe comprehension of the evolved concepts. Results obtained from a prototype are presented. Keywords: Pocket Milling, Cutting Path, Voronoi Diagram

Borges, J . E. B. T., Gato, L. M. C. and Pereira, R M. R. J., "Oesign of Turbine BJade Rows Using a Tlme-March ing Method", RBCM - J. of the Braz. Soe. Mechanical Sciences - Vol. 18, No. 4, pp 395-402 (ln Portuguesa). An inverse based on lhe iterative use ofa direct (analysis) timc-marching code is described. The design input specification used consists of an imposition of the mean tangencial velocity distribution and the desired distribution of blade thiclrness normal to lhe blade camber given as a function of the axial distance. The desc ribed inverse method was applied to thc derange of a turbine cascade typical of a rotor blade row representing a high aerodynamic load. This example brings out clearly the relationship between the blade pressure loading an the derivation oflhe mean velocity distribution used as design input. Keyworcb: CFD, Turbomachine. lnvcrse Melhods, Turbine Cascade

Nota for Contributora: Artlcln on Olak • Authora are aneou,..gad to aubmlt tha final •ccapted manuacrtpta on dlak, ualng taxt editora for

Windows or F,..me Makar

• Tha dlak must ba martted wlth tha papar ldentlfleatlon numbar and softwa,.. uaed. TWo eoplu of tha prlntout ahould ba lneluded.

R8CM - J . of Ule Braz. Soe. Mechanlcal Sclences Vo/. XVIII - No.4- 1996

ISSN 0100-7386 Printed ln Brazil

JOURNAL OF THE BRAZILIAN SOCIETY OF MECHANICAL SCIENCES REVISTA BRASILEIRA DE CI~NCIAS MECANICAS VOLUME XVIII - 1996

Referees - 1996

Abelardo Alves de Queiroz Almir Monteiro Quites Alvaro Costa Neto Américo Scotti Amílcar Porto Pimenta Antonio Arlindo Guidettl Porto Antonio C. P. Brasil Júnior Armando Albertazzi Atair Rios Neto Átila Pantaleao da Silva Freire Aureo Campos Ferreira Benedito Di Giacomo Benedito Manoel Vieira Breno Pinheiro Jacob Carlos A. Flesch Carlos Alberto de Almeida Carlos Alberto Schneider Carlos Frederico Bremer Carlos Henrique Ahrens Celso K. Morooka César Philippi Celso Pupo Pesce Cláudio S. Kiminami Clóvis Raymundo Maliska Edgar Nobuo Mamiya Edison Hiroshi Tamai Edson Bazzo Edson Gomes Eduardo Lobo Lustosa Cabral Eduardo Morgado Belo Euclides Mesquita Neto Francisco Domingues Alves de Souza Hans lngo Weber Jair Carlos Outra João Andrade Carvalho Júnior João Batista de Aguiar José Arnaldo Barra Montevechi José Francisco Ribeiro José Guilherme Senna José Jaime da Cruz

UFSC UFSC

USP-EESC UFU

CTA-ITA UNICAMP-CT

UnB UFSC

INPE-LAC UFRJ-COPPE

UFSC USP-EESC

CTA-ITA UFRJ-COPPE

UFSC PUC-RJ

UFSC USP-EESC

UFSC UNICAMP - FEM

UFSC USP-EP

UFPB UFSC

UnB USP-EP

UFSC USP-EP USP-EP

USP-EESC UNICAMP-FEM

IPT-DME UNICAMP-FEM

UFSC INPE

USP-EP E FEl UFU

CTA-IEA USP-EP

J . ofthe Bnlz. Soe. Mechanical Scienc::a - Vol 18, Oecember 1996

José M. Saiz Jabardo José Roberto de França Arruda Kazuo Nishimoto lourival Boehs luis Carlos Sandoval Góes luis Felipe Mendes de Moura Marcelo J. S. de lemos Márcio luiz de Souza Santos Marcos Ribeiro Pereira Baretto Miguel Cezar Santoro Miguel H. Hirata Newton Ribeiro Rocha Nivaldo lemos Coppini Olivío Novaski Oswaldo Horikawa Paulo Eigi Miyagi Paulo Tadeu de Mello lourenção Ricardo de Andrade Medronho Rogério Tadeu da Silva Ferreira Sebastião Carlos da Costa Sérgio Coite Silvia Azucena Nebra de Perez Vilson Carlos da Silva Ferreira Walter lindolfo Weingaertner

USP-EESC UNICAMP-FEM

USP-EP UFSC

CTA-ITA UNICAMP

CTA-ITA UNICAMP-FEM

USP-EP USP-EP

UFRJ-COPPE FMG-DEMEQ

UNICAMP-FEM UNICAMP - FEM

USP-EP USP-EP

EMBRAER UFRJ-Escola de Qulmica

USFC E FEl

UFSC UNICAMP-FEM

UFRGS UFSC

RBCM - J . of the Bru.. Soe. Mechanlcal SclencH Vol. XVIII - No.4- 1996

ISSN 0100-7386 Prlnted ln Brazil

JOURNAL OF THE BRAZILIAN SOCIETY OF MECHANICAL SCIENCES REVISTA BRASILEIRA DE CI~NCIAS MECÂNICAS VOLUME XVIII - 1996

Author Index

A

AI-Qureshi, Hazin Ali- Plastic lnstability of Anisotropic Thin Plates in the Biaxial Stress State (ln Portuguese). No. 2, pp. 174-181. Almeida, Sérgio F. Muller de- Pre1iminary Considerations on a New Design Concept for Be1ow Knee Energy Storing Prosthesis. No. I, pp. 1-7. Alvaro Toubes Prau- Analysis ofthe HFC-134a Refrigerant Flow Through Capillary Tubes Using the Two-Fiuid Model (ln Portuguese). No. I , pp. 50-58. Arruda, José Roberto de França -A Data Compression Method for the Modal Analysis of Spatially Dense Laser Vibrometer Measurements. No. I, pp. 17-23. Augusto, Oscar Brito - About Shear Stress in Ship Structures. No. 2, pp. 163-173. Azevedo, Joio Luiz Filgueiras de - The Wake Behind a Backward-Facing Step ·A Numerical Study. No. 3, pp.248-262.

B

Baek, Nelson- Analysis ofthe Structural Static Response, Eigenvalues and Eigenvectors ofManipula· tor With Flexible Joints and Links (ln Portuguese). No. I, pp. 40-49. Balestieri, José Antonio Perrella - Cogeneration System Design Optimization. No. 4, pp. 330-337. Barbosa, Roberto Splnola- Longitudinal Train Dynamics. No. 2, pp. l07-116. Barbosa. Roberto Splnola- Railway Wheelset Dynamics (ln Portuguese) .. No. 4, pp. 318-329. Berni, Mauro Donizeti • Transport and C02 Emission: A Multiobjective Ana1ysis for Energy and Environmental P1anning (ln Portuguese ). No. 2, pp. 117-126. Borges, Joio Eduardo de Barros Teixeira- Design ofTurbine 81ade Rows Using a Time-Marching Method (ln Portuguese). No. 4, pp. 395-402. Braga, Arthur M. B. - Approximate Theories for Vibration of Laminated PI ates. No. 2, pp. 198-206. Button, Sirgio Tonini- Cylinders and Rings Hot Upset Forging Analysis by the Finite-E1ement Method (ln Portuguese). No. I , pp. 24-32.

c Correia, Paulo de Barros -Transport and C02 Emission: A Multiobjective Analysis for Energy and Environmental P1anning ( ln Portuguese). No. 2. pp. 11 7-126. Correia, Paulo de Barros- Cogeneration System Design Optimization. No. 4, pp. 330-337 . Costa Neto, Álvaro- Railway Whee1set Dynamics (ln Portuguese) .. No. 4, pp. 318-329. Cruz, José J aime - Position Control of Mechanical Manipulators Using the Optimal Root Locus. No. 4, pp. 308-317.

D

Diniz:., Anselmo Eduardo - Monitoring the Tool Wear in the Tuming Process Using Acoustic Emis­sion (ln Portuguese). No. 3, pp. 227-238. Duduch, Jaime Gilberto - Model of Brittle Materiais Single Point Machining With High Remova! Rates. No. I, pp. 33-39.

J . of the Btaz. Soe. Mechanical Sciences • Vol. 18, Deoember 1998

E 408

Esperança, Paulo de Tarso T.- Added Mass and Damping on Botton-Mounted Rectangular Cylin­ders. No. 4, pp. 355-373.

F

Ferreira, Maria do Carmo- Fluid Dynamics Bchavior of a Pneumatic Bed With a Spouted Bed Type Solid Feeding Systcm. No. I , pp. 67-73. Forcellinl, Fernando Antonio - Analysis of lhe Structural Static Response, Eigenvalues and Eigen­vectors of Manipulator Wilb Flexible Joints and Links (ln Portuguese). No. I, pp. 40-49. França, G. A. C. - Turbulent lncompressible Flow Within a Channel Wilh Transpiration: Solution of lhe Coupled Problem Porous Waii-Main Flow. No. I , pp. 88-96. Freire, Atila P. Silva- On Kaplun Limits and lhe Asymptotic Structure of lhe Turbulent Boundary Layer. No. I , pp. 80-87. Freire, José Telleira - Fluid Dynamics Behavior of a Pneumatic Bed With a Spouted Bed Type Solid Feeding System. No. I , pp. 67-73. Frey, Sérgio- On lhe Numerical Heat Transfer Based Upon Mixture Theory Space Dynamics. No. 3, pp. 282-296. Frota, Maurkio N.- Wall Effects for a Sphere Falling in a Non-Newtonian Fluid. No. I, pp. 97-102.

G

Gama, Rogérío M. Saldanha da • On the Numerical Heat Transfer Based Upon Mixture Theory Space Dynamics. No. 3, pp. 282-296. Gato, Luis Manuel de Carvalho - Design ofTurbine Bla.de Rows Using a Time-Marching Method (ln Portuguese). No. 4, pp. 395-402. Gee, Anthony E. - Model of Brittle Materiais Single Point Machining With High Remova! Rates. No. I , pp. 33-39. Gomes, Marcos Sebastilo de Paula- lnertial EfTects on the Retention ofParticles in the Near Wake of Blunt Obstacles. No. I, pp. 74-79.

llkiu, Anselmo Monteiro- Plastic lnstability of Anisotropic Thin Plates in the Biaxial Stress State (ln Portuguese). No. 2, pp. 174-181.

J

Jasinevicius, Renato Goulart- Model of Brittle Materiais Single Point Machining With High Remova! Rates. No. 1. pp. 33-39.

K

Knudsen, J. Claudio da C. M. - Analysis of Approxirnations to lhe Delay Dynamics in the Optimiza· tion ofControl Systems (ln Portuguese). No. 2, pp. 182-197.

L

Lallemand, A.- Turbulent lncompressible Flow Wilhin a Channel Wilh Transpiration: Solution ofthe Coupled Problem Porous WalJ-Main Flow. No. I, pp.88-96. Lallement, Girard - Added Mass and Damping on Bonon-Mounted Rectangular Cylinders. No. 4, pp. 374-382. Lima, Washington J. N. de- Approximate Theories for Vibration ofLaminated Plates. No. 2, pp. 198-206.

409

M

Marczak. Ro&~rio J.- Analyticallntegration of Membrane-Bending Coupling Tenns in Non-Linear Boundary Element Analysis of Mindlin-Reisner Plates. No. 3, pp.218-226. Martins, Fábio -Use of Modified lmplant Testas a Manner of Detennirung Susceptibility to Reheat Cracking (ln Portuguese). No. 3, pp. 273-281. Martins-Costa, M.lria Laura - On the Numerical Heat Transfer Based Upon Mixture Theory Space Oynamics. No. 3, pp. 282-296. Matuob, Carlos Alberto- Position Control of Mechanical Manipulators Using the Optimal Root Locus. No. 4, pp. 308-3 17. McConnel, Kennetb G. - Oynamic Response of a Damped Cable with a Compliant Oamped Support System. No. I, pp. 8-16. Medeiros, Carlos Antonio de- Selection and Assignement ofMachines (ln Portuguese). No. 3, pp. 239-247. Melo, Claudio- Analysis ofthe HFC-1 34a Refrigerant Flow Through Capillary Tubes Using lhe Two-Fluid Model (ln Portuguese). No. I, pp. 50-58. Milioli, Fernando Eduardo - Atmospheric Bubbling Fluidi.zed Bed Combustion: Application to Hig.h Ash Coais and Approach to Scientific Research. No. 2, pp. 127-142. Morgenstern Jr., Algacyr- Supersonic Flow Past Pressure Vent Orifices on Satellite Launcher V chi­cles. No. I , pp. 59-66. Morooka, Celso Kuuyuki- Ocean Waves: Simulation and Spectral Analysis. No. 2, pp. 143-162. Moscato, Lucas Antonio- Use of the Voronoi Diagram for Oefining lhe Cutting Path (ln Portu­guese). No. 4, pp. 383-394.

N

Nicoletti, Rodrigo - Self-Excited Vibrations in Active Hydrodynamic Bearings. No. 3, pp. 263-272. Nogueira, Felipe- Preliminary Considerations on a New Design Concept for Below Knee Energy Storing Prosthesis. No. I, pp. 1-7.

o Oliveira, Mareio Mattos Borges de- Selection and Assignement of Machines (ln Portuguese). No. 3, pp. 239-247. Ortega, Marco• Aur~lio- The Wake Behind a Backward-Facing Step- A Numerical Study. No. 3, pp.248-262.

p

Pereira, Rui Manuel Roma de Jesus- Design ofTurbine Blade Rows Using a Time-Marching Method (ln Portuguese). No. 4, pp. 395-402. Pigari, Almir- Monitoring the Tool Wear in the Tuming Process Using Acoustic Emission (ln Portu­guesc). No. 3, pp. 227-238. Porto, Arthur José Vieira - Model of Brittle Materiais Single Point Machining With High Remova! Rates. No. I, pp. 33-39. Prado, Antonio Fernando Bertachini de Almeida- Optimal Rendezvous Maneuveurs for Space V chicles. No. 3, pp. 297-301

R

Ribeiro, Cassilda Maria- Selection and Assignement ofMachines (ln Portuguesc). No. 3, pp. 239-247. Ribeiro, José Francisco Ferreira - Selection and Assignement of Machines (ln Portuguesc). No. 3, pp. 239-247. Rade, Domingos Alves - Analysis of Structures Subjected to Boundary Condition Modification

J. of the Braz. Soe. Mechanieal Scienoes • Vol. 18, Oeoember 1996 -410

Using Experimental Data. No. 4, pp. 374-382. Rigbetto, Antonio Marozzi- Recent Advances and Application ofTurbulence Depth-lntegrated Two­Equation Closure for Modeling Contaminant Discharges. No. 4, pp. 338-354.

s Santos, limar Ferreira- Self-Excited Víbrations in Active Hydrodynamic Beanngs. No. 3, pp. 263-272. Santos, Luiz Antonio S. B. -A Data Compression Method for the Modal Analysis of Spatially Dense Laser Vibrometer Measurements. No. I, pp. 17-23. Seixlatk, Andre Luiz- Analysis ofthe HFC-134a Refrigerant Flow Through Capillary Tubes Using the Two-Fluid Model (ln Portuguese). No. I , pp. 50-58. Silva, Elcio M. V.- Fluid Dynamics Behavior of a Pneumatic Bed With a Spouted Bed Type Solid Feeding System. No. I, pp. 67-73. Silva, Evaodro Cardozo da- Cylínders and Ríngs Hot Upset Forgíng Analysis by the Finite-Eiement Method (ln Portuguese) . No. I , pp. 24-32. Sotelo Júnior, Josl- Analysís of Approxímatíons to the Delay Dynamícs ín tbe Optimizatioo of Coo­trol Systems (ln Portuguese ). No. 2, pp. 182-197. Sphaier, Sérgio Hamilton - Added Mass and Damping on Botton-Mounted Rectangular Cylinders. No. 4, pp. 355-373.

T

Tedescbi, G.- Turbuleot Incomprcssíble Flow Within a Channel With Transpíration: Solutíon of the Coupled Problem Porous Waii-Main Flow. No. I, pp.88-96. Teixeira, Marco André O. M.- WaJI Effects for a Sphere FaJiing in a Non-Newtonian Fluid. No. I, pp. 97-102. Teodoro, Elias Bittencourt- Dynamíc Rcsponse of a Damped Cable wíth a Complíant Damped Sup­port System. No. I, pp. 8-16. Trevisan, Roseana da Eultaçlo- Use of Modífied lmplant Testas a Manner of Determining Suscep­tibility to Reheat Cracking (ln Portuguese). No. 3, pp. 273-281 . Tsuzuki, Marcos S.G.- Use ofthe Voronoi Diagram for Defining the Cutting Path (ln Portuguese). No. 4, pp. 383-394.

w Weber, lbns lngo- Longitudinal Train Dynamics. No. 2, pp.107-116.

y

Yokoo,lrineu Hlroshi- Ocean Waves: Simulation and Spectral Analysis. No. 2, pp. 143-162. Yu, Liren- Recent Advances and Application ofTurbulence Depth-lntegrated Two-Equation Closure for Modeling Contaminant Discharges. No. 4, pp. 338-354.

RBCM. J. ofthe Bru. Soe. Mechanleal Sclences Vol. XVIII· No.4 • 1996 • pp. 411-414

FORMULÁRIO DE AFILIAÇÃO

ASSOCIAÇÃO BRASILEIRA OE CiêNCIAS MECÂNICAS Av. Rio Branco, 124 -18" andar· 20040..001 Rio de Janeiro· RJ Tet.: (021) 221-0438 ·Fax.: (021) 222·7128 e-mail: abcmalfs4» omega.lncc.br CGC 83.431.593/0001-78

INDIVIDUAL

Por favor, preencha os dois lados do formulário.

ISSN 0100-7386 Printed ín Brazil

Nome __________________________________________________________ ___

Endereço Residencia1----------------------------------------------------

CEP --------------------- Cidade Estado Pais ----------- Tel.: ( ) ______ Fax: ( ) __________ ___ E-mail _________________________________________________________ _

Empre~------------------------------------------------------------Dept./Divisào ____________________________________ Posição ----------

Endereço Comercial --------------------------------------------------CEP --------------------- Cidade Estado País Tel.: ( ) _______ Fax: ( ) _____ _

E-mail-------------------------------------------------------------

Candidato-me a: Na categoria de:

O Admissão O Sócio efetivo

Solicito enviar correspondência para o seguinte endereço:

O Comercial

O Mudança de Categoria O Sócio Aspirante

O Residencial

Data------------------------ Assinatura------------------------------

Para uso da ABCM

Aprovado-------------------- Data---------------- Sócio 0° -----------

J. ofthe Staz. Soe. Mechanlcal Sciences- Vol. 18, Oecember 1996

FORMAÇÃO E EXPERIÊNCIA PROFISSIONAL

Por favor, liste em ordem cronológica os dados completos de sua formação e experiência profissional A falta desses dados impedirá o processo de admissão. Obrigado.

FORMAÇÃO ACADÊMICA

Graduação- Área Anos a lnstituiç Pais

Mestrado- Área Anos a Instituição Pais

Doutorado- Área Anos a Instituição País

Outro- Área Anos 8

Instituição Pais

EXPERIÊNCIA PROFISSIONAL

Empres.a_ _____________ _ Anos a Natureza da atividadle----------- Posição

Anos a Posição

Anos 8

Posição

Indique atê um máximo de 8 áreas de acordo com os códigos numéricos do Anexo.

Áreas de Especializ.aç!l,o--------------------------

AplicaçãO-----------------------------

Comen~ios-------------------------------

J. ot lhe 8nlz.. SOe:. Mec:hanic:al Sciences • Vol. 18, Oec:ember 1996 413

Área de Especialização

Especifique no Fonnulário do Afiliação os códigos numéricos das suas Áreas de Especializaçao e de Aplicação (verso).

1000 Funct.nentoee MModoe B61Jlc:oe em M•dnlc:a Taóric:a • Apllc:8da

1010 MKAníca do Continuo 1 110 ~doo Elemenlo. F1nilos 1 120 ~doo Elemenlo. de ConlOmO 1130 "'*"'*'-~ 1 1«l ~du ~ Fl'lb.s 1 1 !iO Oúro M6lodoe em Mec:. Comll'Aecionel 1210 M6liOdos Ellac6alioot e E-lsticos 131 o Modelegem 1 ~ 10 Fundementoe de M61in ExperimenUII 1 51 o Metrologill 1 S 1 O a.tro::la de Projeloo 2000 Dln6mlc:e. ~ 21 10 Clrwn6tice • Oinama 2210 Vll:nçOes de Sólido.· Fundemenlos 2310 Vt>t8ç6eo ·E~ de E.ltnJluru 2320 lllbteç6es • Eoti\Jiunoo 2330 Propegaçloo de Onde• em Sólidoo 2340 l"'l*:lo em Sólldol 2350 Ide! lllfic8çto de Pnnetroe 2A20 Propegeç:to de Ond .. em F\Jdoo 2510 ~ F\Jiclo..€11n.Cun1 2610 Aatroniulíca • Mec:. Celelte e Ol1lit8l 2710 Explodo. Befiltlce 2810 kústicll 3000 Controle • Ollml.uçlo 3110 Projeto e Teorie de Sletemaa Mecllnicoo 3210 s~ de Conlrole Otlmo 3220 Silllemao de Controle Adoptatlvo 3230 ~em SmemM • Controle 331 o Robólíca 3A10 Otr'llizaç6o de Sillemu e ~'roc::eS$0~ 4000 .......... ~ 1 1 o 8lonlalerlelo A 1 20 Metetiolia Mdiooe ~ 130 Mater\111$ c.rtmiooe ~1«l Meterillio Pol~ ~1!10 ,........~

"'210 Co111a ro.-çAD Mec6nice OlO Cal.,..iz8çlo • Controle Míc:roestn.Aunll ... 10 Comp. Mec*lico doo ~ .,_20 Camp. Mec:. Mal. • 8ecxu TemperllutM ~ Camp. Mec. Mil. ·Aliei Temperlt\nl ~ Comp. Mec. Mil. • C..-regemlo. Veri6vll ~ Comp. Mec. Mil. • Cerregemto. O!Nrn. .SOO Mecaniomoa de F11111111 ~ Mednice de Frat ... ~710 EnMios 0-.Aivol 4720 EIIMioe Nlo-OIÀWYOI oiiiOO eor-ao 1000 Mec:tnlc:a dDe S611doe 5010 Eleltlc:icllodl LW.. !1020 Elelliciclecle Nlo-L"-5030 Vleooelaolidd8d 50«l Pluticidede 5050~ 50110 .......,. doo Mllerillo Corfuoedoa 5070 Mlc:*'ka doo Meóot Pcroeoe 5110 Reologja 5210 C.., Hatee. Vlgu 5220 MembriiMI, P- • c-5230 EltnJILno • Geral 52«l EltnJtLno • Cont.to com o Solo 5250 Elln.Cunla • SubmerMa/SemhsubmerNII 5260 EatnAunll • Móotete 5270 Eltnanl • Vuot 1 Contençio 5310 Mectnica doo Soioe • Búoc:o 5320 Mec6nk:a doo Soioe • Aplicaçi>et 5330 Medníca dei RochM

~ 1 o Efeitoe Elelr<>-Magnétic:oa em M. Sólido• ~20 Efeit001 Ténnicoa em M. doo Sólidos 5510 Eltabfliclede de Ea1rutura1 5520 CompclfÚimenlO apóe a Fllwnbagem 5530 Eatldoo lim~u 1 Cá1181 de Colapeo 55<10 Acamodeç6o • AcúmlAo de Oeno se 1 o lolec:6nica de Fr~~~~.n 54560 Trllolcgi8 5855 Atrito e Oeagute 5710 Componentes de MéquineS 5720 Ac:opllr'nenloo e Juntas N6().5oldad .. 5800 An*lile ExperimWilal de Tensões 1000 Mednic:a doe Fluldol 6010 RIOOiogla 811 o HicHuica 6210 E.c:oamento lncrom~lvel 6220 E~o Compnaufvet 8230 EICOIIY*lto RMifllto 82«l Eecoamento em Mlioo PCI'OOIOI 6250 Magnet().Hidrodln•mlea e Plasma 8270 EIQOemento MuHN••ico 8310 Cameda Umfte • Conlomo S6fido 5320 Cemeda Lm«a · Contorno LMe &4 to EICC8metlto lnlemo • Dulos, c-.. etc. &430 E~o oom Superftaellvte 8510 Eatabllidade de E1100emento 6520 TUitYI6nc:ia 6530 Hidrodlntmlc8 • Velc;ulo de EslNt Naval 66-40 Aarodin6mica 6610 Mie:. Fluidoe • ApliaeQ611 em M~onaa 6650 Lubrificaç.io 6710 T~ea em FlUido~ 6810 T6c:. E.xpettal e~ E~ 7000 Tlr'II'IOCitnciH 701 O T IITI!OdNmlc8 7110 Tranap. de Cllio< • Conv.c:. Monof'•ic8 7120 Tranap. de C81or • Convec:. Bifáolea 7130 T~. deCIIIot - Conduç6o 71«l T~. de Celor • Radlaçtor'MOd. Comb 7150 TFW'oll). de c• · Oi~iotamaa 7210 Termodintmíca doo Sólido. 7310 Trenopc~~W de Meue 7 ~ 10 CombuotAo 7 ~20 Combutllo em Le~o Fluidíz.ado 7510 Adonaclotu e Oiapaaitivo4 de PropAalo 1000 Geoc:ltnc:l• 601 O Mictomerítica 8110 Me;os Pcro.oo 8210 Geam-*'""' 8310 Mlc::*líca doo Ab11ot SiiUnÍOOI &410 ~ ~ Mel8on>logl8 1000 Enervla • Melo Alnblenta e110 Com~Mtivm Fôanll 9120 Slttamea Nueleetlla • Fiulo 9125 Sistemas Nuclurlo • Fuslo 9130 SIStemu Geol6rmi<U 91~ SlltemesSolatlll 9150 Slltemaa Eólico• 9160 Slotamu de en.g,. Qceln;ca 921 O AtTMz~Nmento de Enlfgia 9220 o .. lribuiçlo de Energo• 9310 Mec:iníca doo Fluldoe Nr1blental e. 1 O Mee6níca de Oitpo.Kivol de Armuenamento de Reslduo 1 OOOOIIIod6nciH 101108lomec:6nlc:a 1 0210Ergonomia 1031 ORNbilàçloo 1041 0Mec6nica noa Eopor1n

414 J . ol the Braz. Soe. Mechaolcal Sciences- Vol. 18, December 1996

Áreas de Aplicação

Exemplo: um especialista em Mecânica dos Fluidos (famllia 6000) atuando na área de Turbulência (6520), deverá escolher a Área de Aplicação 350, se estiver trabalhando em Propulsao.

01 0 AWIIice e Controle de Ruldo 020 ~~em BIDCiAnclu 030 CAD OoiO CAM 050 Componentes de Máquinas 060 c~ AmbientAl 070 ConltOie de Qualidede 080 CnogeNa 090 Enoeman• e Flaiea de Reatorea

1 00 Engert\llrill de Pelrclleo 110 Engemaria ~Q 120 E~01 de Proceuoa 130 Equlpemenro. lndualrllna 140 Fontea Allemallvea de E'*ll,. 150 Forjamento 160 Funcltçto 170 c.r..m. de Quelidade 180 IndUstrie T 6xt~ e T ecnologie Correlatll 190 lnapeçto e Cert!faoçAo

200 lnlteleçóes lnduslnaia 210 lnJin.mentaçAo 220 Lubnf~ lndul1rilll 230 Mancais e ROiementos 240 Máquinu Fetr..,..entas 250 M~ de Fluxo 260 M6qunas Moem• 270 Mdnica Fina 280 Metalurgia Geral e B<lnllfidwnento de Minério 290 Matrologia

300 Mt1WIOÇ6o e Met81urgoe Extr8\MI 310 ()ptoc:a 320 Pontes e B«regen. 330 ProcMsoe de Febricaçlo 340 Projem de EllruCuru 350 Prop.II$Ao 360 ProspeoçAo e PmpuiSio 370 Setvo MecaniSI'Ios e Controle 380 Sld8fl.l'gill 390 Sostemas HICHulicos

400 S11tem111 Pneumilllicos 41 o Soldegem 420 Sohcilaçóea AciOerKais • Efeotos de Vento. Slamo. EJ<p4odo. Fogo • II'IUldaçAo 430 TecnOlogia de Almentoo 440 T ecnologla M...,.., 450 Trenaporte (axckJido valculoo) 460 TrenamissAo de Energia 470 Tratamento T6nnlco e Tennoqulmlco 480 T~ lndullne>S a NuciMr• 490 Usinal~

500 Uainaa Termoa"lticas 510 VKuc> 520 Vuo. de P,.ado, Troc:adorea de Calor a Equipamentos P.-530 Vaioulos • T.,...nt. Eapacoait a Morillmoe

RBCM - J. of the Braz. Soe. Meehanleal Sclenc" Vol. XVIII - No. 4 - 111116

.......... Nadoullleat aad Mali ~nufer Coa,._ee ... nini ISHMTf&SM& Heat aad Maa Trusfel' C........_

P DeeembS29-31,1997 ~ rlnt A .. ou~He~M~tt ud Calfltr Papen ~-

lndian lnstitute ofTecbnology - Kanpur (lndia)

Organi:ud under the auspices of

lndian Society of Mecbanical Engineering

Co--sponsored by

SESI, fiChE, AESI, Comb lnst (I)

in coorperation wltb

ASME - American Society of Mecbanical Engineers

Deadlines:

15 December, 1996 - 3 copies of abstract due

IS J anuary, 1997 - Abstract Acceptance

IS April, 1997 - Fulllengtb manuscript due

IS August, 1997 - Acceptance after revision

IS September, 1997 - Final Mats due at lSHMT HQ

ISSN 0100-7388 Printed ln Brazn

Prospective authors should send three copies of lhe abstracts not exceeding 500 words to the coordjnating scientist in their zone. The work must be original and unpublished. The abstract should clearly state the objectives of lhe work, results and conclusions.

The covering letter should contain: (i) Fivc keywords to describe and categorize the work easily (ii) Name, address, phone number, fax number ande-mail address ofthe author to whom correspondence should be directed.

3rd Interoational 'rltfrmal Eaeify Conareu Kitakyusbu, Japaa - July l8-A~ut 1, 1997

Sponsored by Kyushu Universlty, City of Kitakyusbu a.nd

Fukuoka Science & Tecbnology Foundation

Submission of Papers: Three copies of papers, maximum 6 single spaced A4 pages should be submitted to the Congress Coordination prior to December I, 1996 for reviews. Notification of acceptance will be mailed by January 15, 1997. Author with papers accepted for presentation are requested to submit a manuscript on mats by February 28, 1997 for inclusion in the Congress Proceedings, which will be available at the Congress.

Coordinating Scientist for 8oth Conferentes in Brazil and South Amer ican Countries

Prof. Álvaro T . Prata Dept. of Mecha nicai Engineering

Federal University of Santa Catar ina Florianópolis, SC 88040-900, Brazil

Tel.: (+55) (48) 234-~166 Fax: (+SS) (48) 134-1519

e-mail: prata@nrva.ufsc.br

Journal of the ~raziUan Society of echanical Sciences

SUBMISS10N

FORMAT

ILLUSTRATlONS AND TABLES

REFERENCES

Information for Authors

• The purpose of the Journal ol the Brazlllan Society of Mecha nicai Sclences Is to publish papers of permanent lnterest <leallng wlth research. development and desion related to sclence and technology ln Mechanlcal Engineering, encompasslng interf~ees wltlt Civil, Electrlcal. Chemlcal. Naval. Nuclear. Aoricullural, Materiais. Petroleum. Aerospace, Food. System Englneenng. etc .• as well as wlth Physlcs and Applled Mathematk:S. • The JoumaJ publlshes FuHength Plj)efS, Revlew Papers and Letters to lhe Editor. Authors must agree not publish elstwhere a paper submitted to and accepted by lhe Joumal. Exceptioo can be ma<le for papers previousJy pubOshed ln proceedongs of conferences. ln thls case " should be cKed as a tootnote on lhe mte page. Copie$ oC the conterence reterees reviews should be liso locluded Revlew artlcles should constltute a criticaJ appralsal of publlshe<l rnfomlllion. • lhe declslon o! acceptance !()( publicallon lles with lhe Eelí!ors and Is baStd on lhe recommendatlons of at least two ad hoc revlewtrs. and of lhe EdkorlaJ Board. n necessary. • Manuscrlpts and ali the correspondence shout<l be sent to the Edhor or, allernabvely, to 11\e approprlate Associate E<litor. • four (4) copies of the manuscr1pt 3/t requlre<l. lhe aUthor shoutd submll tne original noures. wlllch wlll be retumed ff lhe paper is not accepted a1ter the review process. • Manuscripts should be submltted ln English or Port~guese. Spanlsh wlll also be consldered • A manuscrlpt submltted lor publlcatlon should be accompanied by a cover lener containlng lhe full name(s) of authOr(s), malling adresses, the author for contact. lncludlno phone and fax nuonber, and. n the authors so wish, lhe names of up to flva persons who could act as rel&rees. • Manuscripts should begln wHh the tille. lnctudlng the engllsh title. lhe abslract and up to tive key words. llthe pape(s language Is not Engllsh, an extended summary of abo~t 500 words should be lncluded. lhe manuscript should not comaln the aulhor(s) nal!lf(s). • ln researth papers. s~fflclenllnlormatlon should be provl<led ln lhe text or by rtletlng to papers ln generally avaitable Joumals to permh lhe worlt to be rtpeated. • Manuscrlpts should be typed doubte·spaced. on one slc!e of the page, uslng A·4 sizl!d paper, with 2 cm margtns. The pages shollld be numbertd and not to exti!S$ 24 pages, lnclu<llng llbles and fiQ111es. lhe tead aUthor of a RBCM paper whtch exceeds 11\e standat<l rtngm ot pages wll be assessed a ex~s page cha!Qe. • All symbols s/loul<l be dellne<l ln the text A separate nomendature sectlon s/lould 1151. ín alphabetlcal arder. the symbots used ln the I!X! an<l lhelr definitions. lhe greek symbots follow lhe Eng llsh symbols. and are foltowed by the subscrlpts and superscrlpts. Each dimensional symbol mu$1 hive St (Metrlc) unlls mentloned at lhe eoo. ln addltton. Enollsh unlls may be lncl~de<l parenthetlcally. Oimenslonless groups and coefficients must be so indlcated IS dlmenslonless aner thelr deflnitlon. • Uncertalntles shOuld be speclfled lor experimental and numerlcal resuHs. • Figures and lables should be relerred ln conseculíve arablc numerais. They shoulo have a caption and be place<l IS close as posslble to the teiCI flrst relerence. • LI~ drawings should be prepared on traclng paper or vellum. uslng lndla lnk; tine work muSI be even and btack. Laser prtnt output Is acceptable. lhe drawings wlth technlcal datalresulls should have a boundaoy Õn ali four sides wlth scale lndlcators (licl< marks) on ali four sl<les. The legend for lhe data symbols should be put ln the figure as wen as labels for each tuM wnerever possible. • IUIISiratlons should not be llrger lhan 12 x 17 cm. lettering shoUI<I be !aoge enough to be cleal1y legible (1.5-2.0 mm). • Photographs must be glossy prlnts. • References sl\outd be clled ln lhe text by gMng tlle ta.st rwne of lhe aulhor(s) and lhe year ol publlcallon ai the merence· efthe< 'flecelll wonc (Smith and Jones, t985) ... ar 1\ecently Smith and Jones (1985) Wlth fout or more names. use tlle form 'Smfth et ai (1985)" ln lhe text When two or more rtferences would have lhe sarne le.xt l<lentlftcaUon, dostlnguish them by appending •a·. "b', etc .• to lhe year of publlcatíon • Acceptablt references include· Joumal artlcles.<lissertations. published conffrence proceedings, numbered paper preprlnts from conf1rences, books. submltled artlctes 11 the Journal Is ldentifled. and prlvate communlcaUons. • Referentes Should be llsled ln atphabetlcal order. accordlng to lhe last name ot lhe flrst author. at the end of pape r. Some sample referentes loltow:

Bordalo, S.N., Ferzlgeo. J.H and Kllne. S.J., 1989, '1lle Oevelopment of Zonal Models for Turbulence' , Proceeóings. 1 o• ABCM • Mechanlcal Englneering Corrlerence. Vol. I, Rio de Janeiro. Brazll, pp. 4t ·44.

Clarlt, J.A.. 1986, Prtvate Communtcatlon, Univershy of Michigan, Ann Arbor, MI. Coimbra. A. L. 1978, ·ussons of Contlnuum Me<;hanlcs', Editora Edgard Biucher Uda, Slo PaulO, Brazil. Kan<lllbr. S G ano Shah, R.K., 1989, 'Asymptotic Effer::tiveness • NTU Formulas for Multiphase Plale Hell

E<changers' . ASME Joumal ol HUI T11nster, Vot t 11, pp. 3 14·321 McGormak. R.W, 1988, 'On lhe Dewlopment of El1iclent Algorilhms for Three 01menslorW Ruld floW", JoumaJ

of The Brulltan Society Of Mel;llanlcal Scieras, Vot. 10. pp. 323·346 SilVa, LH.M .. 1988 'New Integral Formulalion for Ploblems ln Mechanlcs', (ln portov~.ae), Ph.O. Thesls,

feeleraJ Unlwrsity o1 Santa Catarina, Florianópolls, SC, Brlllt Sparrow, E.M .• 19801, 'Forced-convectton Heat Transler ln a Ouct Havtng Spanwtse-Perlodic RectangUlar

ProtubertnCes'. Numerlc.tl Heat Transffr. Vai. 3. pp. 149-167. Sparrow, E.M .. t980b, 'flultHo-Fiuld Conjugale Heat Transfer for a Vertical Plpt~nternat forced Con'ler::llon

and Ext!mal Natural Convtctton·. ASMEJoumal of Hut Transler, Vol 102, pp. 402-407.

JOURNAL OFTHE BRAZILIAN SOCIETY OF MECHANJCAL SCIENCES REVISTA BRAS/LEJRA DE CIÊNCIAS MECÂNICAS

VOL XVIII- Nt..4- OECEMBER 1996

Robust Contrai • Postlton Conlrul ui Mechantca Mampulalors

Using the Opttmat Aoot Locus

Dynamlcs Mod•llng • Rililw~yWheetsel Dy'latrucs

(ln PortUQUt:SC!

Cogtmeratlon • Ca~"'lera:ton Systcru Design Opltmtzat.Jn

Turbulence • Retem Ailvances ano Apphcahon ol T :Jrhulenc.e

Dep~h-l~legrated Two EQUillio'l Cl·:>sure lor Modelrng Conlaminanl Disc~.arges

Motlon of Flo6tlng Bodies ln Wares • Added Mass and Dafl'lptng on Bonon· Mounled

Rectangular Cyhnders

Analysls ot Modlfled Structures • Vibra!ion AnalyStS ol Slruclures Su01e:.ted

lo Boun!lilry Condtlton Modtlícatton Using Exoewnema' Dala

PockiJt Ml/llng • Use Jlltle Vorcnot Dtagram lor Oeltntng ltte

CJlltt'Q Path (ln PortuguesCJ

Turblne Deslgn • Destyn ol Turbtne l:llade RJws Ustng a

Time-Marchtng MP.lhod (ln Portuguese)

Abstrscts • V o/. 18- No. 4- December 1996

RllfBfBBS • V o/. 18 · 1996

Author JndeJ· Vol. 18 · 1996

J::JseJaune Ctunmd Carios Alberlo Maluoka

Roherto Spinora BarboSil and •\I varo Costa Neto

José Anlonto Per relia &lesllen anó Paulo de Barros Correta

Lirer Yu an:l Anlonto Marozzi Rtytlelto

Sérgio Hamilton Sphater anc P"dulo de Tarso T Esperança

Oorrtngos Al>es Rade Jl'd Gerar:! L arrement

Marcos S. G hu111K1 ;:me Lucas Antonto Moscato

Jo~c; Eouardo de Barros Teixetra Bor·)es. LUIS Ma'luei de CarvalhO Gato and HUt Manue! R:m1a ôe Jesus Pereu~

308

318

330

338

355

374

383

395

403

405

4D/