28 Izvirni znanstveni članek TEHNIKA – zobniška gonila Datum prejema: 7. december 2019 ANALI PAZU 10/2020/1-2: 28-35 www.anali-pazu.si Optimalno določanje stabilnih območij vibracij zobniških gonil Optimal determination of stability domains of vibrations in gear drives Rudolf Pušenjak 1 * 1 Fakulteta za industrijski inženiring Novo mesto/Šegova ulica 112, Novo mesto E-mails: rudolf.pusenjak@fini-unm.si * Rudolf Pušenjak Povzetek: V članku je obravnavana optimizacija stabilnih območij vibracij dvostopenjskega zobniškega gonila s parametričnim vzbujanjem kot posledice periodično spremenljive zobne togosti med ubiranjem. Da bi zagotovili mirnejši tek zobniškega gonila in zmanjšali hrup, želimo v procesu snovanja zobniškega gonila območja nestabilnosti vibracij čimbolj zmanjšati. Na zmanjšanje posameznega območja nestabilnosti v največji meri vplivata kontaktno razmerje posamezne zobniške dvojice in njun fazni premik. Na žalost se ob zmanjšanju enega nestabilnega območja lahko povečajo ostala nestabilna območja. S tem se pojavi potreba po iskanju optimalne rešitve, ki jo lahko zagotovimo z rešitvijo večkriterijskega optimizacijskega problema. V članku je predstavljena večkriterijska optimizacija območij stabilnosti vibracij zobniških gonil z metodo utežne vsote in z metodo ε – omejitev. Rezultate večkriterijske optimizacije območij stabilnosti vibracij v članku analitično pojasnimo. Obe večkriterijski optimizacijski metodi sta v članku prikazani pri dvostopenjskemu zobniškemu gonilu s tremi zobniki, pri katerem sta frekvenci ubiranja enaki. Ključne besede: vibracije zobniških gonil, spremenljiva zobna togost, kontaktno razmerje, fazni premik, optimizacija območij parametrične nestabilnosti. Abstract: The article deals with the optimization of stable vibration zones of a two-stage gearbox with parametric excitation as a consequence of periodically variable tooth stiffness during meshing. In order to ensure a calm running of the gearbox and to reduce noise, we want to minimize the vibration instability zone in the gearbox design process. The reduction of individual instability zone is largely influenced by the contact ratio of each gear pair and their phase shift. Unfortunately, as one unstable area decreases, other unstable areas can increase. This raises the need to find the optimal solution that can be provided by solving the multi- criteria optimization problem. The article presents the multi-criteria optimization of gear vibration stability zones by using the Weighted Sum Method and the ε-Constraint Method, respectively. The results of multicriteria optimization of stability zones in the article are analytically explained. Both multi-criteria optimization methods are presented in the article for a two-stage, three-gear gearbox with the same meshing frequencies. Key words: vibrations of gear drives, variable mesh stiffness, contact ratio, phase shift, optimization of parametric instability domains. 29 1. Uvod Zmanjševanje amplitud vibracij in hrupa zobniških gonil je eno izmed glavnih vodil v snovanju pogonskih sklopov z zahtevo po mirnejšem in tišjem teku. Osnovni vzrok vibracij in hrupa zobniških gonil predstavlja parametrično vzbujanje kot posledica spremenljive zobne togosti med ubiranjem. Zobna togost, povezana z elastičnim upogibom zob se periodično spreminja zaradi sprememb števila zob, ki so med ubiranjem v kontaktu. Parametrično vzbujanje zobniškega gonila povzroča resonančne pojave s povečevanjem vibracij in nestabilnost gonila [1], [5]. 2. Matematični model vibracij dvostopenjskega zobniškega gonila Dvostopenjsko zobniško gonilo lahko izvedemo v dveh konfiguracijah, prikazanih na sliki 1, ki vsebujeta tri oziroma štiri zobnike. V obeh konfiguracijah upoštevamo le relativne zasuke θ 1 , θ 2 in θ 3 , ki jih merimo glede na zasuke togega telesa na mestu vpetja. Togo telo je lahko pogonski motor, ki se vrti s konstantno vrtilno hitrostjo in je povezan z vhodno gredjo zobniškega gonila. Torzijsko togost vhodne gredi označimo s k L0 . V dvostopenjskem zobniškem gonilu na sliki 1.b predpostavimo, da je vmesna gred, ki povezuje 2. zobnik s 4. zobnikom, toga. S tem se model najbolje približa realnim razmeram, ko se 2. in 4. zobnik v praksi pogosto neposredno dotikata ali sta izdelana v enem kosu. Ozobji zobnikov, ki se ubirata, modeliramo kot linearne vzmeti in njihove togosti označimo s k L1 za prvo, oziroma s k L2 za drugo stopnjo. Te togosti niso konstantne, temveč periodično nihajo okrog določene srednje vrednosti v odvisnosti od tega, kolikšno število zob je trenutno v kontaktu. V nadomestnem modelu so linearne vzmeti pripete na oba vznožna kroga zobnikov, katerih radije označimo z r i , i = 1,2,3,4. S pomočjo momentov vztrajnosti mase I i posameznih gredi s pripadajočimi zobniki tvorimo ekvivalentni masi prvih dveh gredi m i = I i /r i 2 (i = 1,2) ter ekvivalentno maso tretje gredi s pripadajočim zobnikom m 3 = νI 3 /r 3 2 , kjer je ν = 1 pri zobniškem gonilu s tremi in ν = r 4 /r 2 pri zobniškem gonilu s štirimi zobniki. Ekvivalentne togosti izrazimo kot k 0 = k L0 /r 1 2 , k 1 = k L1 in k 2 = ν 2 k L2 . Ekvivalentni zobni togosti k i , i=1,2 izrazimo v obliki vsote k i (t) = k gi + k vi (t), kjer sta k gi srednji vrednosti zobne togosti, k vi (t) pa pripadajoči variabilni komponenti. Njun časovni potek prikazuje slika 2. V nadomestnem modelu vpeljemo namesto zasukov deformacije, ki jih računamo na radijih vznožnih krogov s pomočjo transformacij x 1 =θ 1 r 1 , x 2 = θ 2 r 2 in x 3 = θ 3 r 3 /ν. Z združitvijo deformacij v vektor q = [x 1 , x 2 , x 3 ] T , uvedbo vztrajnostne masne matrike M = diag(m 1 ,m 2 ,m 3 ) in uvedbo togostne matrike K(t) = K 0 + K v (t) kot vsote iz konstantnega dela oziroma matrike srednjih vrednosti togosti K 0 ter variabilnega dela K v (t), lahko zapišemo diferencialno enačbo prostega nihanja zobniškega gonila: , ,(1) kjer je K 0 matrika srednjih vrednosti togosti: , (2) K v (t) pa matrika variabilnega dela zobne togosti: (3) Slika 1. Konfiguraciji dvostopenjskega zobniškega gonila. a) gonilo s tremi zobniki, b) gonilo s štirimi zobniki. Variabilne komponente zobne togosti k vi (t) so v odvisnosti od časa periodične s frekvenco ubiranja Ώ i . Potek zobne togosti lahko dobimo eksperimentalno, vendar bomo v tej raziskavi privzeli pravokotno nihajno obliko kot jo prikazuje slika 2. Pravokotno nihajno obliko variabilne komponente k vi (t) izrazimo s pomočjo okrnjene Fourierove vrste z N členi: (4) V praksi zadostuje, če v Fourierovi vrsti obdržimo le prve tri ali štiri člene. S pomočjo slike 2 lahko izrazimo Fourierove koeficiente a i s in b i s : (5) kjer c i pomenijo kontaktna razmerja, produkti d i T i so fazni premiki, d i označujejo fazne koeficiente, T i = 2 π/ Ω i pa so periode variabilnih komponent zobne togosti. Slika 2. Periodična zobna togost k i (t)=k gi +k vi (t). Simbol c i pomeni kontaktno razmerje, k gi so srednje vrednosti zobne togosti, d i T i pa so fazni koti. 01 1 0 1 122 22 0 0 gg g ggg gg kk k k kkk kk  +  = +     K ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) 11 1122 22 0 0 vv v vvvv vv kt kt t ktktktkt kt kt   = +    K ( ) ( ) 1 2 cos sin N ss vi ai i i i i s kt k a stb st  = = +  ( ) ( ) ( ) ( ) 2 cos 2 sin , 2 sin 2 sin , s i ii i s i ii i a s c d sc s b s c d sc s      = −−   = −−  ( ) 0 v t + +  =  Mq K K q 0 28 Izvirni znanstveni članek TEHNIKA – zobniška gonila Datum prejema: 7. december 2019 ANALI PAZU 10/2020/1-2: 28-35 www.anali-pazu.si Optimalno določanje stabilnih območij vibracij zobniških gonil Optimal determination of stability domains of vibrations in gear drives Rudolf Pušenjak 1 * 1 Fakulteta za industrijski inženiring Novo mesto/Šegova ulica 112, Novo mesto E-mails: rudolf.pusenjak@fini-unm.si * Rudolf Pušenjak Povzetek: V članku je obravnavana optimizacija stabilnih območij vibracij dvostopenjskega zobniškega gonila s parametričnim vzbujanjem kot posledice periodično spremenljive zobne togosti med ubiranjem. Da bi zagotovili mirnejši tek zobniškega gonila in zmanjšali hrup, želimo v procesu snovanja zobniškega gonila območja nestabilnosti vibracij čimbolj zmanjšati. Na zmanjšanje posameznega območja nestabilnosti v največji meri vplivata kontaktno razmerje posamezne zobniške dvojice in njun fazni premik. Na žalost se ob zmanjšanju enega nestabilnega območja lahko povečajo ostala nestabilna območja. S tem se pojavi potreba po iskanju optimalne rešitve, ki jo lahko zagotovimo z rešitvijo večkriterijskega optimizacijskega problema. V članku je predstavljena večkriterijska optimizacija območij stabilnosti vibracij zobniških gonil z metodo utežne vsote in z metodo ε – omejitev. Rezultate večkriterijske optimizacije območij stabilnosti vibracij v članku analitično pojasnimo. Obe večkriterijski optimizacijski metodi sta v članku prikazani pri dvostopenjskemu zobniškemu gonilu s tremi zobniki, pri katerem sta frekvenci ubiranja enaki. Ključne besede: vibracije zobniških gonil, spremenljiva zobna togost, kontaktno razmerje, fazni premik, optimizacija območij parametrične nestabilnosti. Abstract: The article deals with the optimization of stable vibration zones of a two-stage gearbox with parametric excitation as a consequence of periodically variable tooth stiffness during meshing. In order to ensure a calm running of the gearbox and to reduce noise, we want to minimize the vibration instability zone in the gearbox design process. The reduction of individual instability zone is largely influenced by the contact ratio of each gear pair and their phase shift. Unfortunately, as one unstable area decreases, other unstable areas can increase. This raises the need to find the optimal solution that can be provided by solving the multi- criteria optimization problem. The article presents the multi-criteria optimization of gear vibration stability zones by using the Weighted Sum Method and the ε-Constraint Method, respectively. The results of multicriteria optimization of stability zones in the article are analytically explained. Both multi-criteria optimization methods are presented in the article for a two-stage, three-gear gearbox with the same meshing frequencies. Key words: vibrations of gear drives, variable mesh stiffness, contact ratio, phase shift, optimization of parametric instability domains. 29 1. Uvod Zmanjševanje amplitud vibracij in hrupa zobniških gonil je eno izmed glavnih vodil v snovanju pogonskih sklopov z zahtevo po mirnejšem in tišjem teku. Osnovni vzrok vibracij in hrupa zobniških gonil predstavlja parametrično vzbujanje kot posledica spremenljive zobne togosti med ubiranjem. Zobna togost, povezana z elastičnim upogibom zob se periodično spreminja zaradi sprememb števila zob, ki so med ubiranjem v kontaktu. Parametrično vzbujanje zobniškega gonila povzroča resonančne pojave s povečevanjem vibracij in nestabilnost gonila [1], [5]. 2. Matematični model vibracij dvostopenjskega zobniškega gonila Dvostopenjsko zobniško gonilo lahko izvedemo v dveh konfiguracijah, prikazanih na sliki 1, ki vsebujeta tri oziroma štiri zobnike. V obeh konfiguracijah upoštevamo le relativne zasuke θ 1 , θ 2 in θ 3 , ki jih merimo glede na zasuke togega telesa na mestu vpetja. Togo telo je lahko pogonski motor, ki se vrti s konstantno vrtilno hitrostjo in je povezan z vhodno gredjo zobniškega gonila. Torzijsko togost vhodne gredi označimo s k L0 . V dvostopenjskem zobniškem gonilu na sliki 1.b predpostavimo, da je vmesna gred, ki povezuje 2. zobnik s 4. zobnikom, toga. S tem se model najbolje približa realnim razmeram, ko se 2. in 4. zobnik v praksi pogosto neposredno dotikata ali sta izdelana v enem kosu. Ozobji zobnikov, ki se ubirata, modeliramo kot linearne vzmeti in njihove togosti označimo s k L1 za prvo, oziroma s k L2 za drugo stopnjo. Te togosti niso konstantne, temveč periodično nihajo okrog določene srednje vrednosti v odvisnosti od tega, kolikšno število zob je trenutno v kontaktu. V nadomestnem modelu so linearne vzmeti pripete na oba vznožna kroga zobnikov, katerih radije označimo z r i , i = 1,2,3,4. S pomočjo momentov vztrajnosti mase I i posameznih gredi s pripadajočimi zobniki tvorimo ekvivalentni masi prvih dveh gredi m i = I i /r i 2 (i = 1,2) ter ekvivalentno maso tretje gredi s pripadajočim zobnikom m 3 = νI 3 /r 3 2 , kjer je ν = 1 pri zobniškem gonilu s tremi in ν = r 4 /r 2 pri zobniškem gonilu s štirimi zobniki. Ekvivalentne togosti izrazimo kot k 0 = k L0 /r 1 2 , k 1 = k L1 in k 2 = ν 2 k L2 . Ekvivalentni zobni togosti k i , i=1,2 izrazimo v obliki vsote k i (t) = k gi + k vi (t), kjer sta k gi srednji vrednosti zobne togosti, k vi (t) pa pripadajoči variabilni komponenti. Njun časovni potek prikazuje slika 2. V nadomestnem modelu vpeljemo namesto zasukov deformacije, ki jih računamo na radijih vznožnih krogov s pomočjo transformacij x 1 =θ 1 r 1 , x 2 = θ 2 r 2 in x 3 = θ 3 r 3 /ν. Z združitvijo deformacij v vektor q = [x 1 , x 2 , x 3 ] T , uvedbo vztrajnostne masne matrike M = diag(m 1 ,m 2 ,m 3 ) in uvedbo togostne matrike K(t) = K 0 + K v (t) kot vsote iz konstantnega dela oziroma matrike srednjih vrednosti togosti K 0 ter variabilnega dela K v (t), lahko zapišemo diferencialno enačbo prostega nihanja zobniškega gonila: , ,(1) kjer je K 0 matrika srednjih vrednosti togosti: , (2) K v (t) pa matrika variabilnega dela zobne togosti: (3) Slika 1. Konfiguraciji dvostopenjskega zobniškega gonila. a) gonilo s tremi zobniki, b) gonilo s štirimi zobniki. Variabilne komponente zobne togosti k vi (t) so v odvisnosti od časa periodične s frekvenco ubiranja Ώ i . Potek zobne togosti lahko dobimo eksperimentalno, vendar bomo v tej raziskavi privzeli pravokotno nihajno obliko kot jo prikazuje slika 2. Pravokotno nihajno obliko variabilne komponente k vi (t) izrazimo s pomočjo okrnjene Fourierove vrste z N členi: (4) V praksi zadostuje, če v Fourierovi vrsti obdržimo le prve tri ali štiri člene. S pomočjo slike 2 lahko izrazimo Fourierove koeficiente a i s in b i s : (5) kjer c i pomenijo kontaktna razmerja, produkti d i T i so fazni premiki, d i označujejo fazne koeficiente, T i = 2 π/ Ω i pa so periode variabilnih komponent zobne togosti. Slika 2. Periodična zobna togost k i (t)=k gi +k vi (t). Simbol c i pomeni kontaktno razmerje, k gi so srednje vrednosti zobne togosti, d i T i pa so fazni koti. 01 1 0 1 122 22 0 0 gg g ggg gg kk k k kkk kk  +  = +     K ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) 11 1122 22 0 0 vv v vvvv vv kt kt t ktktktkt kt kt   = +    K ( ) ( ) 1 2 cos sin N ss vi ai i i i i s kt k a stb st  = = +  ( ) ( ) ( ) ( ) 2 cos 2 sin , 2 sin 2 sin , s i ii i s i ii i a s c d sc s b s c d sc s      = −−   = −−  ( ) 0 v t + +  =  Mq K K q 0 OPTIMALNO DOLOČANJE STABILNIH OBMOČIJ VIBRACIJ ZOBNIŠKIH GONIL 30 V enačbi (1) lahko masno matriko M odpravimo z rešitvijo časovno invariantnega primera, ki ga opisuje pridruženi problem lastnih vrednosti in lastnih vektorjev ,(6) kjer je Λ = diag(ω 1 2 ,ω 2 2 ,ω 3 2 ) diagonalna matrika lastnih vrednosti, ω i , i = 1,2,3 so lastne (naravne) frekvence, Φ = [Φ 1 ,Φ 2 ,Φ 3 ] pa je matrika lastnih vektorjev, ki ustrezajo naravnim načinom nihanja Φ i , i = 1,2,3. Naravne načine nihanja normiramo z izpolnitvijo zahteve Φ T MΦ = I. Če uporabimo modalno transformacijo vektorja deformacij q = Φu, dobimo z normalizacijo lastnih vektorjev Φ T MΦ = I matrično diferencialno enačbo .(7) Matrični produkt variabilnega dela zobne togosti Φ T K v (t) Φ je ugodno razstaviti na dva dela tako, da posamezni del ustreza i-ti stopnji zobniškega gonila, i = 1,2: (8) kjer sta C 1 in C 2 povezovalni matriki, ki sta določeni na naslednji način: . (9) Z vstavljanjem enačbe (4) v enačbo (8), nato pa z vstavljanjem dobljenega izraza v enačbo (7) dobimo: (10) Za reševanje matrične diferencialne enačbe bomo v nadaljevanju uporabili perturbacijsko metodo z več časovnimi skalami. Ker so amplitude variabilne zobne togosti k ai v praksi majhne nasproti srednjim vrednostim zobne togosti k gi , lahko v splošnem vpeljemo dva majhna, pozitivna perturbacijska parametra ε i =k ai /k gi , i = 1,2 in enačbo (10) prepišemo na obliko: (11) Dobljena enačba je matrična Hillova diferencialna enačba parametrično vzbujenih vibracij (dvostopenjskega) zobniškega gonila [4]. Pomembna lastnost matrik Φ T C i Φ, i = 1,2 je, da sta simetrični. 3. Reševanje matrične diferencialne enačbe parametrično vzbujenih vibracij Reševanje enačbe (11) s perturbacijsko metodo z več časovnimi skalami bomo zaradi prostorskih omejitev v nadaljevanju prikazali za primer ε 1 =ε 2 =ε, pri katerem sta širina zob in material zobnikov enaka in sta tudi amplitudi variabilne zobne togosti obeh zobnikov med ubiranjem enaki. Razlikujeta pa se lahko kontaktni razmerji in fazna kota obeh ozobij. Pri reševanju primera ε 1 ≠ ε 2 so lahko vsi parametri dveh ozobij v posamezni stopnji zobniškega gonila različni. Perturbacijsko metodo uporabimo tako, da poleg običajne časovne skale t vpeljemo še počasno časovno skalo τ = εt, običajni drugi časovni odvod zamenjamo z diferencialnim operatorjem enačbo (11) pa pretvorimo v Hillovo parcialno diferencialno enačbo. V perturbacijski metodi je ugodneje, če uporabimo komponentno obliko enačbe: (12) Izraz (Φ T C i Φ) nr v enačbi (12) pomeni n,r-ti element matrike Φ T C i Φ (to je element na presečišču n-te vrstice in r-tega stolpca matrike). Če vpeljemo naslednje elemente: (13) lahko enačbo (12) zapišemo v obliki: (14) Rešitev enačbe (14) iščemo v obliki potenčne vrste: (15) Z uvrstitvijo nastavka (15) v enačbo (14) ter zbiranjem členov ob istih potencah parametra ε dobimo zaporedje enačb: (16.a) ( ) T v t ++ = u Λu Φ K Φu 0 12 110 000 110 , 011 000 011   = =    CC ( ) ( ) ( ) 22 TT T 11 , 1, 2. v vi vi i ii t t kt i = = = = =  Φ K Φ Φ K Φ Φ CΦ ( )  ( ) 2 T 1 2 T 11 2 cos sin vi i i N ss ai i i i i i is kt k a stb st  = = =  ++ =++     +=    u Λ Φ CΦ u u Λ Φ CΦ u 0 ( ) ( ) 222 22 22 32 T 1 11 2 2 cos sin 0, 1,2,3 nnn nn N ss i i i i gi i r nr r is uuu u t t a stb st k u n      = = =  ++++    +     = = Φ CΦ ( ) ( ) ( ) ( ) TT 1 11 1 11 TT 2 22 2 22 ,, ,, ss ss nr g nr g nr nr ss ss nr g nr g nr nr D ak E bk F ak G bk = = = = Φ CΦ Φ CΦ ΦCΦ ΦCΦ ( ) 222 22 22 3 11 11 22 2 2 cos sin cos sin 0, 1,2,3 nnn nn N ss nr nr rs ss nr nr r uuu u t t D stE st F st G st u n      = =  ++++    ++     += =   ( ) ( ) ( ) 01 0 ,, , , 1, 2, 3 nn n k nk k uut ut ut n     = =++ = =   ( ) 2 T 11 2 cos sin . N ss i i i i i gi i is a stb st k  = = ++  +=    u Λ Φ CΦ u 0 22 2 2 22 2 2 d d 2 t tt       =++ 2 2 d d () t  = 0 = K Φ MΦΛ 2 2 0 0 2 0, n nn u u t   +=  31 (16.b) Splošna rešitev enačbe (16.a) je: (17) Uvrstitev rešitve (17) v enačbo (16.b) nam da: (18) Enačbo (18) lahko ob navedenih pogojih ε 1 = ε 2 = ε uporabljamo pri dvostopenjskih zobniških gonilih s tremi ali s štirimi zobniki. V gonilih s štirimi zobniki sta frekvenci ubiranja Ώ 1 in Ώ 2 različni in medsebojno povezani z relacijo Ώ 1 =RΏ 2 , kjer je R razmerje števila zob 2. in 4. zobnika, R = Z 2 /Z 4 (glej sliko 1.b). Iz slike 1.a pa je razvidno, da je pri gonilih s tremi zobniki R = 1 in sta zato frekvenci ubiranja enaki, Ώ 1 = Ώ 2 = Ώ. V nadaljevanju bomo zaradi prostorske omejitve po- drobno obravnavali samo ta primer, vse ostale primere pa lahko obravnavamo na podoben način. 3.1 Dvostopenjsko gonilo s tremi zobniki in enakimi variacijami zobne togosti Uporaba dveh časovnih skal v perturbacijski meto- di nam omogoča, da raziščemo splošni primer ko se sΏ le malo razlikuje od kombinacijske frekvence ω p + ω q . V ta namen postavimo sΏ = ω p + ω q + εσ, kjer je σ parameter razglasitve, ki ga je potrebno določiti. Ko vstavimo to zamenjavo v desno stran enačbe (18) z Ώ = Ώ 1 = Ώ 2 , nastanejo sekularni členi, ki jih je potrebno eliminirati, da zagotovimo enolično rešitev, Eliminaci- jo sekularnih členov zagotovimo s pogoji rešljivosti: (19) Enačbi (19) imata trivialni rešitvi A p (τ) ≡ 0, A q (τ) ≡ 0. Netrivialni rešitvi iščemo z nastavkoma: , (20) kjer sta a p , a q kompleksni konstanti, λ pa sta korena pridružene karakteristične enačbe: (21) Kadar je v enačbi (21) sta korena karakteristič- ne enačbe realna in negativna, A p in A q pa omejena. Nasprotno sta A p in A q neomejena, kadar je Meji stabilnosti sta tako pri izpolnjenem pogoju ,kar lahko zapišemo v eksplicitni obliki: (22) Če je p = q, imamo primer nestabilnosti vibracij v enojnem načinu nihanja, pri katerem preide enačba (22) na obliko: (23) Enačba mej stabilnosti (23) omejuje območje pri- marne nestabilnosti v primeru, ko je s = 1, območje sekundarne nestabilnosti pa tedaj, ko je s=2. Če v ena- čbi (22) p≠q, vendar je pri tem s = 1, govorimo o kombinacijski nestabilnosti. Z uporabo enačb (22) in (23) izračunamo meje stabilnosti zobniškega gonila s tremi zobniki, ki jih prikazuje diagram na sliki 3. Diagram je izračunan pri naslednjih podatkih: m 1 = 1; m 2 = 0.3; m 3 = 4.0; k 0 = 0.5; k g1 = k g2 = 1; c 1 = c 2 = 1.5; d 1 = d 2 = 0. Analitične rešitve so za meje primarnih nestabilnosti 2ω 1 in 2ω 3 ter kombinacijske nestabilnosti ω 1 + ω 3 preverjene s pomočjo numerične integracije enačbe (11) in Floque- tove teorije [3], pri čemer se rezultati po obeh meto- dah zelo dobro ujemajo na celotnem področju razmer- ja togosti ε = k a /k g . Slika 3. Meje stabilnosti dvostopenjskega zobniškega gonila s tremi zobniki ‒ analitična rešitev s perturbacijsko metodo, *, o numerična rešitev. ( ) ( ) 0 e e , 1,2,3 nn it it nn n uA A n   − =+= ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) 11 11 22 22 2 2 1 1 2 3 11 d 2e d d 2e d n rr rr rr rr n it n n nn n N is t is t s r nr rs is t is t is t is t ss nr nr is t is t it s n nr n A u ui t A De e iE e e F e e A iG e e i                 + −− = = + −− + −− +−− −  += − −    +−          −++          −−+       − ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) 11 11 2 2 22 3 11 rr rr r r rr N is t is t s r nr rs is t is t is t is t ss nr nr is t is t s nr A De e iE e e F e e iG e e          − −+ = = −+ − − −+ −+ −   +          +−++          +−       ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) d d d d 2 0, 2 0 p q A p s s s si q pq pq pq pq A q s s s si p qp qp qp qp i A D F iE G e i A D F iE G e           +  +−+ =   +  +−+ =   ( ) ( ) 22 11 1,2 42 22 1 0, , . pq ss pq pq s ss ss pq pq pq pq pq DF EG           + + = = − −    = +++   2 s pq   1 s p q pq s      = +   ( ) ( ) 1 22 1 2 2. p s p pp s ss ss p pp pp pp pp s DF EG        =     = +++    ( ) 22 3 2 10 11 2 11 1 2 20 2 2 cos sin cos sin N s nn n n nr rs sss nr nr nr r uu u D st t t E stF st G st u    = =   += −−       +++   2 s pq   s pq  =  ( ) ( ) ( ) , i i pp qq A ae A ae     + − = = Anali PAZU, 10/2020/1-2, str. 28-35 Rudolf PUŠENJAK 30 V enačbi (1) lahko masno matriko M odpravimo z rešitvijo časovno invariantnega primera, ki ga opisuje pridruženi problem lastnih vrednosti in lastnih vektorjev ,(6) kjer je Λ = diag(ω 1 2 ,ω 2 2 ,ω 3 2 ) diagonalna matrika lastnih vrednosti, ω i , i = 1,2,3 so lastne (naravne) frekvence, Φ = [Φ 1 ,Φ 2 ,Φ 3 ] pa je matrika lastnih vektorjev, ki ustrezajo naravnim načinom nihanja Φ i , i = 1,2,3. Naravne načine nihanja normiramo z izpolnitvijo zahteve Φ T MΦ = I. Če uporabimo modalno transformacijo vektorja deformacij q = Φu, dobimo z normalizacijo lastnih vektorjev Φ T MΦ = I matrično diferencialno enačbo .(7) Matrični produkt variabilnega dela zobne togosti Φ T K v (t) Φ je ugodno razstaviti na dva dela tako, da posamezni del ustreza i-ti stopnji zobniškega gonila, i = 1,2: (8) kjer sta C 1 in C 2 povezovalni matriki, ki sta določeni na naslednji način: . (9) Z vstavljanjem enačbe (4) v enačbo (8), nato pa z vstavljanjem dobljenega izraza v enačbo (7) dobimo: (10) Za reševanje matrične diferencialne enačbe bomo v nadaljevanju uporabili perturbacijsko metodo z več časovnimi skalami. Ker so amplitude variabilne zobne togosti k ai v praksi majhne nasproti srednjim vrednostim zobne togosti k gi , lahko v splošnem vpeljemo dva majhna, pozitivna perturbacijska parametra ε i =k ai /k gi , i = 1,2 in enačbo (10) prepišemo na obliko: (11) Dobljena enačba je matrična Hillova diferencialna enačba parametrično vzbujenih vibracij (dvostopenjskega) zobniškega gonila [4]. Pomembna lastnost matrik Φ T C i Φ, i = 1,2 je, da sta simetrični. 3. Reševanje matrične diferencialne enačbe parametrično vzbujenih vibracij Reševanje enačbe (11) s perturbacijsko metodo z več časovnimi skalami bomo zaradi prostorskih omejitev v nadaljevanju prikazali za primer ε 1 =ε 2 =ε, pri katerem sta širina zob in material zobnikov enaka in sta tudi amplitudi variabilne zobne togosti obeh zobnikov med ubiranjem enaki. Razlikujeta pa se lahko kontaktni razmerji in fazna kota obeh ozobij. Pri reševanju primera ε 1 ≠ ε 2 so lahko vsi parametri dveh ozobij v posamezni stopnji zobniškega gonila različni. Perturbacijsko metodo uporabimo tako, da poleg običajne časovne skale t vpeljemo še počasno časovno skalo τ = εt, običajni drugi časovni odvod zamenjamo z diferencialnim operatorjem enačbo (11) pa pretvorimo v Hillovo parcialno diferencialno enačbo. V perturbacijski metodi je ugodneje, če uporabimo komponentno obliko enačbe: (12) Izraz (Φ T C i Φ) nr v enačbi (12) pomeni n,r-ti element matrike Φ T C i Φ (to je element na presečišču n-te vrstice in r-tega stolpca matrike). Če vpeljemo naslednje elemente: (13) lahko enačbo (12) zapišemo v obliki: (14) Rešitev enačbe (14) iščemo v obliki potenčne vrste: (15) Z uvrstitvijo nastavka (15) v enačbo (14) ter zbiranjem členov ob istih potencah parametra ε dobimo zaporedje enačb: (16.a) ( ) T v t ++ = u Λu Φ K Φu 0 12 110 000 110 , 011 000 011   = =    CC ( ) ( ) ( ) 22 TT T 11 , 1, 2. v vi vi i ii t t kt i = = = = =  Φ K Φ Φ K Φ Φ CΦ ( )  ( ) 2 T 1 2 T 11 2 cos sin vi i i N ss ai i i i i i is kt k a stb st  = = =  ++ =++     +=    u Λ Φ CΦ u u Λ Φ CΦ u 0 ( ) ( ) 222 22 22 32 T 1 11 2 2 cos sin 0, 1,2,3 nnn nn N ss i i i i gi i r nr r is uuu u t t a stb st k u n      = = =  ++++    +     = = Φ CΦ ( ) ( ) ( ) ( ) TT 1 11 1 11 TT 2 22 2 22 ,, ,, ss ss nr g nr g nr nr ss ss nr g nr g nr nr D ak E bk F ak G bk = = = = Φ CΦ Φ CΦ ΦCΦ ΦCΦ ( ) 222 22 22 3 11 11 22 2 2 cos sin cos sin 0, 1,2,3 nnn nn N ss nr nr rs ss nr nr r uuu u t t D stE st F st G st u n      = =  ++++    ++     += =   ( ) ( ) ( ) 01 0 ,, , , 1, 2, 3 nn n k nk k uut ut ut n     = =++ = =   ( ) 2 T 11 2 cos sin . N ss i i i i i gi i is a stb st k  = = ++  +=    u Λ Φ CΦ u 0 22 2 2 22 2 2 d d 2 t tt       =++ 2 2 d d () t  = 0 = K Φ MΦΛ 2 2 0 0 2 0, n nn u u t   +=  31 (16.b) Splošna rešitev enačbe (16.a) je: (17) Uvrstitev rešitve (17) v enačbo (16.b) nam da: (18) Enačbo (18) lahko ob navedenih pogojih ε 1 = ε 2 = ε uporabljamo pri dvostopenjskih zobniških gonilih s tremi ali s štirimi zobniki. V gonilih s štirimi zobniki sta frekvenci ubiranja Ώ 1 in Ώ 2 različni in medsebojno povezani z relacijo Ώ 1 =RΏ 2 , kjer je R razmerje števila zob 2. in 4. zobnika, R = Z 2 /Z 4 (glej sliko 1.b). Iz slike 1.a pa je razvidno, da je pri gonilih s tremi zobniki R = 1 in sta zato frekvenci ubiranja enaki, Ώ 1 = Ώ 2 = Ώ. V nadaljevanju bomo zaradi prostorske omejitve po- drobno obravnavali samo ta primer, vse ostale primere pa lahko obravnavamo na podoben način. 3.1 Dvostopenjsko gonilo s tremi zobniki in enakimi variacijami zobne togosti Uporaba dveh časovnih skal v perturbacijski meto- di nam omogoča, da raziščemo splošni primer ko se sΏ le malo razlikuje od kombinacijske frekvence ω p + ω q . V ta namen postavimo sΏ = ω p + ω q + εσ, kjer je σ parameter razglasitve, ki ga je potrebno določiti. Ko vstavimo to zamenjavo v desno stran enačbe (18) z Ώ = Ώ 1 = Ώ 2 , nastanejo sekularni členi, ki jih je potrebno eliminirati, da zagotovimo enolično rešitev, Eliminaci- jo sekularnih členov zagotovimo s pogoji rešljivosti: (19) Enačbi (19) imata trivialni rešitvi A p (τ) ≡ 0, A q (τ) ≡ 0. Netrivialni rešitvi iščemo z nastavkoma: , (20) kjer sta a p , a q kompleksni konstanti, λ pa sta korena pridružene karakteristične enačbe: (21) Kadar je v enačbi (21) sta korena karakteristič- ne enačbe realna in negativna, A p in A q pa omejena. Nasprotno sta A p in A q neomejena, kadar je Meji stabilnosti sta tako pri izpolnjenem pogoju ,kar lahko zapišemo v eksplicitni obliki: (22) Če je p = q, imamo primer nestabilnosti vibracij v enojnem načinu nihanja, pri katerem preide enačba (22) na obliko: (23) Enačba mej stabilnosti (23) omejuje območje pri- marne nestabilnosti v primeru, ko je s = 1, območje sekundarne nestabilnosti pa tedaj, ko je s=2. Če v ena- čbi (22) p≠q, vendar je pri tem s = 1, govorimo o kombinacijski nestabilnosti. Z uporabo enačb (22) in (23) izračunamo meje stabilnosti zobniškega gonila s tremi zobniki, ki jih prikazuje diagram na sliki 3. Diagram je izračunan pri naslednjih podatkih: m 1 = 1; m 2 = 0.3; m 3 = 4.0; k 0 = 0.5; k g1 = k g2 = 1; c 1 = c 2 = 1.5; d 1 = d 2 = 0. Analitične rešitve so za meje primarnih nestabilnosti 2ω 1 in 2ω 3 ter kombinacijske nestabilnosti ω 1 + ω 3 preverjene s pomočjo numerične integracije enačbe (11) in Floque- tove teorije [3], pri čemer se rezultati po obeh meto- dah zelo dobro ujemajo na celotnem področju razmer- ja togosti ε = k a /k g . Slika 3. Meje stabilnosti dvostopenjskega zobniškega gonila s tremi zobniki ‒ analitična rešitev s perturbacijsko metodo, *, o numerična rešitev. ( ) ( ) 0 e e , 1,2,3 nn it it nn n uA A n   − =+= ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) 11 11 22 22 2 2 1 1 2 3 11 d 2e d d 2e d n rr rr rr rr n it n n nn n N is t is t s r nr rs is t is t is t is t ss nr nr is t is t it s n nr n A u ui t A De e iE e e F e e A iG e e i                 + −− = = + −− + −− +−− −  += − −    +−          −++          −−+       − ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) 11 11 2 2 22 3 11 rr rr r r rr N is t is t s r nr rs is t is t is t is t ss nr nr is t is t s nr A De e iE e e F e e iG e e          − −+ = = −+ − − −+ −+ −   +          +−++          +−       ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) d d d d 2 0, 2 0 p q A p s s s si q pq pq pq pq A q s s s si p qp qp qp qp i A D F iE G e i A D F iE G e           +  +−+ =   +  +−+ =   ( ) ( ) 22 11 1,2 42 22 1 0, , . pq ss pq pq s ss ss pq pq pq pq pq DF EG           + + = = − −    = +++   2 s pq   1 s p q pq s      = +   ( ) ( ) 1 22 1 2 2. p s p pp s ss ss p pp pp pp pp s DF EG        =     = +++    ( ) 22 3 2 10 11 2 11 1 2 20 2 2 cos sin cos sin N s nn n n nr rs sss nr nr nr r uu u D st t t E stF st G st u    = =   += −−       +++   2 s pq   s pq  =  ( ) ( ) ( ) , i i pp qq A ae A ae     + − = = OPTIMALNO DOLOČANJE STABILNIH OBMOČIJ VIBRACIJ ZOBNIŠKIH GONIL 32 Iz slike 3 je razvidno, da je primarno nestabilno območje z ubirno frekvenco Ω v okolici dvakratne naravne frekvence ω 3, Ω ≈ 2ω 3 bistveno večje kot sta primarni nestabilni območji Ω ≈ 2ω 2 in Ω ≈ 2ω 2 , medtem ko je območje kombinacijske nestabilnosti Ω ≈ ω 2 + ω 3 zanemarljivo majhno. Iz enačb (5), (11), (13), (21) in (22) nedvomno izhaja, da imata na veli- kost nestabilnih območij ključen vpliv kontaktno raz- merje c i in fazni koeficient d i , i = 1,2. Slika 4 prikazuje vpliv kontaktnega razmerja in faznega koeficienta na velikost območja primarne ne- stabilnosti Ω ≈ 2ω 3 in kombinacijske nestabilnosti Ω ≈ ω 2 + ω 3 . Najpomembnejša ugotovitev, ki iz te slike izhaja, je, da ima zmanjšanje primarnega nestabilnega območja Ω ≈ 2ω 3 za posledico povečanje območja kombinacijske nestabilnosti Ω ≈ ω 2 + ω 3 . Slika 4. Vpliv kontaktnih razmerij c 1 in c 2 ter faznih koefici- entov d 1 in d 2 na velikost nestabilnih območij 4 Optimalno določanje stabilnih območij 4.1 Minimizacija posameznega nestabilnega območja Zmanjševanje ploščine trikotnika dosežemo, če zmanjšujemo smerni koeficient mejnih premic. Ker je pri posameznem območju nestabilnosti parameter s določen in tudi parameter ε v postopku konstruiranja fiksno izbiramo, lahko v problemu optimizacije oba parametra vzamemo za konstantna. S tem se izkaže, da je za namensko funkcijo najbolje vzeti kar koeficient . Če dosežemo, da postane smerni koeficient enak nič, s tem pa enak nič tudi koeficient: ,(24) posamezno nestabilno področje izgine. V primeru enojnega načina nihanja, to je v primeru p=q, bo nestabilno območje izginilo, če dosežemo, da velja Elemente (25) lahko s tvorbo matričnih produktov iz- razimo v eksplicitni obliki. Najprej zapišemo matriko lastnih vektorjev s pomočjo elementov: (26) s čimer dobimo matrični produkt v obliki: (27) s pq  — 3 12 12 2 ,0 cc dd = = = = , - - - - 3 1 12 1 2 22 , 0, cc d d = = = = , - . - . - 1 2 12 1.1, 1.9, 0, 0.9 c c dd = = = = . ( ) ( ) 2 22 1 0 p s ss ss pp pp pp pp pp DF EG    = +++=   T , 1,2 i i = Φ CΦ V elementu ( ) T 1 pq Φ CΦ matričnega produkta je z 112 ppp  = + označen relativni odklon zob prvega ozobja pri p-tem modalnem načinu nihanja, z 223 ppp  = + pa relativni odklon zob drugega ozobja pri p-tem modalnem načinu nihanja. Elemente s pq D , s pq E , s pq F in s pq G lahko sedaj s pomočjo formul (25), (27) in Fourierjevih koeficientov (5) izrazimo v obliki: V slikah 3 in 4 moramo upoštevati, da je razmerje togosti ag kk  = v skladu s predpostavko o majhni pozitivni vrednosti perturbacijskega parametra omejeno z določeno vrednostjo parametra max  , npr. max 0.3  = . Glede na to je vsako posamezno območje nestabilnosti geometrično gledano trikotnik z ogliščem v točki ( ) ( ) 1 0, pq s  = = + , mejnima premicama s smernima koeficientoma s pq s    in stranico v smeri vertikale pri max  . ( ) ( ) 22 1 pq s ss ss pq pq pq pq pq DF EG    = +++   ( ) ( ) ( ) ( ) ( ) TT 1 11 1 11 TT 2 22 2 22 ,, ,, , 1,2,3 ss ss pq g pq g pq pq ss ss pq g pq g pq pq D ak E bk F ak G bk pq = = = = = Φ CΦ Φ CΦ ΦCΦ ΦCΦ   11 12 13 1 2 3 21 22 23 31 23 33 ,,         = =    Φ T , 1,2 i i = Φ CΦ ( ) ( ) ( ) ( ) T 1, 1, , , 1,2,3; 1,2 i ip i p iq i q pq ip iq pq i    ++ = ++ = = = Φ CΦ 33 (28) za primer, ko Oglejmo si sedaj, v katerem primeru lahko dosežemo ničelno vrednost koeficienta za in s tem celo odpravo nestabilnega podro- čja. V ta namen enačbo (24) zapišemo v eksplicitni obliki: (29) Podroben pregled enačbe (29) pokaže, da lahko koefi- cient postane enak nič v dveh primerih, ki sta po- membna tudi za večkriterijsko optimizacijo, v kateri želimo zmanjšati ploščino vseh nestabilnih območij naenkrat in ki si jo bomo ogledali kasneje. S tem se je povsem teoretično izkazalo, da popolna odprava (eliminacija) nestabilnega območja določene- ga tipa lahko povzroči maksimalno povečanje nesta- bilnih območij drugih tipov. 4.2 Večkriterijska minimizacija nestabilnih območij Istočasno minimizacijo vseh nestabilnih območij lahko izvedemo, če uporabimo metode večkriterijske optimizacije [2]. Med mnogimi evolucijskimi algorit- mi večkriterijske optimizacije sta najenostavnejša me- toda utežne vsote (Weighted Sum Method) in metoda z ε- omejitvami (ε- Constraint Method). Pri tem je potrebno opozoriti, da ε – omejitev ne smemo zame- njevati s perturbacijskim parametrom ε=k a /k g . Na vsakem posameznem območju nestabilnosti pozna- mo namensko funkcijo ki jo skušamo čim bolj zmanjšati tako, da iščemo optimalne vrednosti kon- taktnih razmerij c i in faznih koeficientov d i , i = 1,2. Pri tem lahko vselej brez škode vzamemo d 1 = 0, tako da se nabor prostih parametrov skrči na konstrukcijske spremenljivke c 1 , c 2 in d 2 . Vzemimo, da imamo mno- žico nestabilnih območij z M elementi, kjer z indek- som m = 1,2,…,M označimo posamezna območja nestabilnosti. Da moremo namenske funkcije med seboj razlikovati, v ta namen vpeljemo novo oznako V splošnem namenske funkcije sestavimo v vektor namenskih funkcij Γ=[Γ 1 , Γ 2 ,…, Γ m ] T . Namenske funkcije niso linearne, niti konveksne, ampak so transcendentalne, zaradi česar se moramo zateči k bolj splošnim metodam nelinearne optimizacije. 4.2.1 Metoda utežne vsote V metodi utežne vsote celotno množico namenskih funkcij skaliramo v eno samo skupno namensko funk- cijo Γ(c 1, c 2, d 2 ) tako, da seštejemo posamezne namen- ske funkcije , , pomnožene z utežmi w m . pq  pq  ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) 22 2 222 1 11 1 22 12 1 2 1 2 2 222 1 21 1 2 2 2 2 2 2 1 4 sin 2cos 2 sin sin sin s ss ss pq pq pq pq pq pq g pq pq gg pq pq g p q DF EG sck s s c c d d sc sc kk sc k             = +++    =   + −− −    +  ( ) 12 2 ,, , s m pq ccd     ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) 1 1 1 11 1 1 1 1 11 1 2 2 2 22 2 2 2 2 22 2 2 cos 2 sin , 2 sin 2 sin , 2 cos 2 sin , 2 sin 2 sin , s pq g p q s pq g p q s pq g p q s pq g p q D s c d sck s E s c d sck s F s c d sc k s G s c d sc k s                 = −−  = −−  = −−  = −−  ( ) ( ) 22 1 pq s ss ss pq pq pq pq pq DF EG    = +++   a) Prvi primer. 0 s pq  = , ko je ( ) ( ) 12 sin sin 0 sc sc  = = , oziroma ko sta 12 , cc celi števili, ki sta enaki 1 ali 2. Ta primer nastopi, ko imamo enojni ali dvojni kontakt parov zob med ubiranjem zobnikov vzdolž celotne kontaktne premice. Hkrati je to edini primer, ko ima koeficient s pq  vrednost nič pri vseh vrstah nestabilnosti, ali povedano drugače, vsa nestabilna območja pri vrednostih c 1 in c 2 enakim 1 ali 2, izginejo. b) Drugi primer. s pq  je enak nič, 0 s pq  = tudi v primeru, ko 12 , cc nista celi števili, vendar pa sta med seboj enaki tako da velja 1 12 ,2 s s cc s + = =  . Ta primer nastopi pri nestabilnostih višjega reda, na primer pri sekundarni nestabilnosti, s=2, terciarni nestabilnosti, s=3, itd. Pri sekundarni nestabilnosti, s=2, je 0 s pq  = , ko imata 12 in cc vrednost 3 12 2 cc = = . Obenem pa lahko pri primarni ali kombinacijski nestabilnosti, pri katerih je s=1, koeficient s pq  zavzame maksimalno vrednost, če je le znak produkta 112 2 pq pq  pozitiven. (Znak produkta 112 2 pq pq  je npr. pozitiven pri podatkih za območja nestabilnosti na sliki 3). s pq  s pq  ( ) 12 2 ,, , s m pq ccd     ( ) 12 2 ,, , s m pq ccd     Anali PAZU, 10/2020/1-2, str. 28-35 Rudolf PUŠENJAK 32 Iz slike 3 je razvidno, da je primarno nestabilno območje z ubirno frekvenco Ω v okolici dvakratne naravne frekvence ω 3, Ω ≈ 2ω 3 bistveno večje kot sta primarni nestabilni območji Ω ≈ 2ω 2 in Ω ≈ 2ω 2 , medtem ko je območje kombinacijske nestabilnosti Ω ≈ ω 2 + ω 3 zanemarljivo majhno. Iz enačb (5), (11), (13), (21) in (22) nedvomno izhaja, da imata na veli- kost nestabilnih območij ključen vpliv kontaktno raz- merje c i in fazni koeficient d i , i = 1,2. Slika 4 prikazuje vpliv kontaktnega razmerja in faznega koeficienta na velikost območja primarne ne- stabilnosti Ω ≈ 2ω 3 in kombinacijske nestabilnosti Ω ≈ ω 2 + ω 3 . Najpomembnejša ugotovitev, ki iz te slike izhaja, je, da ima zmanjšanje primarnega nestabilnega območja Ω ≈ 2ω 3 za posledico povečanje območja kombinacijske nestabilnosti Ω ≈ ω 2 + ω 3 . Slika 4. Vpliv kontaktnih razmerij c 1 in c 2 ter faznih koefici- entov d 1 in d 2 na velikost nestabilnih območij 4 Optimalno določanje stabilnih območij 4.1 Minimizacija posameznega nestabilnega območja Zmanjševanje ploščine trikotnika dosežemo, če zmanjšujemo smerni koeficient mejnih premic. Ker je pri posameznem območju nestabilnosti parameter s določen in tudi parameter ε v postopku konstruiranja fiksno izbiramo, lahko v problemu optimizacije oba parametra vzamemo za konstantna. S tem se izkaže, da je za namensko funkcijo najbolje vzeti kar koeficient . Če dosežemo, da postane smerni koeficient enak nič, s tem pa enak nič tudi koeficient: ,(24) posamezno nestabilno področje izgine. V primeru enojnega načina nihanja, to je v primeru p=q, bo nestabilno območje izginilo, če dosežemo, da velja Elemente (25) lahko s tvorbo matričnih produktov iz- razimo v eksplicitni obliki. Najprej zapišemo matriko lastnih vektorjev s pomočjo elementov: (26) s čimer dobimo matrični produkt v obliki: (27) s pq  — 3 12 12 2 ,0 cc dd = = = = , - - - - 3 1 12 1 2 22 , 0, cc d d = = = = , - . - . - 1 2 12 1.1, 1.9, 0, 0.9 c c dd = = = = . ( ) ( ) 2 22 1 0 p s ss ss pp pp pp pp pp DF EG    = +++=   T , 1,2 i i = Φ CΦ V elementu ( ) T 1 pq Φ CΦ matričnega produkta je z 112 ppp  = + označen relativni odklon zob prvega ozobja pri p-tem modalnem načinu nihanja, z 223 ppp  = + pa relativni odklon zob drugega ozobja pri p-tem modalnem načinu nihanja. Elemente s pq D , s pq E , s pq F in s pq G lahko sedaj s pomočjo formul (25), (27) in Fourierjevih koeficientov (5) izrazimo v obliki: V slikah 3 in 4 moramo upoštevati, da je razmerje togosti ag kk  = v skladu s predpostavko o majhni pozitivni vrednosti perturbacijskega parametra omejeno z določeno vrednostjo parametra max  , npr. max 0.3  = . Glede na to je vsako posamezno območje nestabilnosti geometrično gledano trikotnik z ogliščem v točki ( ) ( ) 1 0, pq s  = = + , mejnima premicama s smernima koeficientoma s pq s    in stranico v smeri vertikale pri max  . ( ) ( ) 22 1 pq s ss ss pq pq pq pq pq DF EG    = +++   ( ) ( ) ( ) ( ) ( ) TT 1 11 1 11 TT 2 22 2 22 ,, ,, , 1,2,3 ss ss pq g pq g pq pq ss ss pq g pq g pq pq D ak E bk F ak G bk pq = = = = = Φ CΦ Φ CΦ ΦCΦ ΦCΦ   11 12 13 1 2 3 21 22 23 31 23 33 ,,         = =    Φ T , 1,2 i i = Φ CΦ ( ) ( ) ( ) ( ) T 1, 1, , , 1,2,3; 1,2 i ip i p iq i q pq ip iq pq i    ++ = ++ = = = Φ CΦ 33 (28) za primer, ko Oglejmo si sedaj, v katerem primeru lahko dosežemo ničelno vrednost koeficienta za in s tem celo odpravo nestabilnega podro- čja. V ta namen enačbo (24) zapišemo v eksplicitni obliki: (29) Podroben pregled enačbe (29) pokaže, da lahko koefi- cient postane enak nič v dveh primerih, ki sta po- membna tudi za večkriterijsko optimizacijo, v kateri želimo zmanjšati ploščino vseh nestabilnih območij naenkrat in ki si jo bomo ogledali kasneje. S tem se je povsem teoretično izkazalo, da popolna odprava (eliminacija) nestabilnega območja določene- ga tipa lahko povzroči maksimalno povečanje nesta- bilnih območij drugih tipov. 4.2 Večkriterijska minimizacija nestabilnih območij Istočasno minimizacijo vseh nestabilnih območij lahko izvedemo, če uporabimo metode večkriterijske optimizacije [2]. Med mnogimi evolucijskimi algorit- mi večkriterijske optimizacije sta najenostavnejša me- toda utežne vsote (Weighted Sum Method) in metoda z ε- omejitvami (ε- Constraint Method). Pri tem je potrebno opozoriti, da ε – omejitev ne smemo zame- njevati s perturbacijskim parametrom ε=k a /k g . Na vsakem posameznem območju nestabilnosti pozna- mo namensko funkcijo ki jo skušamo čim bolj zmanjšati tako, da iščemo optimalne vrednosti kon- taktnih razmerij c i in faznih koeficientov d i , i = 1,2. Pri tem lahko vselej brez škode vzamemo d 1 = 0, tako da se nabor prostih parametrov skrči na konstrukcijske spremenljivke c 1 , c 2 in d 2 . Vzemimo, da imamo mno- žico nestabilnih območij z M elementi, kjer z indek- som m = 1,2,…,M označimo posamezna območja nestabilnosti. Da moremo namenske funkcije med seboj razlikovati, v ta namen vpeljemo novo oznako V splošnem namenske funkcije sestavimo v vektor namenskih funkcij Γ=[Γ 1 , Γ 2 ,…, Γ m ] T . Namenske funkcije niso linearne, niti konveksne, ampak so transcendentalne, zaradi česar se moramo zateči k bolj splošnim metodam nelinearne optimizacije. 4.2.1 Metoda utežne vsote V metodi utežne vsote celotno množico namenskih funkcij skaliramo v eno samo skupno namensko funk- cijo Γ(c 1, c 2, d 2 ) tako, da seštejemo posamezne namen- ske funkcije , , pomnožene z utežmi w m . pq  pq  ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) 22 2 222 1 11 1 22 12 1 2 1 2 2 222 1 21 1 2 2 2 2 2 2 1 4 sin 2cos 2 sin sin sin s ss ss pq pq pq pq pq pq g pq pq gg pq pq g p q DF EG sck s s c c d d sc sc kk sc k             = +++    =   + −− −    +  ( ) 12 2 ,, , s m pq ccd     ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) 1 1 1 11 1 1 1 1 11 1 2 2 2 22 2 2 2 2 22 2 2 cos 2 sin , 2 sin 2 sin , 2 cos 2 sin , 2 sin 2 sin , s pq g p q s pq g p q s pq g p q s pq g p q D s c d sck s E s c d sck s F s c d sc k s G s c d sc k s                 = −−  = −−  = −−  = −−  ( ) ( ) 22 1 pq s ss ss pq pq pq pq pq DF EG    = +++   a) Prvi primer. 0 s pq  = , ko je ( ) ( ) 12 sin sin 0 sc sc  = = , oziroma ko sta 12 , cc celi števili, ki sta enaki 1 ali 2. Ta primer nastopi, ko imamo enojni ali dvojni kontakt parov zob med ubiranjem zobnikov vzdolž celotne kontaktne premice. Hkrati je to edini primer, ko ima koeficient s pq  vrednost nič pri vseh vrstah nestabilnosti, ali povedano drugače, vsa nestabilna območja pri vrednostih c 1 in c 2 enakim 1 ali 2, izginejo. b) Drugi primer. s pq  je enak nič, 0 s pq  = tudi v primeru, ko 12 , cc nista celi števili, vendar pa sta med seboj enaki tako da velja 1 12 ,2 s s cc s + = =  . Ta primer nastopi pri nestabilnostih višjega reda, na primer pri sekundarni nestabilnosti, s=2, terciarni nestabilnosti, s=3, itd. Pri sekundarni nestabilnosti, s=2, je 0 s pq  = , ko imata 12 in cc vrednost 3 12 2 cc = = . Obenem pa lahko pri primarni ali kombinacijski nestabilnosti, pri katerih je s=1, koeficient s pq  zavzame maksimalno vrednost, če je le znak produkta 112 2 pq pq  pozitiven. (Znak produkta 112 2 pq pq  je npr. pozitiven pri podatkih za območja nestabilnosti na sliki 3). s pq  s pq  ( ) 12 2 ,, , s m pq ccd     ( ) 12 2 ,, , s m pq ccd     OPTIMALNO DOLOČANJE STABILNIH OBMOČIJ VIBRACIJ ZOBNIŠKIH GONIL 34 Skalarizirana namenska funkcija je tedaj: (30) V literaturi lahko najdemo številne predloge za izbiro uteži, med njimi predlog ki pa ni obvezen. Uteži w m zato navadno izbiramo v skladu s preferenca- mi, oziroma prioritetami posameznih namenskih funk- cij. Če nobeni namenski funkciji ne pripišemo posebne vloge, izberemo uteži, ki so enake w m =1, m=1,2,…,M, kar smo uporabili v zgledu, prikazanem v nadaljeva- nju. Večkriterijski optimizacijski problem za istočasno minimizacijo nestabilnih območij z metodo utežne vsote formuliramo v obliki: (31) kjer so d 1 =0, f(d 2 )=0 in g i (c i )=0, i=1,2 omejitvene ena- čbe, ki jih morajo izpolnjevati konstrukcijske spremen- ljivke c 1 , c 2 , d 1 in d 2 , omejitvene neenačbe 1≤c i ≤2, i=1,2 in 0≤d 2 ≤1 pa določajo intervale dopustnih vre- dnosti konstrukcijskih spremenljivk. Ker namenske funkcije niso konveksne, tudi skalarizirana skupna namenska funkcija ni konveksna. Metoda utežne vsote je enostavna, vendar je njena uporabnost zaradi nekonveksne skalarizirane namen- ske funkcije omejena, ker obstaja možnost, da nekatere optimalne rešitve v Paretovem smislu v nekonveksnem ciljnem prostoru ne bodo izračunane. 4.2.2 Metoda ε – omejitev Za razliko od metode utežne vsote, lahko metodo ε - omejitev uporabljamo tako v konveksnih kot tudi v nekonveksnih problemih. V metodi ε- omejitev za na- mensko funkcijo izberemo eno izmed funkcij m=1,2, …,M, npr. funkcijo vse ostale funkcije m=1,2,…,M, m≠μ pa omejimo z omejitvami ε m , ki jih sami predpi- šemo. V metodi ε- omejitev za namensko funkcijo izbere- mo eno izmed funkcij , m=1,2,…,M, npr. funkcijo , vse ostale funkcije , m=1,2,…,M, m≠μ pa omejimo z omejitvami ε m , ki jih sami predpišemo. Predpisovanje ε – omejitev predstavlja težji del optimizacijskega problema, ker ne obstaja nobeno enotno objektivno merilo, s katerim bi omejitvene vre- dnosti ε m lahko določili. V primeru večkriterijske mini- mizacije nestabilnih območij vibracij zobniških gonil lahko ta problem zadovoljivo rešimo na tak način, da na posameznih nestabilnih območjih izvedemo lokalno optimizacijo in v njej določimo minimalno, oziroma maksimalno vrednost , ,m≠μ. Vrednost ε m nato izračunamo po pravilu: (32) kjer je h nek prost parameter, ki si ga izberemo. Več- kriterijski optimizacijski problem, v katerem nastopajo enake omejitvene funkcije in enaki intervali dopustnih vrednosti konstrukcijskih spremenljivk kot v metodi utežne vsote, z metodo ε- omejitev zapišemo v nasled- nji formalni obliki: (33) 4.2.3 Primerjava rezultatov Kot zgled večkriterijske optimizacije nestabilnih območij z metodo utežne vsote in z metodo ε – omeji- tev si oglejmo minimizacijo nestabilnih območij Ω ≈ 2ω 3 in Ω ≈ ω 2 + ω 3 , ki sta prikazani na sliki 4. Primer- javo rezultatov v obeh metodah prikazujeta Tabela 1 in 2. V obeh tabelah primerjamo rezultate minimizacije območij nestabilnosti v treh primerih, v katerih je vre- dnost faznega koeficienta d 1 vselej enaka nič. V prvem primeru so preostali parametri c 1 , c 2 in d 2 omejeni z intervali dopustnih vrednosti, znotraj teh intervalov pa lahko zavzamejo poljubne vrednosti. Rezultat optimi- zacije je pričakovan, kontaktni razmerji v obeh meto- dah zavzameta celoštevilčni vrednosti 1 oziroma 2, namenski funkciji pa sta enaki nič v skladu s primerom a) v razdelku 4.1. Optimizacija v prvem primeru je popolna, ker obe nestabilni območji v celoti izgineta. Rešitev je optimalna v Paretovem smislu, ker ne obstaja nobena druga rešitev, ki bi nad to rešitvijo dominirala. V drugem primeru razen faznega koefici- enta d 1 =0 predpišemo še vrednost kontaktnega raz- merja c 1 = 1.5, v minimizaciji nestabilnih območij pa ( ) ( ) 12 2 12 2 1 ,, , ,, M s m m pq m ccd w ccd  =  =    ( ) 12 2 ,, , s m pq ccd     ( ) 12 2 ,, , s pq ccd      ( ) 12 2 ,, , s m pq ccd     ( ) ( ) ( ) ( ) 12 2 12 2 1 1 2 2 ,, ,, , , 0, 0, 0, 1, 2 0 1, 1 2, 1, 2 M s m m pq m ii i ccd w ccd d fd gc i d ci   =  =    =   =   = =      =   minimiziraj z omejitvami ( ) 12 2 ,, ccd  ( ) ( ) ( ) ( ) 12 2 12 2 1 2 2 ,, , ,, , , 1,2, , in kjer 0, 0, 0, 1, 2 0 1, 1 2, 1, 2 s pq s m pq m ii i ccd ccd m Mm d fd gc i d ci             =    =   =  = =       =  minimiziraj z omejitvami 11 11 12 in .  1 1, M m m w = =  ( ) 12 2 ,, , s m pq ccd     ( ) 12 2 ,, , s m pq ccd     ( ) ( ) ( ) ( ) 1 12 2 2 12 2 min , , , max , , , , s m m pq s m pq h ccd ccd m     = +            35 iščemo optimalne vrednosti kontaktnega razmerja c 2 in faznega koeficienta d 2 . Rezultat minimizacije se v obeh metodah ujema z izsledki primera b) v razdelku 4.1, to je z enakima vrednostima kontaktnega razmer- ja c 1 = c 2 =1.5 in vrednostjo faznega koeficienta d 2 =0.5 ob maksimalnih vrednostih namenskih funkcij Nestabilni območji se zaradi tega prav nič ne zmanj- šata in rešitev ni optimalna v Paretovem smislu. Tabela 1. Večkriterijska optimizacija po metodi utežne vsote Tabela 2. Večkriterijska optimizacija po metodi ε – omejitev V tretjem primeru predpišemo vrednosti obeh faznih koeficientov in sicer d 1 =0, d 2 =0.5, v postopku mini- mizacije nestabilnih območij Ω≈2 ω 3 in Ω≈ ω 2 + ω 3 pa iščemo optimalni vrednosti obeh kontaktnih razmerij c 1 in c 2 . V obeh metodah večkriterijske optimizacije dobimo optimalni vrednosti c 1 =c 2 =2 in vrednosti na- menskih funkcij Rezultat nam pove, da obe nestabilni območji izgineta in ustreza primeru a) iz razdelka 4.1. Rešitev je opti- malna v Paretovem smislu. 5. Zaključek V raziskavi smo razvili sistematični matematični model za analizo stabilnosti dvostopenjskih zobniških gonil kot posledice časovno spremenljive zobne togo- sti. Meje stabilnosti smo s predstavljenim modelom dobili v analitični obliki. V članku smo raziskali po- goje, ob katerih lahko odpravimo vsa območja nesta- bilnosti, prav tako pa smo našli pogoje, pri katerih odprava enega območja nestabilnosti povzroči pove- čanje drugih nestabilnih območij. S pomočjo teoretič- nih izsledkov smo potrdili rezultate večkriterijske optimizacije z metodo utežne vsote in z metodo ε – omejitev v treh karakterističnih primerih. Zahvala Avtor se zahvaljuje Javni Agenciji za raziskovalno dejavnost Republike Slovenije (ARRS), za podporo tega dela v okviru temeljne raziskave J2-9224, Mini- strstvu za izobraževanje, znanost in šport Republike Slovenije in Evropskemu skladu za regionalni razvoj. Literatura 1. M. Benton, A. Seireg, Simulation of Resonan- ces and Instability Conditions in Pinion-Gear Systems, ASME J. Mech. Des.. 100, 26--30, 1978. 2. K. Deb, Multi-Objective Optimization using Evolutionary Algorithms, John Wiley&Sons, Inc., 2001. 3. P. Friedmann, C.E Hammond, T-H Woo, Effi- cient numerical treatment of periodic systems with application to stability problems, Int. J. numer. Meth. Engng. 11, 1117-1136,1977. 4. R. Pušenjak, Nestabilnost dvostopenjskih zob- niških gonil kot posledica spremenljive zobne togosti v ubiranju, Slovensko društvo za meha- niko, Kuhljevi dnevi 19, Strunjan, 26.-27. sept 2019, 115-122. 5. G. V. Tordion, R. Gauvin, Dynamic Stability of a Two-Stage Gear Train Under the Influen- ce of Variable Meshing Stiffness, ASME J. Eng. Ind. 99, 785--791, 1977. 11 11 12 in  11 11 12 0  = = Anali PAZU, 10/2020/1-2, str. 28-35 Rudolf PUŠENJAK 34 Skalarizirana namenska funkcija je tedaj: (30) V literaturi lahko najdemo številne predloge za izbiro uteži, med njimi predlog ki pa ni obvezen. Uteži w m zato navadno izbiramo v skladu s preferenca- mi, oziroma prioritetami posameznih namenskih funk- cij. Če nobeni namenski funkciji ne pripišemo posebne vloge, izberemo uteži, ki so enake w m =1, m=1,2,…,M, kar smo uporabili v zgledu, prikazanem v nadaljeva- nju. Večkriterijski optimizacijski problem za istočasno minimizacijo nestabilnih območij z metodo utežne vsote formuliramo v obliki: (31) kjer so d 1 =0, f(d 2 )=0 in g i (c i )=0, i=1,2 omejitvene ena- čbe, ki jih morajo izpolnjevati konstrukcijske spremen- ljivke c 1 , c 2 , d 1 in d 2 , omejitvene neenačbe 1≤c i ≤2, i=1,2 in 0≤d 2 ≤1 pa določajo intervale dopustnih vre- dnosti konstrukcijskih spremenljivk. Ker namenske funkcije niso konveksne, tudi skalarizirana skupna namenska funkcija ni konveksna. Metoda utežne vsote je enostavna, vendar je njena uporabnost zaradi nekonveksne skalarizirane namen- ske funkcije omejena, ker obstaja možnost, da nekatere optimalne rešitve v Paretovem smislu v nekonveksnem ciljnem prostoru ne bodo izračunane. 4.2.2 Metoda ε – omejitev Za razliko od metode utežne vsote, lahko metodo ε - omejitev uporabljamo tako v konveksnih kot tudi v nekonveksnih problemih. V metodi ε- omejitev za na- mensko funkcijo izberemo eno izmed funkcij m=1,2, …,M, npr. funkcijo vse ostale funkcije m=1,2,…,M, m≠μ pa omejimo z omejitvami ε m , ki jih sami predpi- šemo. V metodi ε- omejitev za namensko funkcijo izbere- mo eno izmed funkcij , m=1,2,…,M, npr. funkcijo , vse ostale funkcije , m=1,2,…,M, m≠μ pa omejimo z omejitvami ε m , ki jih sami predpišemo. Predpisovanje ε – omejitev predstavlja težji del optimizacijskega problema, ker ne obstaja nobeno enotno objektivno merilo, s katerim bi omejitvene vre- dnosti ε m lahko določili. V primeru večkriterijske mini- mizacije nestabilnih območij vibracij zobniških gonil lahko ta problem zadovoljivo rešimo na tak način, da na posameznih nestabilnih območjih izvedemo lokalno optimizacijo in v njej določimo minimalno, oziroma maksimalno vrednost , ,m≠μ. Vrednost ε m nato izračunamo po pravilu: (32) kjer je h nek prost parameter, ki si ga izberemo. Več- kriterijski optimizacijski problem, v katerem nastopajo enake omejitvene funkcije in enaki intervali dopustnih vrednosti konstrukcijskih spremenljivk kot v metodi utežne vsote, z metodo ε- omejitev zapišemo v nasled- nji formalni obliki: (33) 4.2.3 Primerjava rezultatov Kot zgled večkriterijske optimizacije nestabilnih območij z metodo utežne vsote in z metodo ε – omeji- tev si oglejmo minimizacijo nestabilnih območij Ω ≈ 2ω 3 in Ω ≈ ω 2 + ω 3 , ki sta prikazani na sliki 4. Primer- javo rezultatov v obeh metodah prikazujeta Tabela 1 in 2. V obeh tabelah primerjamo rezultate minimizacije območij nestabilnosti v treh primerih, v katerih je vre- dnost faznega koeficienta d 1 vselej enaka nič. V prvem primeru so preostali parametri c 1 , c 2 in d 2 omejeni z intervali dopustnih vrednosti, znotraj teh intervalov pa lahko zavzamejo poljubne vrednosti. Rezultat optimi- zacije je pričakovan, kontaktni razmerji v obeh meto- dah zavzameta celoštevilčni vrednosti 1 oziroma 2, namenski funkciji pa sta enaki nič v skladu s primerom a) v razdelku 4.1. Optimizacija v prvem primeru je popolna, ker obe nestabilni območji v celoti izgineta. Rešitev je optimalna v Paretovem smislu, ker ne obstaja nobena druga rešitev, ki bi nad to rešitvijo dominirala. V drugem primeru razen faznega koefici- enta d 1 =0 predpišemo še vrednost kontaktnega raz- merja c 1 = 1.5, v minimizaciji nestabilnih območij pa ( ) ( ) 12 2 12 2 1 ,, , ,, M s m m pq m ccd w ccd  =  =    ( ) 12 2 ,, , s m pq ccd     ( ) 12 2 ,, , s pq ccd      ( ) 12 2 ,, , s m pq ccd     ( ) ( ) ( ) ( ) 12 2 12 2 1 1 2 2 ,, ,, , , 0, 0, 0, 1, 2 0 1, 1 2, 1, 2 M s m m pq m ii i ccd w ccd d fd gc i d ci   =  =    =   =   = =      =   minimiziraj z omejitvami ( ) 12 2 ,, ccd  ( ) ( ) ( ) ( ) 12 2 12 2 1 2 2 ,, , ,, , , 1,2, , in kjer 0, 0, 0, 1, 2 0 1, 1 2, 1, 2 s pq s m pq m ii i ccd ccd m Mm d fd gc i d ci             =    =   =  = =       =  minimiziraj z omejitvami 11 11 12 in .  1 1, M m m w = =  ( ) 12 2 ,, , s m pq ccd     ( ) 12 2 ,, , s m pq ccd     ( ) ( ) ( ) ( ) 1 12 2 2 12 2 min , , , max , , , , s m m pq s m pq h ccd ccd m     = +            35 iščemo optimalne vrednosti kontaktnega razmerja c 2 in faznega koeficienta d 2 . Rezultat minimizacije se v obeh metodah ujema z izsledki primera b) v razdelku 4.1, to je z enakima vrednostima kontaktnega razmer- ja c 1 = c 2 =1.5 in vrednostjo faznega koeficienta d 2 =0.5 ob maksimalnih vrednostih namenskih funkcij Nestabilni območji se zaradi tega prav nič ne zmanj- šata in rešitev ni optimalna v Paretovem smislu. Tabela 1. Večkriterijska optimizacija po metodi utežne vsote Tabela 2. Večkriterijska optimizacija po metodi ε – omejitev V tretjem primeru predpišemo vrednosti obeh faznih koeficientov in sicer d 1 =0, d 2 =0.5, v postopku mini- mizacije nestabilnih območij Ω≈2 ω 3 in Ω≈ ω 2 + ω 3 pa iščemo optimalni vrednosti obeh kontaktnih razmerij c 1 in c 2 . V obeh metodah večkriterijske optimizacije dobimo optimalni vrednosti c 1 =c 2 =2 in vrednosti na- menskih funkcij Rezultat nam pove, da obe nestabilni območji izgineta in ustreza primeru a) iz razdelka 4.1. Rešitev je opti- malna v Paretovem smislu. 5. Zaključek V raziskavi smo razvili sistematični matematični model za analizo stabilnosti dvostopenjskih zobniških gonil kot posledice časovno spremenljive zobne togo- sti. Meje stabilnosti smo s predstavljenim modelom dobili v analitični obliki. V članku smo raziskali po- goje, ob katerih lahko odpravimo vsa območja nesta- bilnosti, prav tako pa smo našli pogoje, pri katerih odprava enega območja nestabilnosti povzroči pove- čanje drugih nestabilnih območij. S pomočjo teoretič- nih izsledkov smo potrdili rezultate večkriterijske optimizacije z metodo utežne vsote in z metodo ε – omejitev v treh karakterističnih primerih. Zahvala Avtor se zahvaljuje Javni Agenciji za raziskovalno dejavnost Republike Slovenije (ARRS), za podporo tega dela v okviru temeljne raziskave J2-9224, Mini- strstvu za izobraževanje, znanost in šport Republike Slovenije in Evropskemu skladu za regionalni razvoj. Literatura 1. M. Benton, A. Seireg, Simulation of Resonan- ces and Instability Conditions in Pinion-Gear Systems, ASME J. Mech. Des.. 100, 26--30, 1978. 2. K. Deb, Multi-Objective Optimization using Evolutionary Algorithms, John Wiley&Sons, Inc., 2001. 3. P. Friedmann, C.E Hammond, T-H Woo, Effi- cient numerical treatment of periodic systems with application to stability problems, Int. J. numer. Meth. Engng. 11, 1117-1136,1977. 4. R. Pušenjak, Nestabilnost dvostopenjskih zob- niških gonil kot posledica spremenljive zobne togosti v ubiranju, Slovensko društvo za meha- niko, Kuhljevi dnevi 19, Strunjan, 26.-27. sept 2019, 115-122. 5. G. V. Tordion, R. Gauvin, Dynamic Stability of a Two-Stage Gear Train Under the Influen- ce of Variable Meshing Stiffness, ASME J. Eng. Ind. 99, 785--791, 1977. 11 11 12 in  11 11 12 0  = = OPTIMALNO DOLOČANJE STABILNIH OBMOČIJ VIBRACIJ ZOBNIŠKIH GONIL