let. - vol. 48 (2002) {t. - no. 2 STROJNIŠKI VESTNIK JOURNAL OF MECHANICAL ENGINEERING strani - pages 57 - 118 ISSN 0039-2480 . Stroj V . STJVAX cena 800 SIT 00 ¦* CO o Robno obmo~na integralska metoda a numeri~no modeliranje lebde~ih s--lo--jev Boundary Domain Integral Method for Numerical Modeling of Fluidized Beds Model za analizo in optimiranje vpenjalnih priprav ----A Model for Analysing and Optimazing Fixtures Meritev gibanja kolena z industrijskim robotom - avtomatska kompenzacija gravitacije prijemala ------- Measuring knee movement using an industrial robot gravity compensation Realni preto~ni ~asi operacij in NKP sistema Lead Times of Operations and Efficiency of the PPC System 1. 2. 3. 4. © Strojni{ki vestnik 48(2002)2,57 Mese~nik ISSN 0039-2480 © Journal of Mechanical Engineering 48(2002)2,57 Published monthly ISSN 0039-2480 Vsebina Contents Strojni{ki "vestnik" - Journal of Mechanical Engineering letnik - volume 48, (2002), {tevilka - number 2 Razprave Požarnik, M., Škerget, L.: Robno območna integralska metoda za numerično modeliranje lebdečih slojev 58 Župerl, U., Čuš, F.: Model za analizo in optimiranje vpenjalnih priprav 73 Omrčen, D., Nemec, B.: Meritev gibanja kolena z industrijskim robotom - avtomatska kompenzacija gravitacije prijemala 87 Starbek, M., Grum, J., Kušar, J.: Realni pretočni časi operacij in uspešnost sistema NKP 98 Strokovna literatura Osebne vesti Navodila avtorjem 112 116 117 Papers Požarnik, M., Škerget, L.: Boundary Domain Integral Method for Numerical Modeling of Fluidized Beds Župerl, U., Čuš, F.: A Model for Analysing and Optimazing Fixtures Omrčen, D., Nemec, B.: Measuring knee movement using an industrial robot - gravity compensation for the automatic gripper Starbek M., Grum, J., Kušar, J.: Realistic Lead Times of Operations and Efficiency of the PPC System Professional Literature Personal Events Instructions for Authors grinpaBiJjiMMKni stran 57 glTMDDC © Strojni{ki vestnik 48(2002)2,58-72 © Journal of Mechanical Engineering 48(2002)2,58-72 ISSN 0039-2480 ISSN 0039-2480 UDK 532.5:519.61/.64 UDC 532.5:519.61/.64 Izvirni znanstveni ~lanek (1.01) Original scientific paper (1.01) Robno obmo~na integralska metoda za numeri~no modeliranje lebde~ih slojev Boundary Domain Integral Method for Numerical Modeling of Fluidized Beds Matej Po`arnik - Leopold [kerget V prispevku je prikazano numerično reševanje dinamike dvofaznih dvosestavinskih tokov z robno območno integralsko metodo. Model opisa gibanja sestavin temelji na modelu dveh tekočin s hitrostno-vrtinčno formulacijo dopolnjenih Navier-Stokesovih enačb. Poseben poudarek je namenjen členu medfazne izmenjave gibalne količine. Kot testna primera sta prikazana enofazni tok v kanalu z nenadno simetrično razširitvijo in dvofazni dvosestavinski tok v navpičnem kanalu. © 2002 Strojniški vestnik. Vse pravice pridržane. (Ključne besede: tok dvofazni, modeliranje numerično, metode robno-območne, metode integralske, enačbe Navier-Stokes) The paper deals with the numerical modeling of two-phase two-component flows using the boundary domain integral method. The two-fluid model with the velocity-vorticity formulation of modified Navier-Stokes equations is adopted. Particular attention is given to the interphase momentum exchange term. As test cases a single-phase symmetric sudden expansion flow and two-phase two-component vertical channel flow are investigated. © 2002 Journal of Mechanical Engineering. All rights reserved. (Keywords: two phase flow, numerical modeling, boundary domain integral methods, Navier-Stokes equations) 0 UVOD Lebdeči sloj sestavljata navzgor gibajoča se tekočina, ponavadi plin in na nosilno ploščo nasut sloj polnil. Kljub temu, da se trdni delci, ko hitrost plina preseže najmanjšo hitrost lebdenja, večino časa še vedno dotikajo, se mešanica plina in trdnih delcev obnaša kot tekočina. Tlak v zmesi se povečuje linearno z razdaljo pod površino, težji delci tonejo, lažji se dvigujejo, opaziti je mogoče gibanje v obliki valov. Trdne delce, ki jih imenujemo polnila, lahko stalno dodajamo ali odvzemamo. Vsa drobna polnila imajo izredno veliko specifično površino; 1 m3 delcev premera 104 m ima površino 30 000 m2. Zelo pomembno je nemirno delovanje plinskih mehurčkov ki skrbi za popolno mešanje polnil. Posledica sta veliki toplotni in snovski prestopnosti med površino in lebdečim slojem ter med plinom in polnili, ki povzročita, da je temperaturno in koncentracijsko polje homogeno tako v prečni kakor tudi vzdolžni smeri. Če primerjamo lebdeči sloj z nasutim slojem enakih polnil pri enaki višini sloja in hitrosti plina ugotovimo, da je tlačni padec v lebdečem sloju veliko manjši. Zaradi vsega naštetega so lebdeči sloji oziroma dvofazni dvosestavinski tokovi plin-trdni delci zelo privlačno procesno orodje. 0 INTRODUCTION A fluidized bed is formed by passing a fluid, usually a gas, upwards through a bed of particles that are supported on a distributor. Even though above the minimum fluidization velocity the particles are touching each other most of the time the interparticle friction is so small that the fluid/solid assembly behaves like a fluid. The pressure increases linearly with the distance below the surface: denser objects sink, lighter ones float, and a wave motion is observed. Solids can be removed from or added to the bed continuously, and this provides many processing advantages. All fine powders have a very large specific surface area — 1 m3 of 10-4 m particles has a surface area of about 30000 m2 — but in a fluidized bed the stirring action of the gas bubbles continuously moves the powder around, shearing it and exposing it to the gas. This excellent solids mixing gives the high rates of heat transfer from the surface to the bed and from the gas to particles resulting in isothermal conditions: both radially and axially. Compared with a fixed bed of the same powder operated at the same bed depth and gas velocity, the pressure drop over a fluidized bed is much smaller, and this together with most of the other characteristics make the fluidized bed an attractive choice as a chemical or physical processing tool. grin^(afcflM]SCLD VH^tTPsDDIK stran 58 Po`arnik M. - [kerget L.: Robno obmo~na integralska metoda - Boundary Domain Integral Method Računalniška dinamika tekočin (RDT) postaja Computational fluid dynamics (CFD) is vse pomembnejše orodje za določevanje tokov v becoming an increasingly common engineering tool različnih vrstah industrijskih naprav. Kljub temu, da to predict flows in various types of apparatus on an so programski paketi za modeliranje enofaznih tokov industrial scale. Although the tools for applying široko dosegljivi, modeliranje večfaznih sistemov še single-phase CFD are widely available, the application vedno pomeni veliko težavo, tako s fizikalnega kakor of multiphase CFD is still complicated from both the z numeričnega vidika. Avtorjema znani numerični physical and the numerical points of view. All the algoritmi za modeliranje dvofaznih dvosestavinskih numerical algorithms known to the authors that were tokov brez izjeme temeljijo na postopkih končnih razlik, developed so far to simulate two-phase two- končnih elementov oziroma kontrolnih prostornin. component flows are strictly based on the finite- V prispevku podajamo razvoj alternativne difference method, the finite-element method or the numerične sheme na podlagi robno območne control-volume method. integralske metode (ROIM) za reševanje splošnega In this paper an alternative numerical scheme primera gibanja dvofaznega dvosestavinskega based on the boundary domain integral method sistema. Prispevek predstavlja prvi primer uporabe (BDIM) is presented for the solution of a general tako ROIM kakor tudi hitrostno-vrtinčne formulacije two-phase two-component flow motion problem. This za modeliranje dvofaznih tokov. Kot podlago za is definitely the first attempt to implement the BDIM izpeljavo dopolnjenih sistemov Navier-Stokesovih and velocity-vorticity formulation for modeling two- enačb smo uporabili v literaturi dobro znan model phase flows. The two-fluid model (TFM) is used to dveh tekočin (MDT). Prednost hitrostno-vrtinčne derive two sets of modified Navier-Stokes equations. formulacije omenjenega matematičnega modela The velocity-vorticity formulation of the physical fizikalnih zakonov ohranitve mase in gibalne količine conservation laws of mass and momentum then je, ob določenih dodatnih predpostavkah, numerična follows. The advantages of this approach lie with the ločitev kinematike in kinetike toka obeh sestavin od numerical separation of kinematic and kinetic aspects računanja termodinamičnega tlaka in navideznega of both phases motion from the thermodynamic tlaka trdne snovi. Pomembna za oblikovanje hitrostnih pressure and the solid’s pressure computation. in vrtinčnih polj sestavin je izmenjava gibalne količine Particular attention is given to the drag between the med sestavinama, ki jo opišemo s koeficientom phases, which is described by the interphase medfazne izmenjave gibalne količine. momentum exchange coefficient. 1 DVOFAZNI DVOSESTAVINSKI MODEL 1 TWO-PHASE TWO-COMPONENT MODEL Število delcev v lebdečem sloju realne In spite of increasing computational power velikosti je kljub naraščajočim računalniškim the number of particles in a gas-solid flow in large zmogljivostim preveliko, da bi z Lagrangeovo metodo scale equipment is still much too large to handle each modelirali gibanje vsakega delca posebej. Tak particle separately. Simulating each particle separately postopek omogoči študij mikroskopskih lastnosti is called a Lagrangian method, which can be used to lebdečega sloja. V prispevku prikazana shema temelji study the microscopic properties of fluidized beds. na modelu dveh tekočin, dopolnjenim s teorijo The CFD model used in this work is based on a TFM gnanega toka dvofaznih tokov. V MDT obe sestavini that is extended with the drift-flux theory of a two- obravnavamo kot zvezni in medsebojno popolnoma phase flow. In a TFM both phases are considered to pronicajoči. Model dveh tekočin sta prva predstavila be continuous and fully interpenetrating. The TFM Anderson in Jackson [1]. was first proposed by Anderson and Jackson [1]. 1.1 Zapis za osnovne spremenljivke 1.1 Primitive variables formulation Kontinuitetno enačbo za sestavino p (f za plin The continuity equation or mass balance for in s za trdno snov) zapišemo kot: phase p (f for gas and s for solid) reads: Izmenjava snovi med sestavinama, npr. zaradi Mass exchange between the phases, e.g. due kemijske reakcije ali zgorevanja, ni zajeta. to reaction or combustion, is not considered. Gibalna enačba plinaste sestavine je podana The momentum balance for the gas phase is kot dopolnjena Navier-Stokesova enačba, ki vsebuje given by the Navier-Stokes equation, which is modified člen medfazne izmenjave gibalne količine: to include an interphase momentum transfer term: stran 59 0^BBfirTMK Po`arnik M. - [kerget L.: Robno obmo~na integralska metoda - Boundary Domain Integral Method Gibalna enačba trdne snovi je podana z: The solid-phase momentum balance is given by: kjer je ps* navidezni tlak trdne snovi, ki je določen z uporabo kinetične teorije zrnatega toka. V tenzorju strižnih napetosti moramo v splošnem upoštevati tako strižno kakor tudi normalno viskoznost [2]. Z ef = 1 in fi= 0 postane enačba (2) običajna Navier-Stokesova enačba. Kontinuitetni in gibalni enačbi sta podrobno predstavljeni v [3]. 1.2 Hitrostno-vrtinčna formulacija Za uporabo v predlagani shemi na podlagi ROIM originalna sistema dopolnjenih Navier-Stokesovih enačb nadalje preoblikujemo z uporabo hitrostno-vrtinčne formulacije. Na tak način računsko shemo gibanja plina in trdne snovi razdelimo na kinematični in kinetični vidik. V primeru uporabe hitrostno-vrtinčne formulacije za modeliranje enofaznih tokov iz numerične sheme izločimo računanje tlaka. Postopek vodi v primerjavi z običajnim k preprostejšemu predpisovanju robnih pogojev predvsem na tistih robovih območja, kjer tlak ni znan. Izpeljani algoritem na podlagi hitrostno-vrtinčne formulacije je kljub vsemu še vedno mogoče zapisati popolnoma splošno za modeliranje dvo-in tridimenzionalnih tokov. Računsko shemo gibanja obeh sestavin razdelimo na kinematični in kinetični vidik z vpeljavo vektorjev vrtinčnosti mpi, ki predstavljata rotorja hitrostnih polj. Z uporabo operatorja rotor neposredno na definiciji vrtinčnosti in upoštevanjem preoblikovane kontinuitetne enačbe (1) z ep = konst.: (3), where ps* is the solids pressure, originally obtained from the kinetic theory of the granular flow. Both the shear and bulk viscosities should be used in a viscous strain-rate tensor, in general [2]. With ef = 1 and b = 0 Eqn. (2) becomes the classical Navier-Stokes equation. The mass and momentum balances are discussed in detail in [3]. 1.2 Velocity-vorticity variables formulation In BDIM the original sets of Navier-Stokes equations for the gas phase and the solid particles are further transformed with the use of the velocity-vorticity variables formulation. Within this approach the flow field computation is decoupled into flow kinematics and flow kinetics. The main advantages of this scheme in the case of single-phase flow lie with the numerical separation of the kinematic and kinetic aspects of the flow from the pressure computation. This leads to a simpler way of enforcing the proper boundary conditions than the primitive variables approach whenever the pressure is not specified on the boundary as a known quantity. The developed algorithm can still be written in a general form for both two and three dimensions. With the vorticity vector wpi, representing the curl of the velocity field, the two phases motion computation scheme is partitioned into its kinematic and kinetic aspects. By taking the curl operator directly to the vorticity vector definition, and applying the reformed continuity Eqn. (1) with ep = const.: 'hi lir- M = _10^+1^M (4) izpeljemo kinematiko gibanja plinaste in trdne snovi: the kinematics of both phases motion is carried out: (5). Eqn. (5) represents the kinematics of an incompressible fluid and solid-phase motion or the compatibility of the velocity and vorticity fields at a given point in space and time. The kinetics are governed by the vorticity transport equations obtained as a curl of the momentum balances, Eqns. (2) and (3). In the case of low solid concentrations the approach of Chapman and Cowling [4] with constant viscosities is applied. The vorticity transport equations can be written in the following form: Enačba (5) podaja kinematiko gibanja nestisljive tekočine in trdnih delcev oziroma združljivost hitrostnih in vrtinčnih polj v dani točki prostora in časa. Kinetiko podamo s prenosnima enačbama vrtinčnosti, ki pomenita rotor gibalnih enačb (2) in (3). V primeru majhnih prostorninskih deležev trdnih delcev v toku lebdečega sloja uporabimo postopek Chapmana in Cowlinga [4] z nespremenljivima strižno in normalno viskoznostjo. Prenosni enačbi vrtinčnosti zapišemo kot: (6), grin^(afcflM]SCLD VBgfFMK stran 60 Po`arnik M. - [kerget L.: Robno obmo~na integralska metoda - Boundary Domain Integral Method ki podajata porazdelitev vektorjev vrtinčnosti plinaste sestavine in trdne snovi. Z opisanim postopkom iz izračuna izločimo termodinamični tlak p, vendar enačba (7) še vedno vsebuje gradient navideznega tlaka trdne snovi p . Zaradi tega predlagana shema temelji na tehniki podobmočij v njeni skrajni izvedenki, kar pomeni, da je vsaka celica podobmočje, imenovano makroelement, ki je omejeno s štirimi robnimi elementi. S predpostavljenim nespremenljivim prostorninskim deležem sestavine ep znotraj posamezne iteracije po vsakem izmed makroelementov zagotovimo, da gradient prostorninskega deleža po območju ne obstaja (dep/dxj = 0). Z upoštevanjem primernih vmesnih pogojev na mejah med makroelementi, ki so podrobneje predstavljeni v poglavlju 1.4, enačbe (5), (6) in (7) zapišemo v poenostavljeni obliki. Po predpostavki glede prostorninskega deleža ep kinematiko gibanja obeh sestavin zapišemo v obliki parabolične enačbe z uporabljenim nepravim neustaljenim načinom: Transportni enačbi vrtinčnosti, enačbi (6) in (7), prepišemo v naslednjo obliko: (7), describing the redistribution of the vorticity vector in the fluid and solid particles’ flow field. While the thermodynamic pressure p is not part of the computation, Eqn. (7) is still dealing with the solids pressure p gradient. Therefore, the proposed numerical scheme is based on the subdomain technique in its limit version. Each internal cell represents one subdomain called a macroelement, assumed to be constant within a particular iteration of the numerical algorithm, therefore, in that moment the gradient dep/dxj = 0 does not exist all over the macroelement. Then the Eqns. (5), (6), and (7) can be rewritten in a simple manner, but a lot of the physics is moved to the macroelement interface boundary conditions explained in detail in section 1.4. After the assumption regarding the volume fraction ep the kinematics of both phases motion is written in the sense of the parabolic equation where the false transient approach is implemented afterwards: (8). Using the same principle the kinetics given by Eqns. (6) and (7) is rewritten as: Dr (~y< -i*J«J / U d-UlS: WsiSe, dum P , » Dr "-8xjdxj+^B&F^si8xj+QBTB *"* ~ "*] =[/,_ (9), (10). Enačbi (9) in (10) podajata v primeru tridimenzionalnih tokov časovno spremembo vrtinčnosti delca plina oziroma trdne snovi, podano s Stokesovim odvodom na levi strani enačb, zaradi učinkov viskozne difuzije, nastajanja mehurjev, učinkov deformacije in medfazne izmenjave gibalne količine, podane s členi na desni strani enačb. 1.3 Prostorninski delež Za sklenitev sistema kinematičnih in kinetičnih enačb gibanja sestavin dvofaznega dvosestavinskega toka potrebujemo dodatno enačbo za izračun prostorninskega deleža plinaste sestavine. Izpeljemo jo iz teorije gnanega toka dvofaznih tokov [9]. Teorija je izpeljana splošno in obsega različnost fizikalnih lastnosti in tokovnih polj dvofaznega toka. Omogoča zapis ločenih kontinuitetnih in gibalnih enačb za vsako sestavino posebej. Uporabna je za modeliranje različnih režimov dvofaznih tokov plin -kapljevina kakor tudi za modeliranje dvofaznih dvosestavinskih sistemov tekočina - trdni delci, npr. Eqns. (9) and (10) show that the rate of change of the vorticity in the case of 3D flow as one follows a fluid or solid particle, given by the Stokes derivation on the left-hand side of the equations, is due to the viscous diffusion, bubble formation, vortex twisting and stretching, and interphase-momentum transfer, represented by the terms on the right-hand side. 1.3 Volume fraction To close the system of kinematics and kinetics equations an additional equation to compute the volume fraction of the fluid phase is derived from the drift flux theory of the two-phase flow [9]. This model treats the general case of modelling each phase or component as a separate fluid with its own set of governing balance equations. In general, each phase has its own velocity, vorticity, and temperature. Drift-flux theory has widespread application in the bubbly, slug, and drop regimes of gas-liquid flow as well as to fluid-particle systems such as fluidized beds. It provides a starting point for the extension of the grin^otJJiMiscsni stran 61 glTMDDC Po`arnik M. - [kerget L.: Robno obmo~na integralska metoda - Boundary Domain Integral Method lebdeči sloji. Namenjena je tudi kot izhodišče za preučevanje tokov s prevladujočimi 2D in 3D vplivi. Prostorninski delež sestavine določimo z enačbo: e/ = " kjer je v Stokesova hitrost posameznega trdnega delca v mirujoči tekočini. Določitev vrednosti koeficienta n podajata Richardson in Zaki [7]. Celotna korelacija Richardsona in Zakija za vse vrednosti Reynoldsovega števila podaja vrednosti n med 4,65 in 2,39, upoštevajoč, da so trdni delci toge kroglice majhnega premera v primerjavi z izmerami kanala. Vrednost koeficienta n lahko še povečamo, če prihaja do združevanja delcev. V predstavljeni raziskavi uporabljamo srednjo vrednost koeficienta n za sisteme plin - trdni delci (n=3). Vpeljati je mogoče tudi korekcijski faktor v odvisnosti od razmerja med premerom delcev in izmerami kanala. 1.4 Pogoji vmesnega roba Kadar zapišemo dopolnjen sistem Navier-Stokesovih enačb gibanja dvofaznega dvosestavinskega sistema za nespremenljiv prostorninski delež sestavine v podobmočju v iteraciji numeričnega algoritma, so najbolj kritičen del izračuna pogoji vmesnega roba na mejah med makroelementi. Zaradi tega smo posebno pozornost namenili analizi nezveznosti. Vektorje hitrosti razdelimo na normalno in obodno komponento glede na rob podobmočja. Zaradi predpisanega nespremenljivega prostorninskega deleža znotaj posamezne iteracije je glavna značilnost vmesnega roba H med podobmočjema A in Q2 skok prostorninskega deleža sestavine in posledično tudi skok normalne komponente hitrosti. Obodne komponente vektorjev hitrosti se spreminjajo zvezno. Za gibanje v ravnini x-y velja naslednja povezava med normalnimi in tangentnimi odvodi normalne in obodne komponente vektorja hitrosti: theory to flows in which two- and three-dimensional effects are significant. The volume fraction of the component is determined by following equation: , \V/i-1>ri\ •: (11), Iz ohranitvenih zakonov izpeljemo ustrezne pogoje vmesnega roba kinematike: where v is the terminal speed of a single solid particle in an infinite stationary liquid. The evaluation of index n was shown by Richardson and Zaki [7]. The complete correlation of Richardson and Zaki over the whole range of Reynolds numbers gives a value for the index n between 4,65 and 2,39, assuming that the particles are rigid spheres and small compared to the diameter of the channel. The value of n can also be enlarged when the particles flocculate. In our study only small Reynolds number flow is encountered, therefore, an intermediate value of the index n for fluid-particle systems is taken into account (n=3). A correction factor can also be introduced in terms of the ratio of the particle diameter to the tube diameter. 1.4 Macroelement interface boundary conditions When dealing with the modified Navier-Stokes system of equations written for the constant volume fraction of the component over each macroelement within a particular iteration the most critical parts of the numerical scheme are the macroelement interface boundary conditions. Therefore, particular attention has been given to the analysis of a jump discontinuity in the flow properties. Velocity vectors at the macroelement interface boundaries are split into normal and tangential components. The major interface characteristics are a jump in the volume fraction and the continuity of the tangential velocity (no-slip). Due to the volume fraction jump the jump in normal velocity is also necessary. For the macroelement interface the relation relating the vorticity values with the normal and tangential velocity component fluxes is derived. For two-dimensional motion in the x-y plane it reads: (12). Appropriate macroelement interface boundary conditions derived from the conservation laws for the kinematics are: in kinetike * On and for the kinetics are: (13), (14), (15), (16) (17), grin^(afcflM]SCLD VBgfFMK stran 62 Po`arnik M. - [kerget L.: Robno obmo~na integralska metoda - Boundary Domain Integral Method lL' t\_ i'L- Z I- 8n On S predpisanim ep = 1 na obeh straneh vmesnega roba /I med podobmočjema A in Q2 enačbe od (13) do (18) privzamejo običajno obliko pogojev vmesnega roba za enofazno gibanje viskozne nestisljive tekočine. 1.5 Robni pogoji Za numerično rešitev sistema enačb gibanja dvofaznega dvosestavinskega toka plin - trdna snov je treba predpisati primerne robne pogoje za obe sestavini. Običajni robni pogoji v primeru zapisa za osnovne spremenljivke so znane vrednosti robnih hitrosti obeh sestavin, termodinamični tlak plinaste sestavine in navidezna temperatura trdne snovi, ki je izpeljana iz kinetične teorije zrnatega toka in je namenjena za določitev navideznega tlaka trdne snovi. Pri zapisu za hitrostno-vrtinčno formulacijo potrebujemo za obe sestavini le znane vrednosti robnih hitrosti oziroma vrednosti njihovih odvodov. Robni pogoji običajne ROIM so podrobno predstavljeni v [8]. Na neporozni trdni steni kanala pri toku plina in trdnih delcev predpišemo po navadi brezzdrsne robne pogoje plinaste sestavine. Tega ne moremo vedno storiti tudi za trdne delce, saj pri toku delcev večjega premera v primerjavi s hrapavostjo sten prihaja do zdrsa pri trku ob steno. Mnogokrat kljub temu predpostavimo, da so delci izjemno majhni, kar omogoči predpisovanje brezzdrsnih robnih pogojev tudi za trdno snov. Začetni pogoji prostorninskega deleža sestavin lahko zavzamejo katerokoli fizikalno sprejemljivo vrednost. 2 INTEGRALSKA PREDSTAVITEV HITROSTNO-VRTINCNE FORMULACIJE DVOFAZNEGA DVOSESTAVINSKEGA TOKA 2.1 Kinematika gibanja plina in trdne snovi (18). where' subscript I denotes the interface boundary H between the macroelements a and /22. With ep = 1 at both sides of the boundary 7I the Eqns. from (13) to (18) take on the classical form for single-phase fluid motion. 1.5 Boundary conditions To solve the equations for gas-solid particles flow we need appropriate boundary conditions, not only for the fluid phase but also for the solid particles. Classical boundary conditions in the case of primitive variables formulation are the velocities of the gas phase and the solid particles, the boundary condition for gas-phase pressure and the boundary granular temperature derived from the kinetic theory of granular flow to compute the solids pressure. In the velocity-vorticity variables approach we need only the appropriate boundary conditions for the velocities or velocity fluxes of the gas and the solids. The boundary conditions for the classical BDIM approach are discussed in detail in [8]. For the gas-solid particle flow motion the gas-phase velocities are generally set to zero at an impenetrable rigid wall. This no-slip condition cannot always be applied to solid motion. Since the particle diameter is usually larger than the length scale of the surface roughness of the rigid wall, the particles may partially slip the wall. But it is also important to note that for small particle diameters the boundary condition is close to the no-slip condition. Initially, in a mixture of gas phase and solid particles the volume fraction of the two phases can be set to any physically acceptable value. 2 INTEGRAL REPRESENTATION OF THE VELOCITY-VORTICITY FORMULATION FOR TWO-PHASE TWO-COMPONENT FLOW 2.1 Gas and solid phase kinematics Pri obravnavi kinematike obeh sestavin v Considering first the kinematics of both robno-območni integralski predstavitvi phases motion in the boundary domain integral upoštevamo, da vsaka komponenta vektorjev representation one has to take into account that each hitrosti v enačbi (8) zadošča nehomogeni component of the velocity vector in Eqn. (8) obeys parabolični enačbi: the nonhomogeneous parabolic equation: ^e^&r^+^0 (19). Z enačenjem vektorja navidezne telesne sile bpi z vrtinčnim delom toka obeh sestavin: zapišemo pripadajoč integralski izraz: Equating the body force term bpi to the vortical fluid or solid particles flow term: (20) renders an integral statement : gfin^OtJJlMISCSD stran 63 Po`arnik M. - [kerget L.: Robno obmo~na integralska metoda - Boundary Domain Integral Method kjer je u* parabolična difuzijska osnovna rešitev. Če predpostavimo nespremenljiv potek vseh funkcij polja po posamičnem časovnem koraku, lahko časovne integrale v enačbi (21) rešimo analitično [8], tako da integralski stavek podamo v obliki: (21), where u*p is the parabolic diffusion fundamental solution. Assuming a constant variation of all field variables within the individual time increment, the time integrals in Eqn. (21) may be evaluated analytically [8]. An integral statement can be finally written in the following form: ov ¦+ j etjtiDfj-^-dn ¦+ j VfrF__u^F_ (22), i lil kjer je U*p = apup. 2.2 Kinetika gibanja plina in trdne snovi Integralsko predstavitev kinetike gibanja plina in trdne snovi zapišemo z uporabo difuzivno-konvektivne osnovne rešitve. Časovne odvode vrtinčnosti in prostorninskih deležev aproksimiramo z nesimetričnimi končnimi razlikami. Tak način omogoči zapis enačb (9) in (10) v obliki nehomogene difuzivno-konvektivne parcialne diferencialne enačbe: Diferencialni zapis preoblikujemo v ustrezen integralski stavek: being U*p =apu*p . 2.2 Gas and solid phase kinetics An integral representation describing the kinetics of both components motion is formulated by using the fundamental solution of a steady diffusion-convective partial differential equation with a reaction term. In the case of kinetics the volume fraction and the vorticity time derivatives are approximated by using a non-symmetric finite-difference approximation. Therefore, Eqns. (9) and (10) can be rewritten as a non-homogeneous diffusion-convective partial differential equation: (23). The above differential formulation can be transformed into an equivalent integral statement: c (O «V (O + j Upi-^-dT = — j [ Vp-rsr - <*W J U^dT + — j bpib^dtl (24) z U*=np0-up in v=vpj-nj. u* je difuzivno-konvektivna osnovna rešitev. Izraz navideznih telesnih sil bpi vsebuje konvekcijo zaradi spreminjajočega se dela vektorja hitrosti vj, deformacijo, začetne pogoje, spremembo vrtinčnosti zaradi oblikovanja mehurjev in medfazne izmenjave gibalne količine: with U*=np0-up and vpn=vpj-nj. up is the fundamental solution of the diffusion-convective equation. The pseudo-body term bpi includes the convective flux for the perturbated velocity field vj, deformation, initial conditions, vorticity change on account of the bubble formation and interphase momentum exchange term, for example: Hi' dx-i At + 7T:(EP:F-EP:fJ -F P{Ufi-Usi) epHr (25). Končni integralski stavek se glasi: The final integral representation is as follows: (26). grin^(afcflM]SCLD VBgfFMK stran 64 Po`arnik M. - [kerget L.: Robno obmo~na integralska metoda - Boundary Domain Integral Method 3 NUMERIČNA REŠITEV 3.1 Diskretizirane enačbe Za numerično približno rešitev obravnavanih funkcij polja, npr. hitrosti in vrtinčnosti obeh sestavin, moramo pripadajoče robno-območne integralske predstavitve zapisati v diskretni obliki. Diskretizacijo integralskih predstavitev hitrosti in vrtinčnosti obeh sestavin izpeljemo iz pripadajočih integralskih enačb (22) in (26), za kinematiko: 3 NUMERICAL SOLUTION 3.1 Discretized equations Consider a discretized equation set for the case of both fluid and solid particles motion. The discretization of the integral representations for velocities and vorticities can be readily obtained from the corresponding integral Eqns. (22) and (26) as follows, for the kinematics: in kinetiko C - ey* E {diri}T i <*>*}" ¦+ C i EIW1^.^—}" and for the vorticity kinetics: (27) c (0o* (0 + E {*«TH>}" = ^ E {C„}T 4 " Src " uJjjif jm 4 V1 i 1 C u„_Hj (28), v„_!\t C T — E{bpJ {">/¦-">«}" kjer nadpis T označuje transpozicijo. ^using superscript T to denote the transposition. Kinematiko ravninskega toka (i,j=1,2) podamo The plane kinematics (i,j=1,2) can be given z naslednjim stavkom po enačbi (27): by the following statement based on Eqn. (27): E E f g„ c (0 V K, Tj?) + E (V}T {^j}" = E {OeJT 1 "gf + eH^vni C EM H}" + (29). Diskretizirane transportne enačbe vrtinčnosti obeh sestavin ravninskega toka so: The plane vorticity kinetics reads: (30). grin^otJJiMiscsni stran 65 glTMDDC Po`arnik M. - [kerget L.: Robno obmo~na integralska metoda - Boundary Domain Integral Method Postopek izračunavanja integralov in zapis diskretiziranih integralskih enačb je običajen in je podrobno opisan v [6] in [8]. 3.2 Postopek reševanja vezanih sistemov enačb Zapisana sistema enačb kinematike, enačba (29), ter kinetike plinaste in trdne sestavine, enačba (30), dvofaznega dvosestavinskega lebdečega sloja, povezana z enačbo za izračun prostorninskega deleža sestavine (11), sestavljata močno nelinearen sistem enačb, katerih rešitev poiščemo na iterativen način. Postopek je naslednji: (F - plinasta sestavina, S - trdna sestavina) 1.FPričnemo z začetnimi predpostavljenimi vrednostmi vrtinčnosti plinaste sestavine. 1.S Pričnemo z začetnimi predpostavljenimi vrednostmi vrtinčnosti trdne snovi. 2.F Kinematika plinaste sestavine: (a) rešitev implicitnega sistema za robne vrednosti hitrosti oziroma vrednosti normalnih odvodov hitrosti plinaste sestavine, (b) transformacija novih funkcijskih vrednosti plinaste sestavine iz vozlišč elementov v vozlišča celic, (c) določitev novih vrednosti vrtinčnosti plinaste sestavine na robu, (d) določitev novih matrik kinetike plinaste sestavine, če se nespremenljivi del hitrostnega vektorja plinaste sestavine spremeni za več od predpisane tolerance. 2.S Kinematika trdne snovi: (a) rešitev implicitnega sistema za robne vrednosti hitrosti oziroma vrednosti normalnih odvodov hitrosti trdne snovi, (b) transformacija novih funkcijskih vrednosti trdne snovi iz vozlišč elementov v vozlišča celic, (c) določitev novih vrednosti vrtinčnosti trdne snovi na robu, (d) določitev novih matrik kinetike trdne snovi, če se nespremenljivi del hitrostnega vektorja trdne snovi spremeni za več od predpisane tolerance. 3.F+S Določitev prostorninskih deležev plinaste in trdne sestavine. 4.FKinetika plinaste sestavine: (a) rešitev implicitnega sistema za neznane vrednosti normalnih odvodov vrtinčnosti plinaste sestavine in notranjih območnih vrednosti vrtinčnosti plinaste sestavine, (b) transformacija novih funkcijskih vrednosti plinaste sestavine iz vozlišč elementov v vozlišča celic. 4.SKinetika trdne snovi: (a) rešitev implicitnega sistema za neznane vrednosti normalnih odvodov vrtinčnosti trdne snovi in notranjih območnih vrednosti vrtinčnosti trdne snovi, (b) transformacija novih funkcijskih vrednosti trdne snovi iz vozlišč elementov v vozlišča celic. The procedure is standard and can be seen in detail in [6] and [8]. 3.2 Solution procedure The kinematics relations, Eqn. (29), and the vorticity transport equations, Eqn. (30), for both phases motion are coupled in two sets of non-linear equations. These two sets are related to each other by Eqn. (11), which is used to compute the volume fraction of the fluid phase knowing the velocity fields of both components. In order to obtain a solution of the fluid and solid particles motion problem a sequential computational algorithm was developed. The main steps in this algorithm are: (F - gas phase, S - solid phase) 1.F Start with initial values for the fluid-phase vorticity distribution. 1.S Start with initial values for the solid-phase vorticity distribution. 2.F Fluid-phase kinematic computational part: (a) solve implicit sets for boundary fluid-phase velocity or velocity normal flux values, (b) transform new values from element nodes to cell nodes, (c) determine new boundary fluid-phase vorticity values, (d) compute new fluid-phase matrices for the kinetics if the constant fluid-phase velocity vector is perturbated more than the prescribed tolerance. 2.S Solid-phase kinematic computational part: (a) solve implicit sets for boundary solid-phase velocity or velocity normal flux values, (b) transform new values from element nodes to cell nodes, (c) determine new boundary solid-phase vorticity values, (d) compute new solid-phase matrices for the kinetics if the constant solid-phase velocity vector is perturbated more than the prescribed tolerance. 3.F+S Determine volume fraction of the fluid and solid-phase. 4.F Fluid-phase vorticity kinetic computational part: (a) solve implicit set for unknown boundary fluid-phase vorticity flux and internal domain fluid-phase vorticity values, (b) transform new values from element nodes to cell nodes. 4.S Solid-phase vorticity kinetic computational part: (a) solve implicit set for unknown boundary solid-phase vorticity flux and internal domain solid-phase vorticity values, (b) transform new values from element nodes to cell nodes. grin^(afcflM]SCLD VBgfFMK stran 66 Po`arnik M. - [kerget L.: Robno obmo~na integralska metoda - Boundary Domain Integral Method 5.FPodsprostitev vrtinčnosti plinaste sestavine in preverjanje konvergence. Če je konvergenčni kriterij izpolnjen z izračuni, prenehamo, sicer se vrnemo na korak 2.F. 5.SPodsprostitev vrtinčnosti trdne snovi in preverjanje konvergence. Če je konvergenčni kriterij izpolnjen z izračuni, prenehamo, sicer se vrnemo na korak 2.S. 4 TESTNI PRIMERI 4.1 Enofazni enosestavinski tok v kanalu z nenadno simetrično razširitvijo Tok v kanalu z nenadno simetrično razširitvijo je podoben toku v kanalu s stopnico, saj tudi tukaj pride do pojava recirkulacije pri večjih vrednostih Reynoldsovega števila. Posebnost problema je obnašanje toka pri večjih vrednostih Reynoldsovega števila, ko se prvotno simetrični recirkulacijski coni začneta razlikovati po dolžini. V našem primeru smo simulirali enofazni enosestavinski tok v kanalu z nenadno simetrično razširitvijo za vrednost Reynoldsovega števila Re=56, ki je dobro dokumentiran z eksperimentalnimi podatki (Durst [5]). Tok v kanalu z nenadno simetrično razširitvijo označuje preprosta geometrijska oblika z vstopnim območjem, omejenim na 1/3 višine kanala. Locirano je na sredini kanala, kar kaže slika 1, kjer so prikazani tudi robni pogoji in izmere kanala. 5.F Relax fluid-phase vorticity values and check the fluid-phase convergence. If convergence criterion is satisfied, then stop; otherwise, go to step 2.F. 5.S Relax solid-phase vorticity values and check the solid-phase convergence. If convergence criterion is satisfied, then stop; otherwise, go to step 2.S. 4 TEST EXAMPLES 4.1 Single-phase symmetric sudden expansion flow Because of the recirculation zones at larger values of Reynolds number flow downstream of the expansion a plane symmetric sudden expansion flow is similar to the backward-facing step flow. The flow is symmetric at sufficiently low values of the Reynolds number based on the step height and maximum inlet velocity. It becomes asymmetric as the Reynolds number is increased beyond the critical value. Our results were obtained at a Reynolds number of 56, which ensures the symmetric flow and enables us to compare the numerical results with the experimental values of Durst [5]. The geometry is simple with the inlet region constituting 1/3 of the channel height based in the middle of the channel. The geometry and boundary conditions are shown in Figure 1. Sl. 1. Enofazni tok skozi nenadno simetrično razširitev. Geometrija in robni pogoji. Fig. 1. Single-phase symmetric sudden expansion flow. Geometry and boundary conditions. Predpisani vstopni profil hitrosti je enak izmerjenemu [5] in nekoliko odstopa od paraboličnega profila razvitega laminarnega toka. Pri izstopu iz kanala smo predpisali normalne odvode hitrosti. Prikazani razultati temeljijo na diskretizaciji računskega območja z 20 x 18 podobmočji z razmerjem najdaljši/najkrajši element = 2 v koordinatnih smereh x in y. Zaradi nenadne razširitve se pojavita recirkulacijski območji ob zgornji in spodnji steni kanala, ki se z večanjem vrednosti Reynoldsovega števila povečujeta. Hitrostno polje je simetrično, kar je razvidno tudi iz eksperimentalnih rezultatov. V začetku kanala je še viden vpliv recirkulacije, medtem ko se proti izstopu vzpostavi razviti profil hitrosti. Hitrostne profile za ustaljeno rešitev smo primerjali z meritvami Dursta, Mellinga in Whitelawa [5]. The prescribed inlet velocity profile is exactly equal to the measured one [5], and differs only a little from the parabolic profile of fully developed laminar flow. At the outlet normal velocity profiles are prescribed. The discretization consists of 20 x 18 subdomains with a ratio of 2 between the longest and the shortest element. The calculated separation regions behind the expansion are, as reported by [5], of equal length, leading to the fully developed, parabolic profile far downstream. The numerical results are compared to the experimental values of Durst, Melling and Whitelaw [5]. gfin^OtJJlMISCSD stran 67 Po`arnik M. - [kerget L.: Robno obmo~na integralska metoda - Boundary Domain Integral Method Sl. 2. Enofazni tok skozi nenadno simetrično razširitev. Primerjava z eksperimentalnimi vrednostimi vzdolž toka na različnih oddaljenostih od razširitve. Fig. 2. Single-phase symmetric sudden expansion flow. Comparison with experimental data at different distances downflow the expansion. Sl. 3. Enofazni tok skozi nenadno simetrično razširitev. Primerjava z eksperimentalnimi vrednostimi vzdolž toka v sredini kanala. Fig. 3. Single-phase symmetric sudden expansion flow. Comparison with experimental data along the horizontal line through the center of the channel. Slika 2 prikazuje primerjavo izračunanih hitrostnih profilov z izmerjenimi na različnih oddaljenostih od nenadne razširitve kanala X/H = 0, 1,5, 2,5, 3,5, 5 in 10, pri čemer H označuje višino stopnice. Na sliki 3 je prikazana primerjava največjih hitrosti vzdolž simetrijske osi kanala. 4.2 Dvofazni dvosestavinski tok v navpičnem kanalu Namen raziskave je s predlagano numerično shemo ROIM simulirati vpliv različnih koeficientov medfazne izmenjave gibalne količine, kombiniranih z različnimi Stokesovimi hitrostmi trdnih delcev na hitrostni polji obeh sestavin. Geometrijska oblika in robni pogoji zastavljenega problema so prikazani na sliki 4. Pri vstopu v cev smo predpisali različne parabolične hitrostne profile. V prvem delu izračunov smo za sestavino 1 (plin) pri vstopu v kanal predpisali parabolo z največjo vstopno hitrostjo vfy maks = 1 m/s Figure (2) shows the comparison at different distances downflow the expansion X/H = 0, 1.5, 2.5, 3.5, 5 in 10, where H denotes the height of the step. In Figure 3 the maximum velocity values are shown along the centreline of the channel. 4.2 Two-phase two-component flow in the vertical channel The aim of the research was to establish with the proposed BDIM numerical scheme the influence of different drag coefficients combined with different terminal velocities of solid phase on the velocity fields of both components. The geometry and boundary conditions for the investigated two-phase two-component vertical channel flow are presented in Figure 4. Different parabolic velocity profiles were prescribed at the inlet. In the first set of calculations for the phase 1 (gas) the parabolic inlet velocity profile with vfy max = 1 m/s was defined. For the phase 2 (solid grin^(afcflM]SCLD VBgfFMK stran 68 Po`arnik M. - [kerget L.: Robno obmo~na integralska metoda - Boundary Domain Integral Method Sl. 4. Ravninski dvofazni dvosestavinski tok med navpičnima ploščama. Geometrija in robni pogoji. Fig. 4. Two-phase two-component vertical channel flow. Geometry and boundary conditions. in za sestavino 2 (trdna snov) parabolo z največjo vstopno hitrostjo vsy maks = 0,5 m/s. V drugem delu izračunov smo zmanjšali največjo vstopno hitrost trdne sestavine na vsy maks = 0,4 m/s in primerjali potek hitrostnih profilov vzdolž navpičnice skozi središče kanala. Zdrs sestavin na trdnih stenah ni bil mogoč. Pri izstopu iz kanala smo predpisali normalne odvode vektorjev hitrosti. Po celotnem računskem območju smo upoštevali realno vrednost koeficienta medfazne izmenjave gibalne količine. Izvedli smo izračune z vrednostmi od fi = 10 kg/m3s do fi = 30 kg/m3s. Ustaljeno analizo smo simulirali s prehodno za zelo velikčasovni korak (Zl= 1015). Konvergenčni kriterij vseh izračunov je bil 10-4. Sliki 5 in 6 prikazujeta navpično komponento vektorjev hitrosti vzdolž navpične črte skozi središče kanala. Oblika hitrostnih profilov je močno odvisna od koeficienta trenjaARezultati izračunov na sliki 5 kažejo odvisnost hitrosti od vrednosti koeficienta trenja/?pri izbrani Stokesovi hitrosti posameznega trdnega delca, medtem ko slika 6 prikazuje primerjavo hitrostnih profilov pri izbranem trenju med sestavinama in spreminjajoči se Stokesovi hitrosti trdnih delcev. V obeh primerih je jasno razvidno, da se z večanjem koeficienta trenja fi hitrostni profili obeh sestavin vedno bolj približujejo drug proti drugemu. Prav tako pričakovano je zmanjševanje največjih hitrosti obeh sestavin vzdolž srednice kanala zaradi povečevanja Stokesove hitrosti trdnih delcev. Simuliranja manjšajoče se hitrosti sestavine 1 in večajoče se hitrosti sestavine 2 smo za oba primera vstopnih pogojev izvedli na različnih računskih mrežah. Konvergenco rezultatov z zgoščevanjem računske mreže podajata preglednici 1 in 2. phase) the maximum inlet velocity was set to the value of vsy max = 0.5 m/s. In the second part of the calculations the maximum inlet velocity of phase 2 was decreased to the vsy max = 0.5 m/s. The velocity profiles along the vertical line through the centre of the channel were compared to each other. No-slip conditions on the rigid walls were prescribed for both phases motion. Normal velocity fluxes at the outlet were given as known quantities. All over the computational domain a realistic value of the interphase momentum transfer coefficient was given. Calculations with different values of drag coefficient between b = 10 kg/m3s and b = 30 kg/m3s were made. Steady-state analysis was simulated by a transient one for one very large time step (Dt = 1015). The convergence criteria for all runs were 10-4. Figures 5 and 6 show the vertical component of the velocity vectors along the vertical line through the centre of the channel. The shape of the velocity profiles strongly depends on the drag coefficient b. The results in Figure 5 show the influence of the drag coefficient on the velocity fields at a fixed Stokes velocity of a single solid particle. On the other hand, Figure 6 shows the results at a fixed drag coefficient and variable Stokes velocity. In both cases one can see that the velocity profiles are moving closer to each other with an increasing interphase momentum exchange coefficient. Also, a decrease of the maximum velocities due to the increasing Stokes velocity was expected. Simulations of the decreasing velocity of phase 1 and the increasing velocity of phase 2 were also carried out on different mesh densities. The convergence of the results is given in Tables 1 and 2. gfin^OtJJlMISCSD stran 69 Po`arnik M. - [kerget L.: Robno obmo~na integralska metoda - Boundary Domain Integral Method ms. ^eia 1.: ¦ Vjji teo= h i— ¦ 1\ beta= :: -»— IWLS *Wla= IS -u— . tr\ ** |i--; KI ¦ o a.f i lj 2 y :< 3 3J 1 / ¦ ' bsttFU -i— ¦ bs(aF2(l -»-hup»-'— DBbHI DQBbIB^BB i ' D O.S I 1.5 I 1.' 1 1.5 A (1 Oj ] I j 1 3.3 3 H 4 Sl. 5. Vpliv koeficienta medfazne izmenjave gibalne količine P na navpično komponento vektorja hitrosti sestavine 1 (zgoraj) in sestavine 2 (spodaj) vzdolž srednice kanala na mreži s 6 x 24 podobmočji (v k = 1 m/s, v = 0,5 m/s (levo), v ks = 0.4 m/s (desno), v„ = 10 m/s). Fig. 5. Two-phase two-component vertical channel flow. Influence of interphase momentum exchange coefficient P on vertical component of phase 1 (upper line) and phase 2 (lower line) along the vertical line through the center of the channel at discretization 6 x 24 subdomains (vinle 1 = 1 m/s, vinle 2 = 0.5 m/s (left), vinle 2 = 0.4 m/s (right), v» = 10 m/s). Sl. 6. Vpliv Stokesove hitrosti sestavine 2 v„ na navpično komponento vektorja hitrosti sestavine 1 (zgoraj) in sestavine 2 (spodaj) vzdolž srednice kanala na mreži s 6 x 24 podobmočji (v ks = 1 m/s, v k = 0,5 m/s (levo), v k = 0,4 m/s (desno), p = 25 kg/m3s). Fig. 6. Two^phase two-component vertical channel flow. Influence of Stokes velocity of phase 2 v 1 on vertical component of phase 1 (upper line) and phase 2 (lower line) along the vertical line through the center of the channel at discretization 6 x 24 subdomains (vinle 1 = 1 m/s, vinle 2 = 0.5 m/s (left), vinle 2 = 0.4 m/s (right),fi = 25 kg/m3s). grin^sfcflMiscsD VBgfFMK stran 70 Po`arnik M. - [kerget L.: Robno obmo~na integralska metoda - Boundary Domain Integral Method Preglednica 1: Ravninski dvofazni dvosestavinski tok med navpičnima ploščama. Primerjava navpične komponente izstopne hitrosti v v odvisnosti od diskretizacije (v k =1 m/s, v k = 0,5 m/s, p = 25 kg/m3s, v„ = 10 m/s). Table 1: Two-dimensional two-phase two-component vertical channel flow. Comparison of vertical component of the outlet velocity v regarding to the mesh density (vinle 1 = 1 m/s, vinle 2 = 0.5 m/s, 0= 25 kg/m3s, v« 1 --¦. ¦ i - . jltSLT J Dil^ DM -. i; Bil D.E43 Eil d.ul; -. i C_H D.B43 . - ixa D.EILJ -. !¦ - d .ut: ¦. Lil 13x4] .:i -. i- Preglednica 2: Ravninski dvofazni dvosestavinski tok med navpičnima ploščama. Primerjava navpične komponente izstopne hitrosti v v odvisnosti od diskretizacije (v ks =1 m/s, v k = 0,4 m/s, p = 25 kg/m3s, v„ = 10 m/s). Table 2: Two-dimensional two-phase two-component vertical channel flow. Comparison of vertical component of the outlet velocity v regarding to the mesh density (vinle 1 = 1 m/s, vinle 2 = 0.4 m/s, p = 25 kg/m3s, v« IIBljil,1' MniLliiL \j ElbfL^' 1 Ix l:] D.BX ii. L L! Eilf D.B34 ¦I. L Hi Eil ¦ i. -1 _; L . l_M D.B34 II.ILL Hil D.B33 II.ILL L 3 411 DBS 11.444 IJ-L Hi DBS 11.444 5 SKLEP V prispevku je prikazan razvoj robno-območne integralske metode za numerično simuliranje dvofaznih dvosestavinskih tokov plin - trdni delci. Uporabili smo hitrostno-vrtinčno formulacijo dopolnjenih Navier-Stokesovih enačb. Vodilne enačbe smo poenostavili z uvedbo predpostavke o nespremenljivem prostorninskem deležu sestavine v podobmočju v iteraciji numeričnega algoritma. Zaradi nezvezne porazdelitve hitrosti na mejah med makroelementi smo predpisali primerne pogoje vmesnega roba. Za sklenitev zapisanega sistema kinematičnih in kinetičnih enačb gibanja sestavin lebdečega sloja, smo uporabili dodatno enačbo za izračun prostorninskega deleža, ki je izpeljana iz teorije gnanega toka dvofaznih tokov. Prednost predstavljene sheme na podlagi robno območne integralske metode v primerjavi z običajnimi, ki brez izjeme temeljijo na postopkih končnih razlik, končnih elementov in kontrolnih prostornin, je zmanjšano število dodatnih modelov za določitev navideznih lastnosti trdne snovi. Pravilnost delovanja razvite sheme smo najprej potrdili na primeru enofaznega toka v kanalu z nenadno simetrično razširitvijo. Vplive 5 CONCLUSION The boundary domain integral method was used to simulate two-dimensional two-phase two-component gas-solid flow. The velocity-vorticity approach in combination with modified Navier-Stokes equations was employed. The set of governing equations was simplified under the assumption that the volume fraction of each component is constant in each macroelement within one iteration. The discontinuous velocity distribution on the interfaces between the subdomains was overcome with the prescription of the appropriate interface macroelement boundary conditions. An additional equation to compute the volume fraction of the fluid phase was obtained from the drift-flux theory. The advantage of the proposed scheme is a reduced number of gas-solid physical models. The numerical model was first validated on single-phase test examples such as single-phase symmetric sudden expansion flow. Finally, the influence of the drag force and the terminal grin^otJJiMiscsni stran 71 glTMDDC Po`arnik M. - [kerget L.: Robno obmo~na integralska metoda - Boundary Domain Integral Method koeficienta medfazne izmenjave gibalne količine in vpliv različnih Stokesovih hitrosti trdnih delcev na hitrostni polji sestavin smo testirali na primeru dvofaznega dvosestavinskega toka v navpičnem kanalu. velocity of the solid phase on both components velocity fields was tested for the case of two-phase two-component vertical channel flow. 6 SIMBOLI 6 SYMBOLS prostorninski delež sestavine ep i-ta trenutna komponenta hitrosti sestavine gostota sestavine pp i-ta koordinata snovski oziroma Stokesov odvod D/Dr tenzor viskoznih napetosti sestavine Spij gravitacijski pospešek gi termodinamični tlak p navidezni tlak trdne snovi p koeficient medfazne izmenjave gibalne količine p i-ta trenutna komponenta vrtinčnosti copi Stokesova hitrost trdnih delcev koeficient sprostitveni parameter v kinematični enačbi ap sestavine nespremenljiva kinematična strižna viskoznost sestavine nespremenljiva kinematična normalna viskoznost Kp0 sestavine permutacijski enotski tenzor eijk volume fraction of the phase i-th instantaneous phase velocity component density of the phase i-th coordinate substantial or Stokes derivative viscous stress tensor of the phase gravity acceleration thermodynamic pressure solids pressure interphase momentum transfer coefficient i-th instantaneous phase vorticity component Stokes velocity of solid particles coefficient relaxation parameter in kinematic equation of the phase constant phase kinematic shear viscosity constant phase kinematic bulk viscosity permutation unit tensor 7 LITERATURA 7 REFERENCES [1] Anderson, T.B., R. Jackson. (1967) A fluid mechanical description of fluidized beds. Ind Engng Chem Fundam 6, 527-539. [2] Bird, R.B., W.E. Stewart, E.N. Lightfoot (1960) Transport phenomena. New York. John Wiley & Sons. [3] Boemer, A., Qi, H., Renz, U., Vasquez, S., Boysan, F. (1995) Eulerian computation of fluidized bed hydrodynamics - A comparison of physical models. Proceedings of the 13th International Conference on Fluidized Bed Combustion, Orlando, 775-781. [4] Chapman, S., T.G. Cowling (1970) The mathematical theory of non-uniform gases. Cambridge, Cambridge University Press. [5] Durst, F., A. Melling, J.H. Whitelaw (1974) Low Reynolds number flow over a plane symmetric sudden expansion. J. Fluid Mechanics, 64(1), 111-128. [6] Požarnik, M. (2000) Doktorska disertacija, Univerza v Mariboru. [7] Richardson, J.F., W.N. Zaki (1954) Sedimentation and fluidization: Part I. Trans Instn Chem Engs 32, 35-53. [8] Škerget, L., Z. Rek (1995) Boundary-domain integral method using a velocity-vorticity formulation. Eng Anal Bound Elem 15, 359-370. [9] Wallis, G.B. (1969) One-dimensional Two-phase Flow. McGraw-Hill. Naslov avtorjev: dr. Matej Požarnik profdr. Leopold Škerget Fakulteta za strojništvo Univerza v Mariboru Smetanova 17 2000 Maribor matej.pozarnik@uni-mb.si leo@uni-mb.si Authors’ Address:Dr. Matej Požarnik ProfDr. Leopold Škerget University of Maribor Faculty of Mechanical Eng. Smetanova 17 2000 Maribor, Slovenia matej.pozarnik@uni-mb.si leo@uni-mb.si Prejeto: 5.9.2001 Received: Sprejeto: 29.3.2002 Accepted: 02-2»ina(B>jiMift»nl VH^tTPsDDGC stran 72 © Strojni{ki vestnik 48(2002)2,73-86 © Journal of Mechanical Engineering 48(2002)2,73-86 ISSN 0039-2480 ISSN 0039-2480 UDK 62-229:621.914:004.8:519.65 UDC 62-229:621.914:004.8:519.65 Pregledni znanstveni ~lanek (1.02) Review scientific paper (1.02) Model za analizo in optimiranje vpenjalnih priprav A Model for Analysing and Optimazing Fixtures Uro{ @uperl - Franci ^u{ V prispevku je predstavljen nevronsko-analitični model za analizo in racionalizacijo vpenjalnih priprav, primernih za vpenjanje tankostenih izdelkov, pri katerih obstaja med obdelavo velika verjetnost deformacije zaradi vpenjalnih in rezalnih sil. Izdelan je program FIXAN, ki ovrednoti vpenjalno shemo in izračuna optimalne vrednosti in položaje vpenjalnih ter podpornih sil, ki so potrebne, da je obdelovanec varno vpet med obdelavo. Model je primeren za analizo vpenjalnih priprav, namenjenih za vpenjanje prizmatičnih in rotacijsko simetričnih izdelkov. Model upošteva trenje, ki se pojavi med obdelovancem in elementi vpenjalne priprave. Hitrost izračuna je zaradi uporabe umetnih nevronskih mrež (UNM) zelo velika, zato je postopek mogoče izvesti v realnem času. Z opisanim postopkom zmanjšamo čas snovanja vpenjalne priprave in preprečimo napake in deformacije med postopkom obdelave. © 2002 Strojniški vestnik. Vse pravice pridržane. (Ključne besede: naprave vpenjalne, frezanje, optimiranje, mreže nevronalne) This paper is about a neural-analytical model for the analysis and rationalization of fixtures that are suitable for clamping thin-wall products likely to undergo deformation due to the clamping and cutting forces that occur during machining. A program called FIXAN was used for the evaluation of the fixturing scheme and for the calculation of the optimum magnitude and positioning of the clamping forces required to enable the workpiece to be safely clamped during machining. The model is suitable for the analysis of fixtures intended for the fixing of prismatic and rotational products. The model takes into consideration the friction occurring between the workpiece and the fixture components. Because of the use of an artificial neural network (ANN) the time needed for the calculation is very short, therefore, the procedure can be carried out in real time. The described procedure ensures a reduction of the fixture-planning time and the prevention of defects and deformations during the machining process. © 2002 Journal of Mechanical Engineering. All rights reserved. (Keywords: fixture analysis, milling, optimization, neural networks) 0 UVOD Snovanje vpenjalnih priprav je zapleten in domiselen postopek, ki zahteva izkušenega tehnologa. Za vsak izdelek obstaja več mogočih izvedb vpenjalnih priprav, zato je obseg mogočih rešitev velik. Razvoj umetne inteligence je prispeval k omejevanju obsega mogočih rešitev in s tem k doseganju bolj ših izvedb. Umetna inteligenca ponuja različne tehnike za rešitev problema snovanja vpenjalnih naprav. Najpomembnejši sta prigodno razsojanje (case-based reasoning) in ekspertni sistemi. Vendar ti dve tehniki ne zagotovita vedno optimalne rešitve. Bolj uspešni so sistemi, ki 0 INTRODUCTION The designing of fixtures is a complex and intuitive process for which an experienced technologist is required. For each workpiece there are several possible solutions of the design for modular fixtures, therefore, the number of possible solutions is large. The development of artificial intelligence has contributed to limiting the number of possible solutions and, consequently, to achieving better designs. Artificial intelligence offers various methods for solving the problem of fixtures’ design. The most important methods are the case-based reasoning and the expert systems. However, these two methods do not always provide the best fixture solutions. The systems that use a grin^otJJiMiscsni stran 73 glTMDDC @uperl U. - ^u{ F.: Model za analizo in optimiranje - A Model for Analysing and Optimizing uporabljajo kombinacijo več tehnik, npr. genetski algoritmi in nevronske mreže. Naloga računalniško podprtega sistema za snovanje vpenjalnih naprav [8] je izbrati pravilno kombinacijo osnovnih modularnih vpenjalnih elementov, jih postaviti in sestaviti na ustrezno mesto, tako da bo obdelovanec pri obdelavi varno vpet. Opisan sistem (sl. 1) vsebuje modul za analizo vpenjalne priprave. Zmogljivejši moduli ponujajo tudi možnost racionalizacije in optimizacije dobljene rešitve. Po navadi izračun temelji na metodi analize sil, opremljeni z optimizacijskim postopkom. combination of several methods, such as genetic algorithms and neural networks, are more successful. The task of the computer-supported system [8] for designing the fixtures is to select the correct combination of basic modular fixtures and to locate and assemble them in an appropriate place so that the workpiece will be safely fixed during machining. Such a system (Figure 1) contains a module for the analysis of the fixture. These highly capable modules also offer the possibility to rationalize and optimize the obtained fixture solution. Usually, the calculation is based on the force analysis method provided with the optimization process. Tehnološka baza podatkov Technological database OBLIKOVANJE IZDELKA (CAD) PRODUCT DESIGN (CAD) NAČRTOVANJE TEHNOLOGIJE DESIGNING OF TECHNOLOGY Baza pravil optimalne izbire Selection & optimal rule database I SNOVANJE VPENJALNE PRIPRAVE CONCEIVING OF FIXTURE Osnovna plošča Podporni element Vpenjalni element Baseplate Locator Clamp |ne ustreza no pass I Baza vpenjalnih komponent Assembly fixture database MONTAŽA VPENJALNE PRIPRAVE FIXTURE ASSEMBLY Ne/No Sl. 1 Zgradba avtomatiziranega sistema za izbiro analizo in montažo vpenjalnih priprav Fig. 1. Structure of the automated system for the selection, analysis and assembly of fixtures 1 DOSEDANJE RAZISKAVE IN UGOTOVITVE Poleg iskanja matematičnih rešitev za pozicioniranje in vpenjanje obdelovancev poteka razvoj v smeri iskanja rešitve z uporabo računalniške rutine ([4] in [5]) Raziskovalci so nedavno predlagali model “obdelovanec-vpenjalna priprava” ki temelji na osnovi “screw theory” in uporabili , metodo linearnega programiranja za določitev vpenjalnih sil [3] Mittal predlaga dinamični model “vpenjalna priprava- 1 CURRENT RESEARCH AND FINDINGS In addition to searching for a mathematical solution for the positioning and clamping of workpieces the development is oriented towards searching for solutions by means of a computer routine ([4] and [5]). Recently, researchers [3] proposed a “workpiece-fixture” model based on the screw theory and used the linear programming method for a determination of the clamping forces. Mittal proposed the dynamic “fixture-workpiece” model to determine 02-2^rsuWI^r W§©™0C€ stran 74 @uperl U. - ^u{ F.: Model za analizo in optimiranje - A Model for Analysing and Optimizing obdelovanec” za določitev potrebnih vpenjalnih sil, ki so potrebne, da je obdelovanec med obdelavo v ravnotežju [1]. Vse naštete metode uporabljajo poenostavljene modele, ki ne upoštevajo trenja pri svojih izračunih. Rezultati so le približni. 2 TEORETIČNA IZHODIŠČA PRI IZDELAVI MODELA Pri izdelavi modela je bilo predpostavljeno, da bo obdelovanec vpet s prilagodljivo modularno vpenjalno pripravo, ki omogoča vpenjanje obdelovancev različnih oblik. Predpostavljeno je, da sta obdelovanec in vpenjalna priprava togi telesi. 2.1 Vpenjalno načelo Pri izdelavi modela je bilo uporabljeno za vpenjanje vpenjalno načelo 3-2-1 (sl. 2), ki zahteva tri lokatorje na osnovni položajni ravnini, dva lokatorja na drugi položajni ravnini in en lokator na tretji položajni ravnini. the required clamping forces needed for equilibrium of the workpiece during machining [1]. All the above-mentioned methods use simplified models that do not take friction into account in their calculations. The results are only approximate. 2 A THEORETICAL ASSUMPTION USED IN MAKING THE MODEL When making the model it was assumed that the workpiece would be fixed by a flexible modular fixture to ensure that workpieces with different shapes could be clamped. The fixture and the workpiece are assumed to be rigid bodies. 2.1 Clamping principle The clamping principle 3-2-1, which requires three locators on the first locating plane, two locators on the second locating plane and one locator on the third locating plane (Figure 2), was used for making the model. Sl. 2. Uporabljen vpenjalno načelo 3-2-1 za vpenjanje prizmatičnih obdelovancev Fig. 2. The 3-2-1 clamping principle used for clamping the prismatic workpieces 2.2 Vpliv vpenjalnih sil na shemo vpetja Vpenjalna sila mora biti dovolj velika in ustrezno usmerjena, da se lega obdelovanca med obdelavo zaradi rezalnih sil ne spremeni. Vpenjalne sile na obdelovancu ne smejo ustvarjati notranjih napetosti in poškodovati oziroma deformirati površine obdelovanca. Dovoljene deformacije smejo biti takšne, da je izdelek po postopku obdelave znotraj predpisanih tolerančnih vrednosti. Vpenjalne sile delujejo vedno proti podporam, na katerih je podprt obdelovanec. 2.3 Rezalne sile Izračun rezalnih sil je izveden po nevronskem modelu za simuliranje rezalnih sil [7]. Analitično modeliranje rezalnih sil [2] je težavno zaradi velikega števila medsebojno odvisnih obdelovalnih parametrov. Namesto poskusov iskanja analitičnih povezav med obdelovalnimi 2.2 The influence of clamping forces on the clamping scheme The clamping force must be large enough and suitably oriented so that the workpiece’s position does not change during machining due to the cutting forces. The clamping forces on the workpiece’s are not allowed to create internal stresses and to damage or deform the workpiece surface. The permissible deformation may be such that after machining the product is within the specified tolerance values. The clamping forces always act towards the locators on which the workpiece is supported. 2.3 Cutting forces The calculations of the cutting forces are made according to the neural network model for the simulation of cutting forces [7]. Analytical cutting-force modelling [2] is difficult due to the large number of interrelated machining parameters. Instead of attempting to find analytical relationships between gfin^OtJJlMISCSD stran 75 @uperl U. - ^u{ F.: Model za analizo in optimiranje - A Model for Analysing and Optimizing parametri s statistiko, je uporabljeno strojno učenje. Izdelava modela temelji na izvedenih eksperimentalnih meritvah rezalnih sil pri frezanju. Za modeliranje rezalnih sil so bile uporabljene usmerjene nevronske mreže (UNM) s tremi nivoji. Vsebovale so 11 nevronov v vhodni ravni in 3 nevrone v izhodni ravni. Pri preskusih se je število nevronov v skritem nivoju spreminjalo. Učenje UNM je bilo izvedeno z naslednjimi parametri: material obdelovanca, trdota obdelovanca, premer orodja, tip ploščice, rezalna hitrost, podajanje, prečna in vzdolžna globina rezanja, obraba orodja, komponente rezalne sile. Med postopkom učenja so bili UNM posredovani tudi želeni izhodi (tri komponente rezalne sile). Učenje UNM je bilo izvedeno z neobdelanimi eksperimentalnimi podatki, pri čemer je bilo uporabljenih 3500 popolnih učnih primerov. 3 MODEL ZA ANALIZO IN OPTIMIRANJE VPENJALNIH PRIPRAV Razvit model je uporaben pri snovanju vpenjalnih priprav, saj lahko v kratkem času rutinsko določi optimalne velikosti, smeri in prijemališča vpenjalnih in podpornih sil za različne vpetostne primere. Naloga modela je, da preveri (analizira) dobljeno rešitev (konfiguracijo vpenjalne priprave), jo potrdi oziroma zavrne, če niso izpolnjeni vsi zastavljeni pogoji. Namen modela je izboljšati izvedbo vpenjalne priprave in s tem povečati geometrijsko natančnost izdelanega tankostenega izdelka. Pomembno je upoštevati rezalne sile in vpenjalne sile ter izmere in razpoložljivost vpenjal, kakor tudi delovni prostor na stroju, ki omejuje možnosti vpetja. Pri izdelavi programa FIXAN je bilo znanje izbrano iz literature (proizvajalci vpenjal) in od izvedencev iz prakse v proizvodnji. Pomemben vir znanja so že izdelani in preverjeni tehnološki postopki. Stroški za načrtovanje in izdelavo vpenjalne priprave znašajo tudi do 15% celotnih proizvodnih stroškov [6]. Nadalje, če se pozicioniranje in fiksiranje izdelka ne da izvesti z zmernimi stroški, ali ne zadostimo zahtevam postopka, potem postopek obdelave ni upravičen. Zmanjševanje stroškov in časa za načrtovanje postopka vpenjanja je največje gonilo za sistematično načrtovanje vpenjalnih priprav. Program FIXAN določi (sl. 3): - minimalno število in položaj podpornih in vpenjalnih elementov, - gibanje, ki ga dovoljujejo podporni elementi, - reakcije na mestih stika “obdelovanec-vpenjalna priprava” (podporne sile), - minimalne vpenjalne sile, potrebne za uravnoteženje rezalnih sil. 02-2»«firaByi*!UBii«iil Mig™« | stran 76 the machining parameters using statistics, machine learning is used. The development of the model is based on the measurements of cutting forces during milling. For modeling the cutting forces, three-layer feed-forward artificial neural networks (ANNs) were used. They contained 11 neurons in the input layer, and 3 in the output layer. The number of neurons in the hidden layer was varied in the experiments. The ANNs were trained with the following parameters: type of machined material, hardness of the machined material, cutting tool diameter, type of insert, cutting speed, feed, radial and axial depth of cutting, tool wear, and components of cutting force. The desired outputs (the three cutting-force components) of the network also being supplied during training. Training of the ANN was made with raw experimental data from 3500 full training examples. 3 PRESENTATION OF THE MODEL FOR THE ANALYSIS AND OPTIMIZATION OF FIXTURES The developed model is useful for designing fixtures since it can routinely determine, within a short time, the optimum sizes, the direction and the application points of clamping and locating forces for different cases of clamping. The model is aimed at verifying (analysing) the obtained solution (configuration of the fixture), confirming or rejecting it, if all the set conditions are not fulfilled. The purpose of the model is to improve the design of the fixture and thus to increase the geometrical accuracy of the thin-wall product. It is important to consider the cutting forces, the clamping forces and the dimensions and availability of fixtures as well as the space on the machine, which limits the possibility of clamping. For working out the programme FIXAN, information was collected from the literature (makers of modular fixtures) and from experts in production practice. The technological procedures already worked out and verified are an important source of information. The designing and manufacturing costs of the fixture amount to as much as 15% of the total production costs [6]. Further more, if the positioning and fixing of the product cannot be carried out with moderate costs or if the requirements of the process are not satisfied, the machining process is not justified. The reduction of the costs and the time to design the fixing process is the greatest motivating force for the systematic design of fixtures. The FIXAN program determines (Figure3): - the minimum number and position of the locating and clamping elements, - the motion allowed by the locating elements, - the reactions at the places of the “workpiece-fix-ture” contact (locating forces), - the minimum clamping forces required for balancing the cutting forces. @uperl U. - ^u{ F.: Model za analizo in optimiranje - A Model for Analysing and Optimizing PODPORNIH ELEMENTOV of LOCATING COMPONENTS VPENJALNIH ELEMENTOV of CLAMPING COMPONENTS ^ :r t NAJPRIMERNEJŠI POLOŽAJI MOST SUITABLE POSITIONS VPENJALNE SILE PROGRAM FIXAN PODPORNE SILE LOCATING FORCES CLAMPING FORCES 1' POTREBNO ŠTEVILO REQUIRED NUMBER r ^ PODPORNIH ELEMENTOV of LOCATING COMPONENTS VPENJALNIH ELEMENTOV of CLAMPING COMPONENTS Sl. 3. Shematski prikaz nalog v programu FIXAN Fig. 3. Representation of the tasks in the FIXAN program 3.1 Teoretična zasnova modela za analizo in 3.1. Theoretical concept of the model for the analysis optimiranje vpenjalnih priprav and the optimization of fixtures Stvari postanejo močno zapletene, kadar upoštevamo sile trenja med obdelovancem in elementi vpenjalne priprave. Obdelovanec je podprt na šestih točkah P -P6 in vpet s tremi vpenjalnimi silami Fvp1, Fvp2, Fvp3 v točkah P, P, P (sl. 4). Kjer so: F1,...,F6 - reakcije, ki delujejo na podporne elemente (N), Fvp1, Fvp2, Fvp3 -vpenjalne sile, ki delujejo v smeri normale na pozicionirne ravnine(N), Ft, F, Ff - komponente rezalne sile F (N), M, My M - komponente rezalnega momenta M (Nm), f, (i = 1...9) - rezultirajoče sile trenja v stičnih točkah (N), F - sila teže obdelovanca (N), /g - koeficient trenja. Fa, Mz The scene is much more complex when frictional forces between the workpiece and fixture elements are taken into account. The workpiece is located on the six points P1-P6 and is held by three clamping forces Fvp1, Fvp2, Fvp3 at points P7, P8, P9 (Figure 4). Where: F1,.,.,F6 - reactions acting on locating elements (N), Fvp1, Fvp2, Fvp3 - clamping forces acting in the direction of the normal on to positioning planes (N), Ft, F, F - components of cutting force F (N), M, My M - components of cutting moment M (Nm), fi, (i = 1...9) - resulting frictional forces in contact points (N), F - force of workpiece weight (N), g - friction coefficient. Mc, Fc F5 Zi Y F2 Sl. 4. Sile na obdelovancu med postopkom frezanja Fig. 4. Forces on the workpiece during the milling process @uperl U. - ^u{ F.: Model za analizo in optimiranje - A Model for Analysing and Optimizing Rezultirajoča sila trenja fi med lokatorjem in obdelovancem je u .Fi in med vpenjalnim elementom in obdelovancem ^ .Fvp j, (j = 1 ..3). Reakcije na podpornih elementih morajo biti pozitivne, ker se v nasprotnem primeru izgubi stik med obdelovancem in elementi vpenjalne priprave. 3.2 Ravnotežne enačbe Za dosego statičnega ravnotežja in dimenzijske natančnosti pri obdelavi morata biti rezultirajoča sila in moment na obdelovancu enaka nič. Ravnotežne enačbe so: The resulting force of the friction fi between the locator and workpiece is m .Fi and between the fixing element and the workpiece it is m .Fvp j, (j = 1...3). The reactions on the locating elements must be positive, because otherwise the contact between the workpiece and the fixture is lost. 3.2 Equations of equilibrium To achieve static equilibrium and dimensional accuracy during machining the resulting force and moment on the workpiece must be zero. The equations of equilibrium are: XFi -Rx= ZFi]-Ry= ZFi -Rz=0 E(Fixri) \-Mx= L(Fi My=\ll(F^rM -Mz=0 (1) kjer so: x y where: r - vektorji, ki definirajo podporne točke, ri - the vectors defining the locating points Rx, R, R - komponente rezultirajoče rezalne sile F. rx, Ry R - components of the resulting cutting force F 3.3 Matrične ravnotežne enačbe za izračun 3.3 Matrix equilibrium equations for the calculation podpornih sil of locating forces Zaradi numeričnega reševanja problema so Because the problem is solved numerically the ravnotežne enačbe zapisane v matrični obliki: equilibrium equations are written in the matrix form: Geometrijska matrika [A] je enaka: lok [] [] +we ] =0 The geometrical matrix [A]lok is as follows: (3) r f1 f 1y [A — -r1 L 2x f 3x 2y f 3y 1 1 -1 f 4y -f4z -1 f 5y -f5z f 6x -f6z 2y f 4y ' r z V f 4y ' r 4y J -r4z +f4z -f .r -\ f-f .r + N 5y 5z 6z 6y -f5z-r5 V 5z 5y 7 +f5z f1x-r1y + y+ f1y' r1x j -f2x-r2y + \+f2y ' r2x J f3x ' r3y + \+f3y ' r3x J r4y + r5y + \+f4y 'r4x J \+f5y 'r5x J f6x ' r6z + \+f6z ' r6x / f -r \ 6x 6y (4) Vektor podpornih sil [F] T: lok Vektor zunanjih sil [w ]: we [] The vector of supporting forces [F]lokT: [F]lokT =[F1 F2 F3 F4 F5 F6 ] The vector of external forces [we]: f7x + f9x +Fvp2 +Rx f7 y + f8 y +Fvp +Ry 3 -f8z - f9z -Fg +Rz (5) -f7y f 8y f8 f 9z z 8y 9z 9y F -r F - -F -r +M vp3 r9z g gyx f7 -r7 + f ¦ r + f9 -r9 + f9 -r9 +F -r7 +F ¦ r +F -r +M xz 8z 8x x z z x vp 1x vp 2 8z g gx y -f7x ' r7y + f7y ' r7x + f8y f 9x ' r 9y Fvp 2-r8y+Fv 3-r9x+Mz (6) grin^(afcflM]SCLD VBgfFMK stran 78 @uperl U. - ^u{ F.: Model za analizo in optimiranje - A Model for Analysing and Optimizing Po vstavitvi geometrijske matrike [A], in After entering the geometrical matrix [A]lo and vektorja zunanjih sil [w ] v enačbo (3) sledi: ok the vector of external forces [w ] into equation (3) the following is obtained: r f1 f 1y 1 — 2x f 3x 2y f 3y 1 1 -1 f 4y -f4z -1 f 5y -f5 5z f 6x -1 -f6z -r 2y f 4y ' r4z f5y ' r5z f 6z ' r6y + Vf 4y ' r 4y J 5y 5z V 5z'r 5y J -r -r 1x 2x 3x i f +f5z y+ f1y ' r1x J I 2x 2y \+f2y ' r2x J 3x 3y \+fy ' r3x J 4y \+f4y 'r4xJ r5y + \+f5y 'r5x J f6x ' r6z + \+f6z ' r6x f -r \ 6x 6 y (7) f7x + f9x +Fvp2 +Rx f7y + f8y +Fvp +Ry 3 -f8z - f9z -Fg +Rz f7y'r7z f8y'r8z f8z'r8y f9z'r9y Fvp1'r7y vp3'r9z g ' rgy + x 1 f7x ¦ r7z + f8z ¦ r8x + f9x ¦ r9z + f9z ¦ r9x + Fvp1 ¦ r7x + F ¦ r8z +F -r +M L -f7x ¦ r7y + f7y ¦ r7x + fy ' r8x - f9x ' r9y - Fvp2 ' r8 y + Fvp3 ¦ r9 x + Mz J Kadar je koeficient trenja med obdelovancem in vpenjalnimi elementi enak nič, se zgornja enačba poenostavi in dobi naslednjo obliko: [F]lo When the coefficient of friction between the workpiece and the clamping elements is equal to zero, the above equation is simplified and assumes the following form: (8) fA ' 0 0 0 -1 -1 0 -1 Fvp 2 + Rx F2 0 0 0 0 0 -1 Fvp3+Ry F3 = - 1 1 1 0 0 0 Fg+Rz F4 r 1y r 2y y 3y 0 0 +r6y -F 1 r7y-F 3 r9 z-Fg-rgy+Mx F5 -r1x -r2x -r3x -r4z -r5z 0 F 1-r7x+Fvp 2-r8z+Fg-rgx+My F6] 0 0 0 r 4y r 5y -r6x _ Fvp2-r8y+Fvp3-r9x+Mz 3.4 Iskanje ustrezne vpenjalne razporeditve in vpenjalnih sil Zaradi sil trenja je število neznank v sistemu večje od števila ravnotežnih enačb, zato sistem ni vedno rešljiv. Sistem ima netrivialno rešitev takrat, ko je determinanta matrike različna od nič. Postopek izračuna vpenjalnih sil se poenostavi z iteracijskim načinom reševanja matrične enačbe. Z iteracijsko metodo rešimo enačbo (3), s čimer izračunamo minimalne potrebne vpenjalne sile. Iteracija se začne z začetno vrednostjo vpenjalnih sil Fvp j = 0; j = 1, 2, 3, potem se ta vrednost postopoma koračno veča, dokler vse sile Fi niso pozitivne. Na tak način pridemo do osnovne rešitve problema. Osnovno rešitev je mogoče optimirati takole: Prvi vpenjalni sili se priredi vrednost osnovne rešitve, preostalim pa se postopoma koračno zvišuje 3.4 Searching for the appropriate fixing configuration and clamping forces Because of frictional forces the number of unknown variables in the system is far larger than the number of equilibrium equations, therefore, the system is not always solvable. The system has a non-trivial solution when the determinant of the system is other than zero. The procedure for calculating the clamping forces is simplified by the iteration method of solving the matrix equation. By using the iteration method equation (3) is solved and the minimum required clamping forces are calculated. The iteration starts with the initial value of the clamping force Fvp j = 0; j = 1, 2, 3, afterwards, this value gradually increases incrementally until all the forces Fi are positive. In this way we reach the basic-fundamental solution of the problem. The basic solution can be optimised in this following way: The value of the basic solution is adapted to the first clamping force whereas the values of the grin^otJJiMiscsni stran 79 glTMDDC @uperl U. - ^u{ F.: Model za analizo in optimiranje - A Model for Analysing and Optimizing vrednost, dokler vse izračunane podporne sile niso pozitivne. Nato se postopek ponovi za vsako vpenjalno silo. others are gradually increased incrementally until all the calculated locating forces are positive. The procedure is then repeated for each clamping force. 3.5 Grafični prikaz poteka dela v programu FIXAN 3.5 The graphical representation of the process of work in the FIXAN program Slika 5 prikazuje zaporedje korakov, ki jih je Figure 5 shows the sequence of steps to be taken treba izvesti pri iskanju optimalne vpenjalne sheme. when searching for the optimum clamping scheme. J Geometrija 1) obdelovanca 1 Workpiece geometry Vrednost in položaj ^ rezalne sile 2 Value and position of cutting force Koeficient trenja Friction coefficient I Začetni položaji vpenjalnih in podpornih elementov 'Starting positions of clamping and locating components Napovedovanje rezalnih sil z nevronsko mrežo Prediction of cutting forces by neural network Da / Yes GRAFIČNI PRIKAZ REZULTATOV GRAPHICAL REPRESENTATION OF RESULTS -vpenjalne sile / clamping forces podporne sile / locating forces -potrebno število vpenjalnih / podpornih elementov required number of clamping / locating elements -najprimernejši položaji : vpenjalnih / podpornih elementov most suitable positions: clamping / locating elements__________________________________ Sl. 5. Definicija programskih korakov v programu FIXAN Fig. 5. The definition of the program steps in the FIXAN program grin^sfcflMiscsD VH^tTPsDDIK stran 80 3 4 @uperl U. - ^u{ F.: Model za analizo in optimiranje - A Model for Analysing and Optimizing Program zahteva v korakih (1 do 4) ročen vnos začetnih (želenih) položajev vpenjalnih/ podpornih elementov ter vrednosti rezalnih sil, koeficienta trenja in teže obdelovanca. Začetne položaje komponent priprave določi operater glede na izkušnje, geometrijo obdelovanca in operacijo obdelave. V primeru nerodne izbire položajev komponent vpenjalne priprave, ko program z variiranjem vrednosti vpenjalnih sil ni zmožen zagotoviti učinkovite vpenjalne sheme, avtomatično spremeni koordinate posameznim vpenjalnim ali podpornim elementom. Elemente priprave razporedi okrog obdelovanca tako, da je ta varno vpet z minimalnimi potrebnimi vpenjalnimi silami. 4 PRIMER ANALIZE IN OPTIMIRANJA VPENJALNE SHEME Na številsko krmiljenem frezalnem stroju je treba izdelati utor, prikazan na sliki 6. Uporabimo frezalo s premerom 16 mm z dvema rezalnima ploščicama (R-216-16 03 M-M) pri naslednjih rezalnih pogojih: rezalna hitrost (v = 25 m/min), podajanje na zob (f = 0,01 mm/zob), globina rezanja (a = 4 mm). Material obdelovanca je jeklo z oznako Ck-45. Na podlagi formul, podanih v literaturi [7], so izračunane komponente rezalnih sil pri podanih rezalnih parametrih. Vrednosti komponent rezalnih sil (F = 450 N, F = 315 N, F = 810 N), položaj orodja, začetne položaje vpenjalnih/podpornih elementov, koeficient trenja (h = 0,4), in težo obdelovanca (Fg = 47 N) vstavimo v okno za vnos vhodnih podatkov. 4.1 Določitev vpenjalnih sil z iteracijo Vpetje izvedemo s tremi vpenjalnimi elementi. Z zgornjim vpenjalnim elementom je vpet In steps (1 to 4) the programme requires manual entering of the starting (desired) positions of the clamping/locating elements and the values of the cutting forces, the friction coefficient and the workpiece’s weight. The starting position of the fixture components are determined by the operator on the basis of experience, workpiece geometry and machining operation. In the case of inappropriate selection of the fixture components, when the program is not capable of ensuring an efficient clamping scheme by varying the values of the clamping forces, it automatically changes the coordinates of the individual clamping or locating elements. It arranges the clamping components around the workpiece so that the latter is safely clamped by the minimum required clamping forces. 4 EXAMPLE OF AN ANALYSIS AND OPTIMIZATION OF THE CLAMPING SCHEME On a NC milling machine it is necessary to make the slot shown in Figure 6. To this end we use a milling cutter of 16-mm diameter with two cutting inserts (R-216-16 03 M-M) and the following cutting conditions: cutting speed (v = 25 m/min), feedrate (fz = 0,01 mm/tooth), cutting depth (a = 4 mm). The workpiece material is the steel Ck-45. The components of the cutting forces with the cutting parameters defined above are calculated on the basis of equations given in literature [7]. The values of the components of the cutting forces (Fa = 450 N, Ff = 315 N, Ft = 810 N), the tool position, the starting positions of the clamping/locating elements, the friction coefficient (h = 0.4) and the workpiece’s weight (Fg = 47 N) are entered into the window for the input data. 4.1 Determination of clamping forces by iteration Clamping is effected by the three clamping elements. With the upper clamping element the Fvp Zl Y X F2 polozaji podpornih elementov: positions for locators: 1 (38, 140, 0) 2 (165, 80,0) 3 (165, 241, 0) 4 (178, 102, 76) 5 (254, 188, 76) 6 (102, 254, 76) polozaji vpenjalnih elementov: positions for clamping elements: 1 (122, 140, 100) 2 (25, 140,76) 3 (102, 25, 76) polozaj rezalne sile: position of cutting force: (76, 127, 96) Sl. 6. Vpenjalna in podporna shema za izdelavo utora na prizmatičnem obdelovancu Fig.6. Clamping and locating scheme for machining the slot on a prismatic workpiece @uperl U. - ^u{ F.: Model za analizo in optimiranje - A Model for Analysing and Optimizing obdelovanec v smeri osi Z, s stranskima pa je pritisnjen ob navpični pozicionirni ravnini. Enačba (3) dobi naslednjo obliko: workpiece is clamped in the direction of the Z axis and with the two side elements it is pressed along the vertical locating plane. Equation (3) assumes the following form: 0,373 0,373 0,373 -1 -1 0,145 0,145 0,145 0,229 0,229 1 1 1 -0,328 -0,328 5,512 1,496 9,488 -2,002 -3,112 -1496 -6,496 -0,696 -0,696 0,385 -1,838 0,384 -2,595 5.623 9,695 0,194 -1 -0,35 -0,505 1,985 -5,958 F1 F (9) =- 0,373-Fvp+Fvp+0,194-Fvp+810 0,145-Fvp+0,328-Fvp+Fvp+315 -0,459-Fvp -6.38-Fvp-2,245-Fvp-7,034-Fvp+3,218-Fvp+1,985-Fvp+6361,4 + 359,05 2 -F 1 -450 + 80 3,336-F -440,95 + 4501,57 1,359 -F vp 1 5,321-F +3,825-F 2910,63 Mogoče rešitve za Fvpi so tiste, ki vodijo do pozitivnih vrednosti Fi; z drugimi besedami, obdelovanec bo ostal med obdelavo v stiku s podpornimi elementi. Determinanta Floc je negativna, zato obstaja netrivialna rešitev problema. Obstoj netrivialne rešitve pomeni, da je konfiguracija vpenjalne priprave sprejemljiva. Za rešitev šestih simultanih linearnih ravnotežnih enačb z devetimi neznankami predpostavimo, da imajo Fvp1, Fvp2 , in Fvp3 enake vrednosti. Vrednosti teh treh spremenljivk so na začetku enake 0, nato se z nespremenljivim korakom v vsaki iteraciji stopnjujejo do vrednosti, pri katerih so vse sile F pozitivne. Dobljene vrednosti sil Fvp1, Fvp2 in Fvp3, ki so enake 440N, bodo prvi niz sprejemljivih rešitev, ki so navedene kot primer (1) v preglednici 1. The possible solutions for Fvpis are those that result in positive values of Fis; in other words, the workpiece will remain in contact with the locators during the entire cutting process. The determinant of Floc is other than zero, therefore, a non-trivial solution exists. The existence of the non-trivial solution implies that the clamping fixture configuration is acceptable. To solve these six linear simultaneous equilibrium equations with nine unknowns, we assume that Fvp1, Fvp2 and Fvp3 have the same magnitude. Their values start from zero and are incremented by a constant value in each iteration until positive values of all Fis are achieved. The obtained values of Fvp1, Fvp2 and Fvp3, which are equal to 440N, will be the first set of possible solutions, which are listed as case (1) in Table 1. Preglednica1. Reakcije na podpornih elementih, ki so posledica različnih načinov vpetja obdelovanca Table 1. Reaction on the locating elements which results from different methods of workpiece clamping Sile N Forces N Fvp1 Fvp2 Fvp3 F1 F2 F3 F4 F5 F6 (1) primer / case 440 440 440 2,14 475,9 693,7 826,0 673,5 962,9 (2) (3) (4) (5) 430 440 410 0 440 440 220 320 440 0 0 0 0,967 125,57 472,12 432,77 2 688,765 613,25 419,354 643,647 823,769 652,30 539,109 618,558 672,009 761,73 590,891 564,291 961,526 522,97 419,862 446,534 7,808 0,44 249,644 396,842 grin^sfcflMiscsD VH^tTPsDDIK stran 82 @uperl U. - ^u{ F.: Model za analizo in optimiranje - A Model for Analysing and Optimizing 4.2 Racionalizacija zasnove vpenjalne priprave 4.2 Rationalization of fixture design Preglednica 1 prikazuje vrednosti reakcij Fi (i = 1 do 6) na podpornih elementih in mogoče rešitve za Fvp pri različnih primerih vpetja. Primer (1) prikazuje prvo sprejemljivo rešitev, pri kateri so vse vpenjalne sile enake (440 N). Pri drugih primerih (2 in 3) priredimo dvema vpenjalnima silama vrednost 440 N, medtem ko vrednost tretje vpenjalne sile postopoma večamo, dokler ne postanejo vsi Fi pozitivni. S takim postopkom zmanjšamo vrednost potrebnih vpenjalnih sil, pri čemer pa ostane obdelovanec vseskozi v stiku z elementi vpenjalne priprave. Vrednosti vpenjalnih sil, ki jih dobimo v primeru (5), so bolj še kakor v primeru (1). V primeru (3) vpenjalni element 3 ni potreben, saj je sila na ta element enaka 0 (Fvp = 0). Program vedno predlaga dve optimalni rešitvi konfiguracije vpenjanja (primer 4 in 5), med katerima lahko uporabnik izbira. Vpliv vpenjalnih sil Fvp1, Fvp2 in Fvp na vrednosti reakcij podpornih elementov med frezanjem prikazuje slika 7. Z večanjem vpenjalne sile, po pričakovanjih, se večajo tudi vrednosti reakcij na podpornih elementih. Kritično mesto vpenjalne priprave je pri podpornem elementu 1. Pri obdelavi bo obdelovanec najprej zgubil stik z vpenjalno pripravo prav na tem mestu (P1). Razpored vpenjalne priprave je mogoče izboljšati: - s postavitvijo dodatnega vpenjalnega elementa nad kritični podporni element 1, - povečati velikost vpenjalne sile Fvp, ali spremeniti položaj vpenjalnega elementa 7. Table 1 lists the reaction forces Fi (i = 1 to 6) on the six locators and the possible solutions for Fvpis under different clamping conditions. The case (1) shows the first acceptable solution where all the clamping forces are equal (440 N). For the other cases (case 2 and 3) two of the clamping forces are set to a value of 440N, while the value of the third clamping force is gradually increased until all Fis are positive. With such a procedure the value of the required clamping forces is reduced, while the workpiece remains in contact with the fixture components at all times. The values of the clamping forces obtained in case (5) are more adequate than in case (1). In case (3) the clamping element 3 is not necessary, since the force acting on that element is equal to 0 (Fvp3 = 0). The program proposes the optimum solutions of the clamping configuration (case 4 and 5) from which the user can choose. The influence of the clamping forces Fvp1, Fvp2, Fvp3 on the values of the reactions of the locating elements during milling is shown in Figure 7. If the clamping forces are increased, the values of the reactions on the locating elements are also expected to increase. The critical point of the fixture is on the locating element 1. During machining the workpiece will first lose contact with the clamping device just at that point (P1). The fixture configuration can be improved by: - placing on additional clamping element above the critical locating element 1, - increasing the value of the clamping force Fvp1 or changing the position of the clamping element 7. Sl. 7. Prikaz reakcij na lokatorjih v različnih razmerah vpenjanja Fig. 7. Representation of the reactions of locators under different clamping conditions @uperl U. - ^u{ F.: Model za analizo in optimiranje - A Model for Analysing and Optimizing Pri vpenjalni sili Fvp =0N do Fvp =440N konfiguracija vpenjalne priprave v danih razmerah obdelave ni primerna za vpenjanje obdelovanca. Obdelovanec bo varno vpet med obdelavo takrat, ko bo vpenjalna sila večja ali enaka 440N (odebeljena črta na sliki 7). Največja obremenitev vpenjalne naprave je na mestu podpornega elementa 6. Ta prevzame med obdelavo v položajni ravnini XZ vse obremenitve, ki delujejo v smeri osi Y. Predlagane in analizirane vpenjalne primere (1 in 4) prikazuje slika 8. With clamping force Fvp=0N up to Fvp=440N the clamping fixture configuration is not suitable for clamping the workpiece in the given conditions of machining. The workpiece will be safely clamped during machining, when the clamping force is higher than, or equal to, 440N (thickened line in Figure 7). The maximum loading of the clamping fixture occurs at the point of the locating element 6, which, during machining in the locating plane XZ, takes all loadings acting in the direction of the Y axis. The proposed and analysed clamping cases (1 and 4) are shown in Figure 8. Sl. 8. Primerjava vpenjalnih primerov Fig. 8. Comparison of clamping cases S programom FIXAN je mogoče s primerno izbiro in postavitvijo vpenjalnih elementov izboljšati razpored vpenjalne priprave. 5 ANALIZA REZULTATOV MODELA S testiranjem je bila potrjena pravilnost rezultatov modela. Odstopanje napovedanih sil od dejanskih je nekoliko večje le pri zelo majhnem koeficientu trenja (0,01 0,3. Povprečna napaka napovedi je A five-layer feed-forward neural network was used. It contains 18 neurons in the input layer, and 9 in the output layer. The input vector consists of: components of the cutting forces, the coordinates of the point of machining, the coordinates of the position of the clamping and supporting parts, the workpiece’s weight and the friction coefficient. The output vector contains nine corrections factors by which the values of the calculated forces are multiplied. The desired outputs (the nine clamping/ locating forces) of the ANN are also supplied during training. Training of the ANN was made with the experimental data of 80 training examples. An additional 40 examples were used to test the trained network. The data for training and testing are obtained from the experimental measurements on the fixtures already made. Due to the introduction of the ANN the accuracy of the predicted forces was improved by 32% in the case of fi < 0.3 and by 2% in the case of ft > 0.3 . The average estimation error is about 7.4%, gfin^OtJJIMISCSD stran 85 @uperl U. - ^u{ F.: Model za analizo in optimiranje - A Model for Analysing and Optimizing približno 7,4%, kar je malo v primerjavi z 12,7-odstotno napako pri analitičnem modelu. Opisan postopek z uporabo UNM je hiter, preprost in učinkovit. 6 SKLEP Z razvitim modelom smo pomembno skrajšali čas snovanja priprav (15%) ter dosegli večjo izdelovalno natančnost. Z opisanim sistemom je mogoče napovedati in preprečiti napake na obdelovancu med postopkom vpenjanja in obdelave. V raziskavi je ugotovljeno, da se z upoštevanjem trenja močno zmanjša vrednost potrebne vpenjalne sile kakor tudi število potrebnih vpenjalnih elementov. Razvit program FIXAN skrajša čas načrtovanja postopka obdelave. Omogoča izdelavo kakovostnih načrtov vpenjanja tudi manj izkušenim operaterjem. which is low compared to the 12,7% estimation error of the analytical model. The described procedure with the use of an ANN is fast, simple and efficient. 6 CONCLUSION With the developed model we have significantly reduced the time of producing the fixture (15%) and we have reached a greater manufacturing accuracy. With the described system it is possible to anticipate and prevent any defects on the workpiece during the clamping and machining process. During the research we found that by taking the friction into account the value of the required clamping force as well as the number of the required clamping elements decreased significantly. The FIXAN program reduces the planning time for the machining process; it even allows inexperienced operators to prepare high-quality clamping drawings. 7 LITERATURA 7 REFERENCES [1] Mittal, R.O. (1991) Dynamic modeling of the fixture-workpiece system. Robotics Computer-Integr Mfg 8, London, 201-217. [2] Čuš, F., J. Kopač (1998) Inclusion of geometrical shape of cutter into optimization of milling process. Metall 10/11, London, 602-610. [3] Ma, W., J. Li, Y. Rong (1999) Developement of automated fixture planning system. Int J Adv Manuf Technol 15, London, 171-181. [4] Senthil, A., V. Subramaniam, K.C. Seow (1999) Conceptual design of fixtures using genetic algorithms. Int J Adv Manuf Technol 15, London, 79-84. [5] Li, Y, S.Y Liang (1999) Cutting force analysis in transient state milling processes. Int J Adv Manuf Technol 15, London, 785-790. [6] Winbourne, J.P., CM. Toolsie (1989) Computer-aided tool cost estimating. In Proceedings of Computers in Engineering ASME, 617-621. [7] Milfelner, M., F. Čuš (2000) System for simulation of cutting process. International Scientific Conference on the Occasion of the 50th Anniversary of Founding the Faculty of Mechanical Engineering, Ostrava. [8] Zuperl, U. (2000) Development of systems for computer-aided design of modular fixtures. Proceedings of the 11th International DAAAM, Opatija. Naslov avtorjev: Uroš Župerl Prof. dr. Franci Čuš Univerza v Mariboru Fakulteta za strojništvo Smetanova 17 2000 Maribor uros.zuperl@uni-mb.si franc.cus@uni-mb.si Authors’ address: Uroš Župerl Prof. Dr. Franci Čuš University of Maribor Faculty of Mechanical Eng. Smetanova 17 2000 Maribor, Slovenia uros.zuperl@uni-mb.si franc.cus@uni-mb.si Prejeto: Received: 30.5.2001 Sprejeto: Accepted: 29.3.2002 grin^(afcflM]SCLD VH^tTPsDDIK stran 86 © Strojni{ki vestnik 48(2002)2,87-97 © Journal of Mechanical Engineering 48(2002)2,87-97 ISSN 0039-2480 ISSN 0039-2480 UDK 611.72:612.75:007.52:004.94 UDC 611.72:612.75:007.52:004.94 Predhodna objava (1.03) Preliminary paper (1.03) Meritev gibanja kolena z industrijskim robotom - avtomatska kompenzacija gravitacije prijemala Measuring knee movement using an industrial robot – gravity compensation for the automatic gripper Damir Omr~en - Bojan Nemec Poškodbe kolenskih vezi pri športnikih so dokaj pogoste. Za uspešno operacijo vezi je potrebno čim natančnejše poznavanje kolena ter kolenskih vezi. V tem prispevku so opisani postopki določitve geometrijskega modela gibanja kolena. Merjenje gibanja kolena smo izvedli z industrijskim robotom. Uporabili smo robota RIKO 106 s šestimi prostostnimi stopnjami, ki je voden s silo. Površine sklepov smo posneli z uporabo koordinatnega merilnika. Na podlagi meritev smo izdelali geometrijski model na osebnem računalniku s programskim paketom Matlab. Med merjenjem robot upogiba koleno v določeni smeri, ne sme pa vplivati na naravno gibanje kolena. Zato moramo minimizirati sile in navore, ki jih ustvari robot v kolenskem sklepu. V ta namen je treba kompenzirati tudi vplive gravitacije prijemala ter vplive merilnih odmikov zaznavala sile/navora. Med meritvijo krmilimo sile v kolenskem sklepu, zaznavalo sil pa je nameščeno v prijemalki robota, zato je treba izmerjene sile prevesti v koordinatni sistem kolenskega sklepa. V prispevku je natančneje opisan avtomatski postopek za kompenzacijo teže prijemala in določitve merilnih odmikov. Opisana pa sta tudi postopka za preslikavo sil/navorov v kolenski sklep in postopek za določitev vrha orodja, kar potrebujemo pri meritvi. © 2002 Strojniški vestnik. Vse pravice pridržane. (Ključne besede: modeli kolena, prijemala, kompenzacija gravitacije, roboti medicinski, določitev vrha orodja) Injuries to the knee ligaments are very common among althetes. Therefore, a thorough understanding of the knee and the knee ligaments is necessary for successful surgical operations on the ligaments. This paper describes procedures for determining of the geometrical model of the knees movement. The movement of the knee was measured with a RIKO 106, force-controlled, six-degrees-offreedom industrial robot. The surface of the knee joint was scanned with a coordinate-measuring machine and a geometrical model of the knee was developed on a PC. For the modelling we used a computer program called Matlab. The robot should only bend the knee in a specified direction, and should not affect the natural movement of the knee. Therefore, we had to minimize the forces and torques in the knee joint that are caused by the robot. In order to do this, we had to compensate for the influence of gravity on the gripper and the sensor offsets. During the measurement we had to control the forces/torques in the knee joint. As the force/ torque sensor was attached on the robot tip the measured forces/torque had to be mapped to the knee joint. This paper more exactly addresses the automatic procedure for the gripper-weight compensation and the offset determination. It also explains the algorithm of the transformation of the forces/torques to the knee-joint coordinate system and the automatic determination of the tool’s centre point. © 2002 Journal of Mechanical Engineering. All rights reserved. (Keywords: knee models, grippers, gravity compensation, medical robots, tool centre point) 0 UVOD Poškodbe kolenskih vezi pri športnikih so dokaj pogoste. Devetdeset odstotkov vseh poškodb kolenskih vezi so poškodbe sprednje križne vezi (SKV) ter medialne kolateralne vezi (MKV) [1]. Poškodbe MKV se največkrat zacelijo same brez kirurškega 0 INTRODUCTION Injuries to the knee ligaments are very common among athletes. Of all knee ligament injuries, 90 % are injuries to the anterior cruciate ligament (ACL) and the medial collateral ligament (MCL) [1]. The MCL heals without intervention; the ACL requires a grin^otJJiMiscsni stran 87 glTMDDC Omr~en D. - Nemec B.: Meritev gibanja kolena - Measuring Knee Movement posega [2]. Pri poškodbah SKV pa je potreben kirurški poseg [2]. Eden glavnih problemov pri vsaditvi nadomestka SKV je, kako določiti primerno mesto za vsadek. Tega ne moremo vsaditi na mestu, kjer je poškodovani SKV. Primerna so tista mesta, kjer ostane dolžina nadomestka SKV med gibanjem enaka. Zato je za uspešno operacijo vezi potrebno čim natančnejše poznavanje kolena ter kolenskih vezi. V zadnjem času poteka vrsta raziskav na to temo ([2] do [6]). Mnogo raziskav poteka v smeri merjenja sil ter raztezkov v vezeh [4], kar je pripeljalo do sklepa, da se vezi kljub dokaj velikim silam le malo raztegujejo. Nekateri raziskovalci so se usmerili v raziskovanje kinematike kolena ([5] in [6]). S tem želijo ugotoviti natančno gibanje površin sklepov v kolenu. Poznavanje natančnega gibanja bo pripomoglo k boljšemu razumevanju kolena in posledično k uspešnim operacijam kolenskih vezi. Ta prispevek opisuje uporabo industrijskega robota za testiranje in vitro kadaverskega kolena. Uporabljen je robot RIKO 106 s šestimi prostostnimi stopnjami, ki je voden s silo. Z robotom smo zelo natančno izmerili kinematiko neobremenjenega kolena. Na vrhu robota je nameščeno univerzalno zaznavalo sil in navorov, ki je zmožno merjenja treh sil in treh navorov. Za zagotavljanje želenih sil in navorov, ki jih ustvari robot v kolenskem sklepu, smo uporabili hibridno krmilno shemo. Sile, ki jih merimo z zaznavalom, moramo preslikati v kolenski sklep. Kompenzirati je treba tudi gravitacijo prijemala, saj teža prijemala povzroči, da sile v sklepu niso enake nič. Tudi merilni odmiki zaznavala sile vnašajo pogreške v meritev in jih moramo prav tako kompenzirati. Preslikava sil/navorov v kolenski sklep ter kompenzacija mase ter odmikov omogočajo večje sklenjenozančno ojačanje. S tem zmanjšamo pogreške pri sledenju želene sile/navora. Površine sklepov smo posneli z uporabo koordinatnega merilnika MicroScribe. Dobljene podatke smo obdelali s programskim paketom Matlab, v katerem je bila izpeljana tudi grafična predstavitev gibanja kolena. 1 METODE Namen postopka je ugotavljanje natančnega gibanja neobremenjenega kolena (slika 1). Gibanje je sestavljeno iz vrtenja ter premika in ga ne moremo zapisati s preprostimi matematičnimi funkcijami. Zato predlagamo eksperimentalno določitev gibanja z uporabo industrijskega robota s šestimi prostostnimi stopnjami, ki pa mora biti voden s silo. Golenica je trdno pritrjena na mizo, stegnenica pa je preko zaznavala sile pritrjena na vrh robota (sl. 2). Naloga robota je, da upogiba koleno od izravnanega položaja (0 °) do 110 °, pri tem morajo biti sile in navori v sklepu enaki nič. Za vodenje robota s silo je uporabljena klasična hibridna krmilna shema sila-lega. Mig™« | stran 88 reconstruction. One of the main problems when implanting the ACL graft is how to determine a convenient place for the graft. Unfortunately, we cannot implant the graft in a place where an injured ACL is located. Convenient places are those where the length of the ACL graft remains constant during motion. Therefore, a through understanding of the knee and the knee ligaments is necessary for a successful surgical operation. A lot of research has been done on this topic recently ([2] to [6]). An analysis of the forces and the extensions in ligaments [4] led to the conclusion that ligaments extend just a little under high forces. Some researchers studied knee kinematics ([5] and [6]) and they determined the exact motion of the knee-joint surfaces. However, knowing the exact motion of the knee would contribute to a better understanding of the human knee and would result in more effective surgical operations on knee ligaments. This paper describes the use of an industrial robot for in vitro tests on a cadaveric knee. We used a RIKO 106, force-controlled, six-degrees-of-freedom industrial robot to measure the exact kinematics of an unloaded knee. The robot tip was equipped with a universal force/torque sensor, which is capable of measuring three forces and three torques. To keep the desired forces and torques in the knee joint caused by the robot a hybrid force-position control method was used. The measured forces/torques had to be transformed to the knee-joint coordinate system. Additionally, we had to compensate for gravitational forces due to the gripper and the sensor offsets. The transformation of force/torque to the knee-joint coordinate system and the compensation of the gripper weight and offset allowed us to achieve a bigger closed-loop gain. This reduced the force/torque tracking errors. We used the Matlab computer program for visualization, and the surface of the knee joint was scanned with the coordinate-measuring machine. 1 METHODS Our aim was to determine the exact motion of the unloaded knee (Figure 1). This motion is a combination of rotation and translation and it cannot be expressed with simple analytical functions. We decided to determine the motion experimentally using a force-controlled robot with six degrees of freedom. As shown in Figure 2 the tibia was attached to the table and the femur was attached to the force sensor on the robot tip. The robot had to move the knee from the straightened position (0 °) to 110 ° and during the motion the forces and torques had to be kept as small as possible. To achieve this we applied the well-known hybrid force-position control method. Omr~en D. - Nemec B.: Meritev gibanja kolena - Measuring Knee Movement Sl. 1. Človeško koleno Sl. 2. Vpetje golenice in stegnenice Fig. 1. Human knee Fig. 2. Attachment of the tibia and femur Ker teža prijemala vpliva na meritev sile, jo As the gripper’s weight influences the force moramo kompenzirati pred začetkom meritve. Nato measurement it had to be compensated before the tibia golenico trdno pritrdimo na mizo. Na golenico in was attached to the table. Next we had to determine the stegnenico pritrdimo po tri kalibracijske točke. Z positions of the three calibration points that were fixed orodjem na vrhu robota izmerimo tri kalibracijske on the tibia. With the tool on the robot tip we measured točke na golenici in jih zapišemo v koordinatnem the coordinates of the three calibration points in the robot- sistemu baze robota. Prijemalo robota pripeljemo do base coordinate system. Next, we had to move the gripper stegnenice in jo pritrdimo. Meritev poteka tako, da to the femur and attach thefemur to the gripper. Then, the koleno najprej upognemo za 7 °, nato minimiziramo robot bent the knee by 7 °. The next step was to minimize sile in navore v kolenu ter izmerimo lego kolena. the forces and torques in the knee joint and to measure Postopek ponavljamo, dokler ne upognemo kolena the position of the knee. The bending phase and force- od povsem izravnane do skrajne želene lege. Pri minimizing phase were repeated until the knee reached upogibu kolena vodimo robota hitrostno v smeri the desired position (110 °). In the bending phase the upogibanja kolena, v preostalih petih smereh pa robot was velocity controlled in the direction of the knee poteka krmiljenje sile. Pri minimizaciji sil in navorov rotation and in all other directions it was force controlled. preklopimo na krmiljenje sil in navorov v vseh In the force-minimization phase the forces were controlled smereh. in all directions. Po končani meritvi gibanja izmerimo še After the measurement cycle we had to measure pozicije kalibracijskih točk na stegnenici glede na the coordinates of the calibration points on the femur in vrh robota. Iz teh podatkov in iz konfiguracij robota the robot-tip coordinate system. These data and the izračunamo položaje kalibracijskih točk med robot configurations data were used to calculate the gibanjem kolena. Nato se lotimo meritve geometrijske positions of the calibration points during the knee oblike kolena. motion. After the motion measurement the knee Pred merjenjem geometrijske oblike površine geometry had to be determined. sklepa najprej odstranimo vse tkivo ter vezi. Nato Before scanning the surface of the knee we izmerimo koordinate mreže točk na površinah sklepov had to remove all tissues and ligaments. The coordi- ter na prijemališčih SKV. Te točke so izražene v nate-measuring machine was used to measure the koordinatnem sistemu koordinatnega merilnika in jih coordinates of points on the knee’s surface and on moramo prenesti v koordinatni sistem baze robota, za the ACL grip. These coordinates were then trans- to potrebujemo kalibracijske točke na golenici in formed to the robot-base coordinate system. For this stegnenici. S koordinatnim merilnikom izmerimo tudi reason calibration points on the femur and tibia were koordinate kalibracijskih točk. Ta izračun poteka v needed, and these were also measured with the coor- programskem paketu Matlab, v katerem je realizirana dinate-measuring machine. The calculations and visu- tudi grafična predstavitev gibanja. alisation were made with the Matlab program. 1.1 Preslikava sil/navorov iz koordinatnega sistema 1.1 Transformation of force/torque from the force (k.s.) zaznavala sile v k.s. kolenskega sklepa sensor coordinate system (c.s.) to the knee joint c.s. S na sliki 3 prikazuje koordinatni sistem Figure 3 shows the coordinate system senzorja sil, K pa koordinatni sistem, ki je v connected to the force sensor S and the coordinate središču kolenskega sklepa. Koordinatna sistema system connected to the centre of the knee joint K. sta med seboj vzporedna in premaknjena za vektor These two coordinate systems should be parallel to p. Zveza med silami F in navori M v obeh each other and translated. The vector p is the translation 02-2 stran 89 lMBuTMHK Omr~en D. - Nemec B.: Meritev gibanja kolena - Measuring Knee Movement Sl. 3. Preslikava sil/navorov iz koordinatnega sistema zaznavala S v koordinatni sistem kolena K Fig. 3. Force/torque transformation from the force sensor coordinate system S to the knee c. s. K koordinatnih sistemih S in K je: between their origins. The relation between the forces F and the torques M in the coordinate system S and K are described by the following equations: (1) kjer se (.)S in (.)K nanašata na izbrani koordinatni sistem, p^ pa je operator vektorskega zmnožka, povezan z vekorjem p: P = Oznake (.), (.) ter (.) se v članku nanašajo na x, y ter z (.) komponente vektorja. 1.2 Avtomatska kompenzacija gravitacije Pri neobremenjenem zaznavalu želimo, da so izmerjene sile/navori enaki nič, kar pa zaradi mase prijemala in zaznavala ni res. Te sile/navore moramo kompenzirati. Masa ter težišče prijemala in zaznavala sta znana le redko. Izračun pa je mnogokrat prezapleten in ni primeren za praktične uporabe. Zato predlagamo določitev kompenzacijskih sil/navorov z meritvami. S predlaganim postopkom določimo maso in težišče prijemala in zaznavala ter tudi merilna odstopanja zaznavala. 0 -Pz Py 1 Pz 0 /'a [~Py Px 0 J (2). where (.)S and (.)K denote the particular coordinate system and p^ is the cross-product operator associated with vector p: (3). The notation (.)x, (.)y and (.)z stands for the x, y and z components of vector (.), respectively. 1.2 Automatic compensation for gravity In the case of an unloaded sensor the measured forces/torques should be equal to zero; however, due to the weight of the gripper and the sensor this is not the case and these forces/torques have to be compensated. The centre of gravity and the weight of the gripper and the sensor are rarely known. We could calculate the weight and the centre of gravity, which requires an exact knowledge of the geometry and materials; however, this calculation is often too complex and is not suitable for practical application. We propose to determine the compensation forces with measurements. Using this procedure we can determine the centre of gravity and the weight of the gripper and the sensor and also the sensor offset. Sl. 4. Vpliv mase prijemala m na meritev sile/navora v zaznavalu S Fig. 4. Influence of the gripper weight on the measurement of force/torque by sensor S M1©T[M]C€ ^jr[^(g)J[R!]OŽi[€D | stran 90 Omr~en D. - Nemec B.: Meritev gibanja kolena - Measuring Knee Movement R naj pomeni robotovo rotacijsko matriko med k.s. S in k.s. okolja W, r pa je vektor med njunima izhodiščema. Kompenzacijsko silo/navor izračunamo: Let R be the robot rotational matrix between the sensor c.s. S and the world c.s. W and r is the vector between their origins. Then the compensation force/torque can be calculated as: Fcamp = RFg =[n o a ] Fg = [n o a] 0 0 LF& aFg2 ***-c!i F„ff (11). Sistem enačb je predoločen. Z uporabo nevidezne obratne vrednosti dobimo rešitev po metodi najmanjših kvadratov: r . n /r ... Ur This system is overdeterminated. Using the pseudo-inverse gives the least-square solution: ¦off oP) ap) rffM (12). Določiti moramo še težišče prijemala ter merilne odmike navora. Pri tem potrebujemo tri meritve. Z uporabo (5) in (10) dobimo: However, we still have to calculate the centre of gravity and the torque offset. For these calculations we need three measurements. Using (5) and (10) we get: (13) grin^otJJiMiscsni stran 91 glTMDDC Omr~en D. - Nemec B.: Meritev gibanja kolena - Measuring Knee Movement in zopet z uporabo psevdoinverza: and using the pseudo-inverse yields: M, off \ p(1) ~-F com p -Tcomp n F r rrUJ rcomp — Fcomp T? (3) Obračanje matrike velikosti 6x6 je problem pri uporabi preprostih računalnikov, ki se uporabljajo za vodenje industrijskih robotov. Zato predlagamo postopek za zmanjšanje velikosti invertirane matrike. Zapišimo (13) v razširjeni obliki: p(0 Fcomp p(2) — rctmtp F(3) fftMlfp t r M(1) ' Mt5> (14). Inversion of a 6x6 matrix can be a problem when a simple computer is used for the robot control. We propose the following procedure for reducing the necessary arithmetic operations. Let us rewrite (13) in its expanded form: 1VJ-metis M^lp + Moff M^l,p + Moff _M^mp + MofL 0 F(i) „(i) i o o" p(l) rcompv o F(i) — "comp^. a 0 1 0 0 0 1 0 _fP) rcompz F(2) 1 0 0 p0) rCOHtp2 p(2) rcompv o _FP) ra>mpz 0 0 1 0 0 0 1 0 p(3) F(3) rcimpy, 1 0 1) p(3) p(3) 0 rcompz p(3) ra>mp:! 0 0 1 0 0 0 ] M, M, -»i r, <-_ »//,. M, 0//J (15). V tej enačbi se pojavlja šest neznank in devet enačb, od katerih je le 6 linearno neodvisnih. Če odstranimo 3., 5. in 7. vrstico, se izkaže, da dobimo 6 linearno neodvisnih vrstic (enačb): Za izračun r in M ff je treba obrniti matriko velikosti 6x6. Postopek poenostavimo, če zamenjamo 3. in 4. vrstico iz enačbe (16) in dobimo: Note that there are six unknowns and nine equations, and only six of them are linearly independent. By removing the 3rd, 5th and the 7th rows (equation) we obtain six linearly independent equations: (16). For the calculation of r and Moff we have to invert a 6x6 matrix. The calculation can be simplified if we swap the 3rd and the 4th rows in (16) as follows: (17). Zapišimo zgornjo enačbo drugače: Rewriting the upper equation yields: Mh Izrazimo težišče in merilne odmike navorov: (18). The centre of gravity and the torque offset are calculated as: (19) (20). grin^(afcflM]SCLD VBgfFMK stran 92 Omr~en D. - Nemec B.: Meritev gibanja kolena - Measuring Knee Movement Za izračun r je treba obrniti le matriko (F F ), ki je velikosti 3x3, kar pa ne dela posebnih težav. Sedaj, ko poznamo vpliv gravitacije ter merilnih odmikov, moramo ta dva prispevka od meritve odšteti, da dobimo pravo vrednost: For the calculation of r we only have to invert the matrix (Fb – Fa) which is of size 3x3. For the calculation of the force/torque without the influence of gravity on the gripper and the offset we have to subtract the gravity and offset influence from measured force/torque: (21) 1.3 Avtomatska določitev vrha orodja Pri mnogih robotskih uporabah je treba določiti lego vrha orodja glede na vrh robota, prav tako potrebujemo znan vrh orodja za določitev koordinat kalibracijskih točk na golenici. To lego vrha lahko izračunamo ali pa izmerimo. Ker je meritev lahko dokaj nenatančna, predlagamo lažji in hitrejši postopek za avtomatsko določitev lege vrha. (22). 1.3 Automatic determination of the tool’s centre point In many robot applications it is necessity to determine the position of the robot tool’s centre point with respect to the robot tip. For the determination of the calibration points on the tibia we also need to know the position of the tool’s centre point. This position can be measured or calculated. We propose a simple automatic procedure for the determination of the tool’s centre point. Sl. 5. Prenos baznega koordinatnega izhodišča v določeno točko v prostoru Fig. 5. Transfer from the robot base to a fixed point in space Z vrhom orodja se trikrat dotaknemo poljubne točke, vendar vsakokrat z drugo usmerjenostjo vrha robota. Pozicijo določene točke lahko izrazimo na dva načina (sl. 5): With the tool tip we have to touch the same point in space every time using different orientations of the robot tip. The position of the point can be written in two ways: TO + RRTR = T (23). S T so označene translacije med posameznimi koordinatnimi sistemi, z R pa rotacije. (.)R se nanaša na robotove transformacije, (.)O na transformacije orodja, (.) pa na transformacije točke. Pozicija neke točke glede na bazo je sestavljena iz translacije robota T ter translacije orodja T , ki je zarotirana za rotacijsko matriko robota R. Pozicijo te iste točke pa lahko opišemo tudi s translacijo TT od baze do točke. V tej enačbi sta znani translacija in rotacija robota, transformacije točke ne poznamo, iščemo pa translacijo orodja T. Za določitev neznank potrebujemo tri meritve (.)(.) tri konfiguracije robota, tako da vrh orodja kaže v enako točko. Zapišimo sistem s tremi konfiguracijami v matrični obliki: T denotes a translation between two coordinate systems, R denotes rotations between them and (.)R relates to the robot transformation, (.)O relates to the tool transformation and (.)T relates to the point transformation. The position of the point with respect to the base is calcualted from the robot translation TR and the tool translation TO, which is rotated with a robot rotational matrix RR. The position of the same point can be described as translation TT from the base to the point. The robot translations and rotations are known, however, we have to obtain the tool translation TO. Equation (23) has three equations and six unknowns. In order to solve eq. (23) we need three robot configurations (.)(.) when the tool centre point points to the same point in space with different orientations. Combining (23) for three measurements we obtain: stran 93 02-2^DCCl Omr~en D. - Nemec B.: Meritev gibanja kolena - Measuring Knee Movement Sistem je mogoče rešiti z uporabo navidezne obratne vrednosti: T0 = ( '< -i" RR] -' .45) -i. i K -i] <2) -¦ .43) -i ) -L K) -I] R -I .4" -i T r tt(i)i T:'2) Zaradi zmanjšanja velikosti obratne matrike uporabimo postopek, ki je bil opisan v prejšnjem poglavju: 4n 4" 4n -10 0 0-10 0 0-1 'T0 = - T(i) % ?f "v T<3) 'R, «<3) 43) «L> «l3' 0<3> *i3) -10 0 0-10 0 0-1 (24). This system could be solved using the pseudo-inverse: (25). However, we propose the same procedure as in the previous section, which results in an easier matrix inversion: (26). Zaradi preprostejšega obračanja matrike, zapišemo sistem enačb: in dobimo: Z uporabo te preproste enačbe določimo pozicijo vrha orodja glede na vrh robota. Vrh orodja je treba trikrat čim natančneje pripeljati v določeno točko v prostoru. Natančnost določitve TO je odvisna od natančnosti robota in natančnosti dotika točke. Brez težav pa lahko dosežemo natančnost pod 0,5 mm. 1.4 Obdelava podatkov ter grafična predstavitev Točke, ki smo jih izmerili na koordinatnem merilniku, moramo prenesti v koordinatni sistem robota. Za ta namen uporabimo po tri kalibracijske točke (A, B, C) na golenici in stegnenici. Iz teh treh točk najprej določimo rotacijsko matriko ravnine, ki poteka skozi te tri točke (slika 6). Rewrite (26) into: and finally we obtain (27) (28). The position of the tool tip with respect to the robot tip is determined using this simple equation. We only need to move the robot tool to same point in space. The accuracy of TO depends on the robot accuracy and the accuracy of the touch. The accuracy can be better that 0.5 mm. 1.4 Data calculating and graphical presentation The points measured with the coordinate-measuring machine had to be transformed to the robot-base coordinate system. To do this three calibration points (A, B, C) on the tibia and femur are needed. Using these three points we define the rotational matrix of the plane, which contains these three points. Sl. 6. Določanje rotacijske matrike ravnine Fig. 6. Defining of the plain rotational matrix grin^(afcflM]SCLD VBgfFMK stran 94 Omr~en D. - Nemec B.: Meritev gibanja kolena - Measuring Knee Movement Med točkama B in A ter B in C določimo enotina vektorja j ter m. Normirani vektorski zmnožek med njima da normalo na ravnino (l). Za določitev vrtilne matrike pa potrebujemo še en pravokotni vektor k: Between points B and A and between B and C we define the unit vectors j and m, respectively. The normalized vector product between j and m gives us the plane normal l. To define the rotational matrix we need one more prependicular vector k: (29) (30) Kot translacijo vzamemo pozicijo točke B. Tako dobimo po dve rotacijski matriki (R) ter dve translaciji (T) za golenico in stegnenico, eno za meritev z robotom (.)R in drugo za meritev s koordinatnim merilnikom (.) . Vse točke na površinah sklepov (P), ki so izmerjene s koordinatnim merilnikom, prenesemo v koordinatni sistem robota: Podatke smo preračunali s programskim paketom Matlab, v katerem je bil izveden tudi izris (slika 9). 2 REZULTATI Meritve smo izvajali na prašičevem kolenu, ki je po zgradbi najbolj podobno človeškemu (sl. 7). Izvedli smo vrsto 10 meritev z zelo dobro ponovljivostjo. Sile ter navori v sklepu so bili znotraj 0,5 N oziroma 0,05 N m. Vzporedno z merivijo z robotom smo izvajali tudi testno merjenje s koordinatnim merilnikom MicroScribe, katerega natančnost je 0,01 mm. Meritev z robotom se je izkazala za dovolj natančno. S koordinatnim merilnikom smo pomerili tudi površine sklepov (sl. 8), in sicer na vsaki površini sklepa po 20 do 30 točk. Grafična predstavitev je potekala v programskem paketu Matlab (sl. 9), pri katerem smo tudi računali dolžino l srednjega vlakna SKV (sl. 10). 3 SKLEP V prispevku smo obravnavali postopek za natančno določitev gibanja kolenskega sklepa. Pri tem smo uporabili industrijskega robota, ki je voden s silo. Zaradi tega je bilo treba kompenzirati vpliv (31) (32) (33). The translation T is the position of the point B. So we have two rotational matrices (R) and two translations (T) for tibia and femur, one for robot (.)R measurement and one for the coordinate-measuring (.) machine measurement. All the points on the knee surfaces (P) that were measured with coordinate measuring machine had to be transformed to the robot-base coordinate system using: -Tcm)) + TR (34). The calculations and graphical presentation were made with the Matlab program (Figure 9). 2 RESULTS The measurements were performed on a porcine knee, the structure of which is similar to a human knee (Figure 7). We performed a series of 10 measurements. Forces and torques in the knee joint were within 0.5 N and 0.05 Nm, respectively. For control purposes we also measured the positions of control points with the MicroScribe coordinate-measuring machine. The accuracy of the MicroScribe machine is better than 0.01 mm. We checked the measurement accuracy of the robot with MicroScribe. The measuring accuracy of the robot was better than 0,5 mm. Using MicroScribe we also scanned the surface of the knee joint (Figure 8). On every surface we measured 20 to 30 points. A graphical presentation of the knee movement was made with the Matlab program package (Figure 9), the length l of the ACL during the motion was also computed (Figure 10). 3 CONCLUSIONS This paper describes procedures for a determination of the geometrical model of knee movement. The movement of the knee was measured with a force-controlled industrial robot. To do this we had to com- gfin^OtJJIMISCSD stran 95 Omr~en D. - Nemec B.: Meritev gibanja kolena - Measuring Knee Movement Sl. 7. Postopek meritve Sl. 8. Snemanje površine sklepa Fig. 7. Measurement procedure Fig. 8. Scanning of knee surfaces Sl. 9. Grafični prikaz površin kolena v Matlabu Sl. 10. Prikaz dolžine SKV med gibanjem kolena Fig. 9. Graphical presentation of the knee surfaces Fig. 10. Length of the ACL during knee movement in Matlab gravitacije prijemala na zaznavalo sile/navora ter pensate the influence of gravity on the gripper and merilne odmike zaznavala. the sensor offsets. Avtomatska kompenzacija gravitacije The automatic gravity and offset compensa- prijemala ter merilnih odmikov se je izkazala kot zelo tion results in a very useful robotic tool. The proce- uporabno orodje v robotiki. Postopek je preprost, dure is very simple and sufficiently accurate, even hiter in dovolj natančen, tudi kadar imamo opravka when we are dealing with small forces and torques. In z majhnimi silami/navori. V našem primeru our case the compensation ensures small forces/tor- kompenzacija zagotavlja majhne sile/navore v ques in the knee joint and therefore a better control sklepu, s tem pa dosežemo boljšo stabilnost stability. An additional contribution to the stability is vodenja. Dodaten prispevek k stabilnosti je prinesla the transformation of the forces/torques from the force/ preslikava sil/navorov v sklep. torque sensor to the knee-joint coordinate system. Razvili smo celovit sistem za meritev in prikaz We developed a system for measuring and gibanja kolen. Ta sistem bi bilo mogoče uporabiti v graphically presentating of human knees. This system can zdravstvene namene za prikaz gibanja kolena. S be used for medical purposes, for the animation of the knee programom lahko izračunamo raztezke vezi, ki so movement. The most important thing in ligaments pritrjene nekje na kosti, to pa je pri zamenjavi vezi implantation is calculating the ligaments extension. This ključnega pomena. calculation can be performed with the developed system. 4 LITERATURA 4 REFERENCES [1] Miyasaka et al. (1991) The incidence of knee ligament injuries in the general population, Am. J. Knee Surg 4, 3-8, 1991. ^^™DCš| stran 96 Omr~en D. - Nemec B.: Meritev gibanja kolena - Measuring Knee Movement [2] Savio L-Y. Woo et al. (1997) Biomechanics of knee ligament healing, repair and reconstruction, Journal of Biomechanics, Vol. 30, No. 5, 431 - 439, 1997. [3] Hiromichi Fujie et al. (1993) The use of robotics technology to study human joint kinematics: A new methodology, Journal of Biomechanical Engineering, Vol. 115, 211- 217, 1993. [4] Hiromichi Fujie et al. (1995) The use of a universal force-moment sensor to determine in-situ forces in ligaments: A new methodology, Journal of Biomechanical Engineering, Vol 117, 1-7, 1995. [5] Bojan Nemec et al. (2000) The use of robotics technology for kinematics test of synovial joints in medicine. V CURK, Boris (ur.), HARNIK, Jože (ur.). 9th International Workshop on Robotics in Alpe-Adria-Danube Region RAAD 2000, Maribor, Slovenia, June 1-3, 2000. [6] Bojan Nemec et al. (1998) The use of force controlled robot for kinematic test in medicine. V DOBROVODSKY, Karol (ur.). 7th International Workshop on Robotics in Alpe-Adria-Danube Region, Smolenice Castle, Slovakia, June 26-28, 1998. Proceedings : RAAD 1998. Bratislava: Institute of Control Theory and Robotics: Slovak Academy of Science, 1998, str. 435-440. [7] Tadej Bajd, Alojz Kralj (1997) Robotika, Fakulteta za elektrotehniko, Ljubljana 1997. Avtorjev naslov: Damir Omrčen dr. Bojan Nemec Institut Jožef Stefan Jamova 39 1000 Ljubljana damir.omrcen@ijs.si Author’s Address:Damir Omrčen Dr. Bojan Nemec “Jožef Stefan” Institute Jamova 39 1000 Ljubljana, Slovenia damir.omrcen@ijs.si Prejeto: Received: 6.8.2001 Sprejeto: Accepted: 29.3.2002 grin^otJJiMiscsni stran 97 glTMDDC © Strojni{ki vestnik 48(2002)2,98-111 © Journal of Mechanical Engineering 48(2002)2,98-111 ISSN 0039-2480 ISSN 0039-2480 UDK 658.53:65.015.3:65.015 UDC 658.53:65.015.3:65.015 Strokovni ~lanek (1.04) Specilaity paper (1.04) Realni preto~ni ~asi operacij in uspe{nost sistema NKP Realistic Lead Times of Operations and Efficiency of the PPC System Marko Starbek - Janez Grum - Janez Ku{ar Uspešnost dela na trgu razpoložljivih komercialnih sistemov NKP (sistemov za načrtovanje in krmiljenje proizvodnje) [1] in [2] je v veliki meri odvisna od realnosti izvedbe pretočnega terminiranja operacij, ki sloni na ocenjenih oziroma ugotovljenih pretočnih časih delovnih mest. V prispevku predlagamo postopek za ugotavljanje srednjih pretočnih časov delovnih mest. V časovnem koraku P ugotovljene vrednosti srednjih pretočnih časov delovnih mest pokažejo dejansko sliko pretoka naročil preko delovnih mest in pomenijo osnovni podatek za izvedbo realnega pretočnega terminiranja operacij v naslednjem, to je časovnem intervalu (P+1). S predlaganim postopkom ugotavljanja pretočnih časov preidemo od statičnega na dinamično pretočno terminiranje operacij naročil. © 2002 Strojniški vestnik. Vse pravice pridržane. (Ključne besede: sistemi PPS, časi pretočni, modeli lijaka, diagrami pretoka) The practical efficiency of commercially available PPC (production planning and control) systems [1] and [2] depends to a great extent on the implementation of operations’ lead times determination which is based on estimated or calculated lead times of workplaces. This paper presents a method for calculating the average lead times of workplaces. The average lead times of workplaces, found in the interval P, gives a realistic picture of the order flow through workplaces, and these data are the basis for the implementation of the realistic operations’ lead times determination in the next (P+1) interval. With the proposed method for calculating lead times a transition from static to dynamic orders’ lead times determination is accomplished. © 2002 Journal of Mechanical Engineering. All rights reserved. (Keywords: PPC-systems, lead times, funnel models, flow diagrams) 0 UVOD Namen sistema NKP je izvedba nalog načrtovanja in krmiljenja proizvodnje [3]. Načrtovanje proizvodnje mora poskrbeti za načrtovanje poteka proizvodnje v prihodnosti in vključuje naloge: načrtovanje primarnih potreb, načrtovanje materialnih potreb, načrtovanja pretočnih časov operacij in izravnavo kapacitet. Krmiljenje proizvodnje pa mora poskrbeti za izvedbo načrtovane proizvodnje, in to kljub neizogibnim spremembam količin, dogovorjenih terminov, pomanjkanju osebja, delovnih sredstev, in vključuje naloge: izpuščanje naročil, fino terminiranje ter monitoring in kontroling proizvodnje. Realno načrtovanje in krmiljenje proizvodnje terjata stalno spremljanje in izvedbo potrebne korekcije toka materiala in informacij, kar je mogoče doseči edino s sistemom NKP. Slika 1 prikazuje mesto sistema NKP v proizvodnem podjetju. 0 INTRODUCTION A production planning and control (PPC) system has to carry out production planning and control tasks [3]. Production planning has to plan the course of future production. It consists of the following tasks: primary needs planning, material needs planning, operations’ lead times determination, and levelling of resources. Production control has to provide for the realisation of the planned production, regardless of the inevitable changes to quantities, agreed terms, loss of staff and production means. It consists of the following tasks: releasing orders, fine termination, and production monitoring and control. Realistic production planning and control requires continuous monitoring and implementation of the required corrective measures in material and information flow, which can only be achieved by a PPC system. The place of the PPC system in a production company is shown in Figure 1. grin^(afcflM]SCLD VBgfFMK stran 98 Starbek M. - Grum J. - Ku{ar J.: Realni preto~ni ~asi - Realistic Lead Times tok informacij > information flow ^^ tok materiala ^* material flow Sl. 1. Mesto sistema NKP v proizvodnem podjetju Fig.1. The position of the PPC system in a production company Informacija, ki jo sistem NKP dobi od prodaje, je informacija o proizvodnem načrtu izdelkov, ki jih bo treba izdelati v opazovanem časovnem koraku (naročila kupcev, načrtovana zaloga). Sistem NKP proizvodni načrt deli na izdelovalna naročila (lastna proizvodnja) in nabavna naročila (nabava na trgu surovin in polizdelkov) in pri tem upošteva skladiščno stanje tako kupljenih kakor gotovih izdelkov. Po vstopu surovin in polizdelkov v vhodno skladišče teče materialni tok preko proizvodnje in skladišča gotovih izdelkov h kupcu. Stalno nadzorovanje naročil in zmogljivosti je namenjeno za oskrbovanje sistema NKP z nujnimi povratnimi informacijami. Za izdelovalna naročila opazovanega časovnega koraka P mora sistem NKP izvesti pretočno terminiranje operacij in za izvedbo te naloge potrebuje podatke o srednjih pretočnih časih delovnih mest v predhodnem časovnem koraku (P-1). Analiza rezultatov pretočnega terminiranja operacij, izvedena v več majhnih in srednje velikih podjetjih, je pokazala, da podjetja pri pretočnem terminiranju delujejo z grobimi izkustvenimi vrednostmi srednjih pretočnih časov delovnih mest in zato ne preseneča dejstvo, da prihaja do velikih odstopanj med načrtovanimi in doseženimi termini. The sales department sends the PPC system the information about the production plan of products that will have to be manufactured in the treated interval (orders of customers, planned stock). The PPC system divides the production plan into production orders (company’s own production) and purchase orders (purchases on the market of raw materials and semi-products), taking into account the bought and manufactured products in stock. After raw materials and semi-products have entered the warehouse of bought parts, the material flows from production through the warehouse of final products to the customer. The continuous control of orders and resources provides feedback information to the PPC system. The PPC system has to perform the operations’ lead times determination for production orders of the treated interval P. In order to accomplish this task it needs the data on average lead times in the previous (P-1) interval. An analysis of the results of the flow termination of operations, which was carried out in several small and medium-sized companies, has shown that when dealing with lead times determination the companies used approximate experience-based values of the average lead times of workplaces and therefore it was not surprising that large discrepancies existed between the planned and actual terms. gfin^OtJJlMISCSD stran 99 Starbek M. - Grum J. - Ku{ar J.: Realni preto~ni ~asi - Realistic Lead Times Proizvodni načrt proizv. P1 Production plan of product P1 -izdelovalna naročila -nabavna naročila -production orders -purchase orders Struktura proizvoda P1 Product P1 structure Strukturna kosovnica proizvoda P1 Product P1 bill of material Ravnina: Level: 0. 1. 2. Tehnološki postopki komponent proizvoda P1 Technology procedures of product P1 components Strukt. kos. za P1 Prod. P1 bill of mat. T.Š. No. Naziv Name Ravnina Level Količ. [kos] Qty [pcs] 1 SD1 . 2 2 SD2 . 5 Ocena srednjih pretočnih časovdelovnih mest [Dd] Estimated average lead times of workplaces _______[Wd]_______ ni Tehnološki pki postopek / Technology SD procedure SD1 No. Oper. TP1 TP*2 TP*j . TP*n Izdelovalno naročilo Production order Izdelovalno naročilo Production order Montažno naročilo Assembly order Naročila Orders A SD SD P1 tp1 tp2 TP3 najkasnejši pretočni termini SD1 latest flow terms of SD1 --------# TP5 TP7, najkasnejši pretočni termini SD2 latest flow terms of SD2 najkasnejši pretočni termini P1 latest flow terms of P1 NAČRTOVANI KONČNI TERMIN PROIZVODNEGA NAROČILA PLANNED FINAL TERM OF A PRODUCTION ORDER Čas [Dd] >Time [Wd] <- flow time of a production order M Sl. 2. Primer retrogradnega pretočnega terminiranja operacij proizvodnega naročila Fig. 2. An example of the reverse-flow order based lead time determination Slika 2 prikazuje pregled dokumentov, ki omogočajo pretočno terminiranje operacij ter prikazuje načelo retrogradnega pretočnega terminiranja izdelka P1, ki je sestavljen iz dveh sestavnih delov SD1 in SD2 (proizvodno naročilo sestoji iz dveh izdelovalnih in enega montažnega naročila). Figure 2 presents an overview of the documents that allow for the flow termination of operations and shows the principle of reverse-flow termination of the P1 product, composed of parts SD1 and SD2 (production order consists of two manufacturing orders and one assembly order). VBgfFMK stran 100 . T te pz . 1 O 11 2 O 12 3 O 13 Starbek M. - Grum J. - Ku{ar J.: Realni preto~ni ~asi - Realistic Lead Times Sistem NKP lahko na temelju podatka o načrtovanem končnem terminu proizvodnega naročila in podatkov o proizvodnem načrtu, strukturi in kosovnici proizvoda, tehnoloških postopkih izvedbe operacij na komponentah proizvoda in ocenjenih srednjih pretočnih časih delovnih mest izvede pretočno terminiranje tehnoloških in montažnih operacij. Da bi v prihodnosti uspešno rešili problem velikega odstopanja načrtovanih terminov od doseženih, smo se raziskovalci Laboratorija za proizvodne sisteme odločili, da sestavimo postopek za trajno ugotavljanje doseženih srednjih pretočnih časov delovnih mest oziroma operacij, saj so v več podjetjih izvedene analize pokazale, da se vrednosti srednjih pretočnih časov delovnih mest stalno spreminjajo. Analiza vzrokov nenehnega spreminjanja pretočnih časov je pokazala, da na velikost srednjih pretočnih časov delovnih mest vpliva množica dejavnikov, naj naštejemo le najpomembnejše: - obremenitev delovnih mest in čakalna vrsta naročil, - način podrobnega razporejanja naročil na delovno mesto, - kakovost strojev, naprav in orodij, - način oskrbe z materiali, orodji in informacijami, - usposobljenost in motiviranost delavcev, ki strežejo strojem, - način organiziranosti transportne in vzdrževalne službe. Na delovno mesto Sj prihajajoča naročila pomenijo obremenitev delovnega mesta, čakajoča naročila pomenijo stanje naročil, z delovnega mesta odhajajoča naročila pa učinek delovnega mesta. V določeni meri spremenljiva odprtina lijaka pomeni razpoložljive oziroma izkoriščene zmogljivosti delovnega mesta. Naloga krmiljenja delovnega mesta je, da po eni strani preprečuje preveliko stanje naročil na delovnem mestu in s tem predolge pretočne čase naročil in po drugi strani premajhno stanje, ki bi lahko povzročilo nedelo delovnega mesta zaradi pomanjkanja naročil. 1 UGOTAVLJANJE SREDNJIH PRETOČNIH ČASOV DELOVNIH MEST Pretok množice naročil N1,...,Ni,...,Nm preko delovnega mesta Sj v opazovanem časovnem koraku P si lahko ponazorimo s poznanim modelom lijaka [4], ki ga prikazuje slika 3. 1.1 Zbiranje podatkov o pretoku naročil Da bi prišli do podatkov o pretočnih časih naročil Ni (1 = i = m), ki so v opazovanem časovnem koraku P prešla preko delovnega mesta S, moramo za vsako naročilo zbrati podatke o: - številki naročila, ki je prišlo oziroma odšlo z The PPC system can carry out the flow termination of technology and assembly operations on the basis of the data of the planned final term of production order and data on the product plan, the structure and the bill of materials of the product, the technology procedures for the execution of operations on the product components, and the estimated average lead times of workplaces. In order to solve the problem of the large discrepancies between the planned and actual terms the Production Systems Laboratory researchers decided to develop a method for continuous monitoring of the actual lead times of workplaces (operations), as the analyses performed in several companies revealed that the values of the average lead times of workplaces were changing continuously. The analysis of the continuous changing of lead times has shown that there were many reasons why these values changed, the most important being: - the load of the workplaces, the queue of orders, - the method of fine allocation of orders to a workplace, - the quality of machines, devices and tools, - the method of material, tool and information supply, - the qualifications and motivation of employees working with machines, - the organisation of transport and support services. In a funnel model the orders that come to the workplace S represents the load of the workplace, the queued orders represent the state of the orders, and the orders departing from the workplace represent the workplace efficiency. The hole of the funnel (variable to a certain extent) represents the available and used resources of the workplace. The control system has two tasks: it prevents too high a state of orders in the workplace (and thus excessive lead times for the orders), and it prevents too low a state of orders which could put the workplace out of operation due to an insufficient number of orders. 1 CALCULATION OF AVERAGE LEAD TIMES OF WORKPLACES The flow of orders N1,...,N,...,N through the workplace S in the treated interval P can be represented by the funnel model [4], as shown in Figure 3. 1.1 Acquisition of orders flow data In order to obtain the data on the lead times of orders Ni (1