Univerza v Ljubljani Fakulteta za elektrotehniko Edi Bulic Izračun kvazistatičnega magnetnega polja ob dolgih kovinskih zaslonih z metodo robnih elementov DOKTORSKA DISERTACIJA Mentor: doc. dr. Anton R. Sinigoj Ljubljana, 2009 Kazalo Povzetek v Abstract vii 1 Uvod 1 1.1 Zaslanjanje magnetnega polja .......................... 1 1.2 Obravnavane zaslonske strukture ........................ 3 1.2.1 Vrednotenje uˇcinka zastiranja ...................... 4 1.3 Metoda robnih elementov ............................ 4 2 Koncept ekvivalentnih ploskovnih virov 6 2.1 Robni integralski enaˇcbi ............................. 6 2.1.1 Matriˇcna oblika robnih integralskih enaˇcb ............... 10 2.2 Sistem sklopljenih integralskih enaˇcb ...................... 11 2.3 Numeriˇcno reˇsevanje sistema integralskih enaˇcb ................ 12 2.3.1 Lastni prispevki ............................. 16 2.4 Izraˇcun magnetnega polja ............................ 19 2.4.1 Izraˇcun vektorskega magnetnega potenciala .............. 20 2.4.2 Izraˇcun vektorja gostote magnetnega pretoka ............. 22 2.4.3 Izraˇcun joulske moˇci v cilindru ..................... 24 3 Verifikacija predlagane metode 25 3.1 Cevasti zaslon kroˇznega preseka ......................... 25 3.1.1 Vodnik v koaksialnem zaslonu ...................... 25 3.1.2 Oklopljen dvovod ............................. 28 3.2 Zelo ˇsirok raven zaslon .............................. 32 3.3 Raven zaslon ................................... 33 3.3.1 Hibridna metoda (HM) ......................... 34 iii Kazalo 3.3.2 Metoda tokovnih niti (MTN) ...................... 34 3.3.3 Metoda ploskovnih tokov (MPT) .................... 35 3.3.4 Meritev .................................. 37 4 Analiza zaslonskih struktur 39 4.1 Enofazno vzbujanje ................................ 40 4.2 Trifazno vzbujanje ................................ 41 4.3 Trifazno vzbujanje in dvozaslonska struktura ................. 42 4.4 Analiza povratnega uˇcinka ............................ 42 5 Sklep 45 Prispevkikznanosti 46 Izjava 48 Zahvala 49 Literatura 57 Izvirni znanstveni ˇclanek 59 iv Povzetek Doktorska disertacija se ukvarja z metodo, ki omogoˇca analizo uˇcinkovitosti zastiranja magnetnega polja s kovinskimi zasloni v 2D strukturah. Viri primarnega polja so harmoniˇcni toki v ravnih vzporednih vodnikih. To polje v delih prostora zastrejo vodnikom vzporedni prevodni in/ali magnetni zasloni. V zaslonih se inducirajo elektriˇcni toki in v njih se odvija magnetenje. Oba procesa imata reakcijski vpliv na primarno polje – da se slednje s sekundarnim dopolni v novo, rezultanˇcno polje. V disertaciji privzemamo vodnike in zaslone kot linearne cilindre znanih specifiˇcnih elektriˇcnih prevodnosti in permeabilnosti, okoliˇski prostor kot prazen, preˇcne izmere zaslonske strukture pa kot majhne v primerjavi z njeno dolˇzino, da je upraviˇcena 2D analiza, in neznatne v primerjavi z valovno dolˇzino, da je zadoˇsˇceno pogoju kvazistatiˇcnosti. Vektorski magnetni potencial v in ob vzporednih prevodnih in/ali magnetnih zaslonih in vodnikih je zapisan kot vsota polj linijskih virov in oblog enojnega in dvojnega sloja na povrˇsinah vodnikov in zaslonov. Ti oblogi se izraˇcuna iz sistema sklopljenih robnih integralskih enaˇcb (2.18), izpeljanih z uporabo lastnosti potencialov enojnega in dvojnega sloja, ter enaˇcb tipa (2.21), (2.22), (2.23) in/ali (2.24), ki zajemajo napetostno-tokovne razmere v vodnikih in zaslonih. Za numeriˇcni izraˇcun sistema enaˇcb sta uporabljeni metoda robnih elementov oziroma momentna metoda. V nadaljevanju je podana vrsta verifikacij predlagane metode. Na konkretnih primerih so predstavljene primerjave rezultatov predlagane metode z rezultati analitiˇcne metode, hibridne metode, metode tokovnih niti pri nemagnetnih zaslonih, metode ploskovnih tokov pri neprevodnih magnetnih zaslonih in primerjave z rezultati razliˇcnih meritev. Pri tem je pomembna predvsem krajevna porazdelitev rezultanˇcnega polja ali faktorja zastiranja, ki je definiran kot razmerje absolutnih vrednosti gostote magnetnega pretoka primarnega in rezultanˇcnega polja. Analitiˇcni reˇsitvi za primer vodnika ali dvovoda v cevastem zaslonu kroˇznega preseka omogoˇcata oceno natanˇcnosti predlagane metode. Iz primerjav na slikah 3.2, 3.3 in 3.6 v Povzetek je razvidno odlično ujemanje rezultatov predlagane metode z analitičnima; edina omembe vredna odstopanja so tik ob površinah zaslona. Analitični rešitvi omogočata tudi oceno vpliva velikosti robnih elementov na natančnost predlagane metode (sliki 3.4 in 3.7). Analitična rešitev obstoja tudi za neskončno širok raven zaslon. Z njo so primerjani rezultati predlagane metode za primer zelo širokega zaslona. Iz diagramov na sliki 3.10 je razvidno dobro ujemanje ob sredini zaslona, odstopanje pa pričakovano narašča s približevanjem krajiščem zaslona. Pri končno širokem ravnem nemagnetnem zaslonu je opravljena primerjava z rezultati meritev iz literature, z rezultati hibridne metode in metode tokovnih niti (slika 3.12). Ujemanje je zelo dobro pri slednji, ujemanje z merilnimi rezultati pa je dokaj dobro. Pri magnetnem zaslonu je primerjava opravljena z rezultati hibridne metode ter metode ploskovnih tokov (slika 3.13). Pri zadnji je ujemanje pri zaslonih, ki so tudi prevodni, frekvenčno zelo odvisno: pri nizkih frekvencah je izvrstno, pri višjih pa ne, saj metoda ploskovnih tokov ne upošteva induciranih tokov. V sklepnem delu je podanih nekaj izračunov zastiranja, ki nakažejo tudi možnosti uporabe predlagane metode. Posebej zanimiva je možnost analize povratnega učinka na primarno strukturo, ki je pri problemih zastiranja v literaturi ni zaslediti. Povratni učinek se kaže predvsem v povečanju notranjih izgub primarnih vodov in je izdatnejši, če so zasloni bliže primarnih virov. Izbrali smo zaslon profila U, ki je prevoden nemagneten ali neprevoden magneten ali prevoden magneten in zastira polje tokovodnika ali dvovoda ali polje trifaznega voda. Na prikazu magnetnih polj okrog zaslona pri tokovodniku na sliki 4.2 in dvovodu na sliki 4.5 sta lepo razvidna osnovna principa zastiranja, indukcija in magnetizacija, ter njun kombiniran učinek. Učinkovitost zaslonov je prikazana na sliki 4.3, v primeru zaslanjanja polja trifaznega voda, ter podobno na sliki 4.4, kjer dodatni vodoravni del sestavljenega zaslona skoraj zapre odprtino profila U. V bližini zaslona je boljši učinek magnetnega, dlje stran prevodnega, povsod pa je od obeh izdatno učinkovitejši prevodni magnetni. Pri zadnjem je še najbolj očitno izboljšanje, ki ga prinese dodatni vodoravni del. vi Abstract In this dissertation we present a method that enables the analysis of the magnetic field shielding efficiency of metal shields in 2D structures. The primary magnetic field source is time-harmonic currents in straight parallel conductors. This field is shielded in some regions by conducting and/or magnetic shields that run parallel to the conductors. Eddy currents are induced in the shields that can also become magnetized. Both processes impact the primary field. The resultant field is a superposition of the primary field and the secondary field. We assume that: 1) the conductors and shields are linear cylinders with known electric conductivities and magnetic permeabilities; 2) the surrounding medium is empty; 3) transverse dimensions of the structure which are small compared to the structure length justify a 2D approach; and 4) these dimensions are negligible compared to the wave length and therefor fulfil quasi-static conditions. The magnetic vector potential inside and outside the parallel conducting and/or magnetic shields and conductors is expressed by a sum of fields caused by line sources and by equivalent single- and double-layer surface sources on the conductors and the shields. A system of coupled integral equations ((2.18), (2.21), (2.22), (2.23), and/or (2.24)) for unknown distributions of these equivalent sources is derived by implementing the properties of single- and double-layer potentials, and by taking into account voltage-current relations in the conductors and the shields. This system is solved by using the boundary element method and the moment method. Several verifications of the proposed method are given afterwards. The results obtained with this method on concrete examples are compared with the results of: 1) analytically solvable problems; 2) hybrid method; 3) multiconductor method in the case of the nonmagnetic shields; 4) surface-current method for analysis of the nonconducting magnetic shields; and 5) measurements. The space distributions of the resultant field and shielding factor (that is defined as a ratio between the primary and the resultant magnetic flux density magnitude) are the focus of these comparisons. vii Abstract In the case of a conductor or bifilar pair of parallel conductors inside a cylindrical shield analytical solution can be used as an estimation of the accuracy of the proposed method. The comparisons in Figures 3.2, 3.3, and 3.6 clearly show that the results obtained with the proposed method are in great accordance with the analytical ones; the only significant deviation is detected close to the shield surfaces. These analytical solutions also enable the estimation of the impact of the length of the boundary elements on the accuracy of the proposed method (Figures 3.4 and 3.7). The analytical solution is also possible in the case of an infinite planar shield. The results of the proposed method in the case of an exceedingly wide shield are compared with the results of the analytical solution. From the diagrams in Figure 3.10 it is clear that: 1) good accordance in the vicinity of the shield centre is obtained, and 2) discrepancy increases rapidly when approaching either end of the shield. In the case of the finite-width flat nonmagnetic shield (Figure 3.12) we compare the results with: 1) measurements reported in literature; 2) the hybrid method; and 3) the multiconductor method. The best agreement is achieved with the multiconductor method. The agreement with the measurements is also good. In the case of the magnetic shield the comparison is made with the results of the hybrid method and the surface-current method (Figure 3.13). The agreement with the surface-current method when the shield is conducting depends on the frequency: at low frequencies it is excellent, but at higher frequencies it is not because this method does not consider the induced eddy currents. Some computations of shielding that demonstrate the possibilities of usage of the proposed method are given in the final part. The possibility of the analysis of the reactive effect on the primary structure, which has not yet been considered in literature for the problems of shielding, can be of great interest. The reactive effect is evident in the increase of the power dissipation and even more if the shields are closer to the conductors. The shield of a U-shape is chosen. It is either nonmagnetic conducting, magnetic nonconducting, or magnetic conducting, and it shields the field of either one conductor, a pair of conductors, or three-phase conductors. The presentation of the magnetic field around the shield in the vicinity of the conductor (Figure 4.2) or the pair of conductors (Figure 4.5) shows two basic principles of shielding, i.e. induction and magnetization, and their combined effect. The efficiency of shielding is demonstrated in Figure 4.3 for the case of three-phase system, while in Figure 4.4 when an additional shield is placed horizontally below the first one, almost closing the aperture of the first one. The magnetic shield has viii Abstract better performance in the vicinity of the shield while the conducting shield has better performance farther away. The magnetic conducting shield is more efficient than the first two in either region. The additional shield improves the shielding factor most significantly in the case of a magnetic conducting type of a shield. ix X Poglavje 1 Uvod Pričujoča disertacija sodi v področji elektromagnetne kompatibilnosti in numeričnih metod v elektromagnetiki, še natančneje, v področji zaslanjanja magnetnega polja in metode robnih elementov. 1.1 Zaslanjanje magnetnega polja Zanimanje za zaslanjanje magnetnega polja presega 100 let [1] in je še vedno aktualno. Pomembnost problematike potrjujeta tudi posebni številki revije IEEE Transactions on Electromagnetic Compatibility na temo elektromagnetne zaščite v letih 1968 in 1988. Danes narekujejo potrebe po zaslanjanju predvsem kriteriji elektromagnetne kompatibilnosti in standardi za neionizirajoča sevanja. Zmanjšanje motečega magnetnega polja v določenem območju se v praksi doseže s prestavljanjem virov polja ali/in z nameščanjem ustreznih kovinskih zaslonov [2, 3, 4, 5]. Kadar na vire polja ni možno vplivati, ostanejo zasloni edini izhod. Učinkovitost zaslonov pogojujejo frekvenca, položaj, oblika, velikost, debelina, elektromagnetne lastnosti in električna povezovanja. Zaradi velikega števila vplivnih parametrov je dobro dizajniranje zaslonov možno le z usrezno metodo, ki omogoča izračun magnetnega polja v njihovi okolici, s tem pa tudi oceno učinkovitosti zaslanjanja. V nadaljevanju so v kratkem opisane značilnosti metod, ki jih zasledimo v literaturi. Najprej so se uporabljale različne analitične ocene učinkovitosti zaslanjanja [1, 6, 7, 8, 9, 10, 11, 12, 13]; te so uporabljive za enostavne geometrije ter postavitve virov in zaslonov. Analitični pristopi kasneje seveda niso zamrli, ampak so se še naprej razvijali in dopolnjevali [14, 15, 2, 16, 17, 18, 19, 20, 21]; njihova prednost je v relativno hitrem izračunu učinkovitosti zaslanjanja. Kot takšni slonijo na določenih privzetkih, poenostavitvah in 1 1.1. Zaslanjanje magnetnega polja omejitvah (npr. okrogli ali ravninski zasloni, idealizirane elektromagnetne lastnosti, nizke frekvence itn.), ki jim v resnici ožajo domeno uporabnosti in mejo natančnosti. Pomanjkljivosti analitičnih pristopov odpravljajo primerne numerične metode. V elek-tromagnetiki nasploh je zelo uveljavljena metoda končnih elementov [22, 23]. Pri uporabi te v problemih zastiranja [24] izstopata težavi odprtih mej in tankih zaslonov, ki posledično terjata razmeroma veliko število elementov, obsežen spomin in dolg računalniški čas [17, 5, 25], zato jo v literaturi zasledimo bolj v vlogi verifikacije kakšnega drugačnega pristopa [17, 20, 26, 5, 27]. Podobne ugotovitve veljajo tudi za metodo končnih diferenc [28]. Težave metode končnih elementov rešijo ali omilijo metode, ki slonijo na konceptu virov. Te spadajo v skupino integralskih metod in so pri problemih zaslanjanja dokaj pogosto uporabljene. Primer takšne je metoda tokovnih niti [29, 30, 31], kije primerna za prevodne, vendar nemagnetne zaslone [3, 26, 32, 33, 34, 27, 35]. Ker se pri njej diskretizira le notranjost zaslonov (ne pa tudi ostalega prostora), je potrebno število neznank precej manjše kot pri metodi končnih elementov. Pri zelo nizkih frekvencah je za magnetne zaslone možno uporabiti tudi magnetostatično aproksimacijo, ki sloni na tokih magnetizacije po površinah zaslonov [36, 31, 32, 33]. Poleg načinov, ki obravnavajo zastiranje polja z dolgimi zasloni, omenimo še nekaj tistih, ki omogočajo analizo 3D problemov. Dela [26, 5, 37] obravnavajo zaslanjanje polja z nemagnetno prevodno ploščo, ki je tako tanka, da se ji sme prirediti žični model. Zastiranje polja s tanko kovinsko steno, na kateri se lahko uporabi impedančni mejni pogoj, obravnava z metodo robnih elementov delo [38], z metodo končnih pa delo [25]. Razen metode končnih elementov, ki ni najbolj primerna za analizo zaslanjanja, je poglavitna pomanjkljivost ostalih pristopov v tem, da so vezani na določene posebne skupine primerov, pri katerih so postavljene kakšne omejitve glede geometrije in/ali snovnih lastnosti zaslonov, medsebojne lege zaslonov in virov polja, frekvence virov in podobno. V praksi se dogaja, da je zaslonov več [2, 17, 39] in da ti ne sodijo nujno v isto skupino, vendar zaenkrat ni zaslediti načina, ki bi bil primeren za vse. Še najbližji temu je kombiniran način [32, str. 31], ki v področju zelo nizkih frekvenc združuje magnetostatično aproksimacijo in koncept tokovnih niti; kot takšen omogoča analizo prevodnih nemagnetnih zaslonov in neprevodnih magnetnih zaslonov. 2 1. Uvod 1.2 Obravnavane zaslonske strukture Problema zastiranja se bomo lotili na splošni 2D strukturi žic, tokovodnikov in zaslonov (slika 1.1). Viri primarnega magnetnega polja so harmonični toki v vzporednih tokovo-dnikih in/ali tankih žicah (žice sprejmimo kot tako tanke, da jih smemo obravnavati kot tokovne premice), njim vzporedni pa so prevodni in/ali permeabilni zasloni, ki primarno magnetno polje delno popačijo. Obravnava tokovodnikov in zaslonov (oziroma valjev ali cilindrov) bo pri predlagani metodi poenotena. Nadalje predpostavimo, da so cilindri iz linearnih snovi, kar pomeni, da bomo analizo elektromagnetnega polja opravili v kompleksnem prostoru, vse poljske količine pa razumeli kot kompleksne. Prečne izmere zaslonske strukture so praviloma neznatne v primerjavi z valovno dolžino in dolžino strukture l0, kar vse govori v prid kvazistatičnega opisa polja v 2D prostoru [40, razdelek 3.3]. £C2 £cnc 7Cnc, MCnc So 7o = 0, /j,(i yk T T Twnv © Iwnw Slika 1.1: Presek splošne strukture zaslonov, tokovodnikov in tankih žic. Pomeni oznak na sliki 1.1 so sledeči: • TiW ~~ StGVllO Z1C, • nc - število cilindrov, • -*W1; -*W2; * * * j ^Wn\v — tOCKG OSI Z1C, • Sei, Sc2, -.., Scnc - preseki cilindrov, • £ci, £c2, • • • > Ccnc - robovi cilindrov, • So - zunanjost valjev, • 7ci, 7C2, • • •, 7cnc in m, m, • • •, ßcnc - specifične prevodnosti in permeabilnosti cilindrov, z x 3 1.3. Metoda robnih elementov • 7o = 0 in ßo - snovni konstanti zunanjosti valjev, • -*W1> -*W2> * * * j ^Wn-yv ~~ t OKI V Z1CB.I1, • ^Cl; ^C2; • • • i ^Cnc ~~ TJOKl V ClllIlQilll, • n eksterna normala cilindra in • t - enotski tangentni vektor cilindra, da velja nxt = ez. Če je T robna točka nekega cilindra, potem oznaki T- in T+ pomenita točki, ki sta njej inflnitezimalno blizu: prva še znotraj in druga že zunaj cilindra. Unijo presekov cilindrov označimo s S = SCI USC2 U • • • USCnc, unijo njih robov pa z C = £C1 U £C2 U • • • U CGnc. 1.2.1 Vrednotenje učinka zastiranja V tem delu je za merljivost učinka zastiranja izbran tako imenovan faktor zastiranja sB, ki je v splošni točki T definiran [27] s kvocientom absolutnih vrednosti vektorjev gostot pretoka primarnega B0 in rezultančnega magnetnega polja B: SB(T) = l^)l. (1.1) Omenjeno razmerje se pogosto podaja v logaritemski skali; v [17] ga imenujejo učinkovitost zaslanjanja. 1.3 Metoda robnih elementov Glede na pristop k problemu moremo numerične metode v elektromagnetiki razdeliti na tiste, ki slonijo na diferencialnem, in tiste, ki slonijo na integralskem pristopu. Najbolj znani predstavnici prvega sta metodi končnih elementov in končnih diferenc, zelo razširjena predstavnica drugega pristopa pa je metoda robnih elementov [41, 42, 30, poglavje 9]. V tem delu je prikazana uporaba metode robnih elementov pri problemih zastiranja magnetnega polja. Prvi korak metode je izpeljava robnih integralskih enačb, temu sledi razdelitev integracijskega območja oziroma robov opazovanega področja na množico kratkih segmentov, imenovanih robni elementi, in zapis aproksimacij iskanih robnih funkcij z interpolacijskimi funkcijami. Z uporabo momentne metode [43] se robno integralsko enačbo nato prevede v sistem algebrskih enačb. Za določitev koeficientov matrične enačbe je potrebno izračunati množico integralov po posameznih elementih. Rešitev matrične enačbe da aproksimacijo 4 1. Uvod iskanih funkcij na robovih opazovanega območja. Iz njih se nadalje izračuna polje kjerkoli v opazovanem območju. V monografiji [44] iz leta 1978 je C. A. Brebbia prvi vpeljal izraz metoda robnih elementov in metodo tudi formaliziral. Od takrat naprej se je začela hitro razvijati in tudi vse bolj uporabljati. Metoda ima v primerjavi z drugimi, ki slonijo na diferencialnem pristopu, kar nekaj prednosti: • dimenzija problema se zmanjša (diskretizirani so le robovi območij), • pisana je na kožo problemov odprtih mej, • natančnost izračuna polja na mejah je ponavadi večja, • ni potrebno računati polja v celotnem prostoru in • napake zaradi aproksimacij se v integralnem smislu delno kompenzirajo, ima pa tudi nekaj pomanjkljivosti: • uporabljane Greenove funkcije so v nekaterih primerih težko določljive, • potrebno je veliko število časovno potratnih numeričnih integracij, • integrali singularnih funkcij terjajo posebno obravnavo, • rezultirajoče matrike so polne in včasih slabo pogojene ter • težje je uporabljiva pri kompleksnejših in nelinearnih problemih. 5 Poglavje 2 Koncept ekvivalentnih ploskovnih virov Poglavje predstavlja implementacijo metode robnih elementov pri problemih zaslanjanja iz razdelka 1.2. V prvih dveh razdelkih je postavljen matematični model, ki opisuje obravnavan elektromagnetni problem, v tretjem je opisan numerični postopek reševanja inte-gralskih enačb oziroma njihovega preoblikovanja v sistem algebrskih enačb, v zadnjem pa je prikazan način izračuna polja v željenem delu prostora in joulskih moči v cilindrih. 2.1 Robni integralski enačbi Začnimo z izpeljavo robnih integralskih enačb za porazdelitve gostot ekvivalentnih virov polja na robovih cilindrov (£). Koordinatni sistem (glej sliko 1.1) orientiramo tako, da je os z vzporedna s cilindri. Električni toki vdolž žic, vodnikov in zaslonov določajo z komponento vektorskega magnetnega potenciala: J = ezJz(x,y) ==> A = ezAz(x,y). V kvazistatičnih razmerah velja za potencial Poissonova enačba [40, razdelek 8.1], {-ß(T)Jz(T), T eS raW _____ , (2.1) - E ßoIwiÖ2D{TTWi) , T G S0 %=\ kjer so ß permeabilnost, /x0 permeabilnost praznega prostora, T(x,y) splošna točka, TT^l = ^(xWi - x)2 + (yWi - y)2 razdalja med točkama T in TWi v (prerezni) ravnini x-y in Č2D Diracova 2D funkcija, za katero velja [42, str. 12] / Ö2v{TT')ds = < , (2.2) Sxy I "i -*■ t1 ^xy Sxy je poljubna ploskev v ravnini x-y, T pa integracijska točka. Tokovno gostoto izrazimo z električno poljsko jakostjo E, to pa z vektorskim magnetnim A ter skalarnim električnim 6 2. Koncept ekvivalentnih ploskovnih virov potencialoma [45], (2.3) kjer so 7 specifična električna prevodnost, j = ^-L imaginarna enota in u krožna frekvenca harmoničnega polja. Zaradi neodvisnosti polja od koordinate z je dV{T)/dz = -U(T)/lQ. U{T) je razlika električnih potencialov med začetkom in koncem cilindra v T G S in je lastna posameznemu cilindru, l0 pa je razdalja med koncema. Za tokovno gostoto dobimo izraz Jz(T) = -jwy(T)Az(T) + l(T)^p. (2.4) Ko enačbo (2.4) zanesemo v Poissonovo enačbo (2.1), dobimo AAZ{T) + k 2{T)Az{T) = -/x(T)7(T)^^, T G S (2.5a) nw AAZ(T) = - J]/xo/wt52D(TIÄv7) , T G So, (2.5b) i=\ kjer je k = V~jw/x7 valovno število cilindra, ki je v vlogi tokovodnika ali zaslona. Reševanja diferencialnih enačb (2.5) se lotimo z metodo robnih elementov [41, 42, 30]. V ta namen uporabimo Greenovi funkciji za notranjost in zunanjost cilindrov1 g{T,T') = \:H(02){kTTTi) (2.7a) go(T,T') = -—lnTT, (2.7b) 2% kjer je H{02) Hankelova funkcija druge vrste ničtega reda. Da bi Greenove funkcije tipa (2.7a) ne pisali posebej za vsak cilinder, smo uporabili kar oznako kT, ki naj pomeni valovno število tistega cilindra, ki mu pripada točka T. Greenovi funkciji v (2.7) sta rešitvi Helmholtzove oziroma Poissonove enačbe za polje linijskega vira ö2d{TT') v točki T' [47, razdelek 1.6]: Ag(T, T') + k2Tg{T, T') = -52D(TT) (2.8a) Ag0(T, T') = -Ö2d(TT> . (2.8b) V primeru zaslonske strukture, ki je vzporedna s površino zemlje, upoštevamo vpliv zemlje z naslednjo Greenovo funkcijo [46], [31, razdelek 30]: gso(T, T ) ~ — In=, (2-6) 2% TT' kjer je de — l,85/|fcg| Carsonova razdalja, ks pa je valovno število zemlje. Kadar de močno presega razdaljo TT', je vpliv zemlje zanemarljiv. Pri nižjih frekvencah je ta pogoj vsekakor izpoljen. 7 2.1. Robni integralski enačbi V drugo Greenovo formulo za območje valjev (S) vstavimo Greenovo funkcijo g in z-komponento Az magnetnega potenciala, za območje zunanjosti valjev (S0) pa funkciji g0 mAz, g(T,T')AAz(T) - Az(T)Ag(T,T') gM°^-MT-)%^äl (2.9a) g0(T,r)AAz(T) - Az(T)Ag0(T,T') = -go(T, T') dA^T+) + AZ(T+) 9go(f' T']čl, (2.9b) c kjer je T integracijska, T' pa splošna točka2. Na levih straneh enačb (2.9) upoštevajmo enačbi (2.5) in (2.8) in dobimo: f (az(T)52d{TV) - p(T)j(T)^pg(T,T')\ (9(T,r0^)_^(T_)2^)c„ (2Jla) 1=1 = -go(T, T') dA^T+) + AZ(T+) 9go(f'T/) čl. (2.11b) c Nadalje upoštevajmo še mejna pogoja za tangencialni komponenti magnetne, Ht = Bt/ß = (-dAz/dn)/[i, in električne poljske jakosti, Ez = -juAz - dV/dz = -juAz + U/h, na robovih valjev, Ht(T+) = HAT.) => - l: ^ =------— z ~ = ai(T) (2.12a) on ßr{T-) on EZ{T+) = EZ{T_) => AZ(T+) = AZ{T_) = a2(T), (2.12b) Greenovi funkciji in njuna normalna odvoda so čez rob L zvezni, lini g(j_,_/ ) = qll,l ) lini qo\l + ,l ) = qo(l,l ) (z.lua) og(l-,l ) og(l,l ) ogo{l + ,l ) ago(l,l J , , lini----------------- =--------------- lini------------------- =-----------------, (z.IObJ t_^t on on t+^t on on kar smo v enačbah tudi upoštevali. _ _ 8 2. Koncept ekvivalentnih ploskovnih virov kjer je /xr relativna permeabilnost, pa dobimo: / (az(T)52D(TT) -^(Th(T)^£(T,T/)) tvQ -ßr(T_)ai(T)g(T, T') - a2(T) ^(J' T'] .;uUiW,i )-ui2K±)----t------dl (2.13a) raw J2^(Twi,T')ß0Iwi+ Az(T)82d{TT>) äs So n/ ai(T)g0(T, TO + a2(T) 9go(f' T']dl. (2.13b) Zatem točko T' v enačbi (2.13a) limitirajmo proti robu C od zunaj (T' -»■ T|), v enačbi (2.13b) pa od znotraj (T' -»■ Ti). Pri tem upoštevajmo rezultate3 iz [47, razdelek 3.1] oziroma lastnosti potenciala enojnega in dvojnega sloja [42, str. 20-21] iz [31, razdelka 11 in 15]. Končno dobimo a2(T') +2 UY(T_)ai{T)g(T,T') + a2{T)dg{TdT,) dl /mt)7(tÄt, T|)ds = 0, T'_£ (2.15a) 2 S a2(T') -2 ^(T^o^TO+aa^^^J'^^d/ c raw = 2^2go(Twi,T')ßoIwi, T'eC. (2.15b) Integrale po presekih valjev (5 = «Sei U 5C2 U • • • U SCrac) preoblikujmo v krivuljne po njih robovih s pomočjo Gaussovega stavka in z upoštevanjem enačbe (2.8a). Za j-ti valj I ßr(T-)ai(T)g(T,T-L.) +ö2(T)------^------- Id/ a2(J ) ,„ ,„ „ „/ rn.og(l,l ) i =---------h fj^ll-)ai(l )g(l, 1 ) + 0.2(1 )------t\--------cii, 1 _ L 2 dn L I ai(T)go(T,T-) + «2(T)------7?------- Id/ a2(-t J J / ,„, ,„ „/ ,_,(»JO U,-l J 1 11 m' =---------------h «i(-t )9o{i , 1 ) + o.2{l )--------------- d/, 1 _ L. 2 @n _ _ A 9 2.1. Robni integralski enačbi dobimo T'+)ds =/xcj7cj^ g{T,T'+)ds 'o 'o =ßGjr)Gjr^ yy ( —^ ' ^9(T,T'+) — Ö2d(TT^) lds 'O Krp ■Sej Ucj-1 i dg(T,T ^ iQ rvr) k k2Cj dn dl =± / ^f^, (2.16) f/cj je razlika električnega potenciala med začetkom in koncem j-tega valja, kCj = -JLü/icjlCj pa njegovo valovno število. Kar velja za en cilinder, velja tudi za presek S z robom C: jlATHT)'mg{T,T>^^j'mM£^E — \ amVm(T), (2.27) ai,m Ctom in sta äitm in a2,m približni vrednosti iskanih gostot enojnega in dvojnega sloja na m-tem robnem elementu. V sistemu integralskih enačb je poleg vektorja a neznana tudi funkcija potencialnih razlik U med konci valjev, ki je (zaradi konstantnih vrednosti na območjih posameznih valjev) na preseku S v resnici že stopničasta funkcija: rac U{T) = ^ UGjfGj{T), Te S, 3 =1 (2.28) 5 ˜ Kadar je rob L poligonalen, se original in aproksimacija ujemata, L = L. _ 13 2.3. Numerično reševanje sistema integralskih enačb kjer so Cj , j = 1,2,...,nc. (2.29) Če v matrični integralski enačbi (2.20) rob C in iskani vektor a nadomestimo z njunima aproksimacijama (2.25) in (2.27) ter upoštevamo (2.28), dobimo približno enačbo: "E p / "-E \ C„2 ' E VmVmiT') + f K(T, T') • ^ am^m(T) d/ m=1 \m=1 C (rac \ J^ UCj - x)2 + (y> - y)2. Določimo najprej normalni odvod te razdalje: VI dP „ „ n,------v>----,------v> {x-x',y-y') n P n • v P = n • V v [X1 — x)z + (y' — y)z = n _ -------------------------------- _--- j(x' - x)2 + (y' - y)2 = n — dn y/{x - x')2 + {y- y')2 P (2.36) kjer je P = TT' = (x' - x,y' - y) distančni vektor od točke T k T'. Analitična izraza za normalna odvoda Greenovih funkcij sta dgo(T,T') , 1 1 „ n-P -------------= n • Vcfnl-/,-/ =-----ii'V tar =------ n • V P = —— (z.cwa) dn y v ' ; 2ti 2%P 2%P2 če je potencialna razlika med koncema kakšnega valja znana ali so te za več električno povezanih valjev enake, se število neznank temu ustrezno zmanjša. 15 2.3. Numerično reševanje sistema integralskih enačb ME^Il = n . Vg(T,T') = -n • VH^\kTP) = --H[2\kTP)n ■ V(kTP) w HY\kTP), pri slednjem smo upoštevali [53, str. 463], daje dH° ^ =-H[2\u). (2.38) Pri določanju koeficientov sistema (2.34) je največ dela z računanjem integralov elementov matrik K (2.19a) in K^ (2.19c) vzdolž posameznih segmentov. Te integrale računamo numerično; uporabljamo adaptivni rekurzivni Simpsonov kvadraturni algoritem [54], kije vgrajen v programski paket MATLAB® [55, str. 2-2694]. Posebno obravnavo terjajo tako imenovani lastni prispevki, ko je m = p. Pri njih so integrandi, ki vsebujejo Greenovi funkciji in njuna normalna odvoda, singularni v središčnih točkah segmentov, zaradi česar numerična integracija odpove. Te singularnosti so integrabilne in o njih bo tekla beseda v nadaljevanju. 2.3.1 Lastni prispevki Lastni prispevki v sistemu enačb (2.34) so (glej enačbi (2.19a) in (2.19c)): Kn(T,Tp)dl = 2 ßr(T_)g(T,Tp)dl K12(T,Tp)dl = 2 dg(T,Tp) dn K2l{T,Tp)dl = -2 g0(T,Tp)dl K22{T,Tp)dl = -2 9go{\,Tp\\l (2 3Q) SjCp S.£p 5£p 5£p |K[,l(T,Tp+)di.||M|Zk>d, Pri prvem je v integrandu razen Greenove funkcije ali njenega odvoda še relativna perme-abilnost /xr, ki je tudi, podobno kot potencialna razlika (2.28), stopničasta funkcija: rac ßr(T) = V ^c,^c,(T_) g(T,Tp)dl ßo J j (2.41) f Kn(T,Tp)dl = — [ iT ßCicXcy)P = *p, (2-45) i=i kjer je fcp okrajšava za valovno število cilindra, ki mu pripada p-ti robni element. Integral Greenove funkcije (2.7a) po bližnjem intervalu 8% je (glej izpeljavo (2.43)) 5%/2 5%/2 /5(T,Tp)d/ = ihm / ^(^de-if^-^hm /in^d^) 2i £^o 2i 2 ti £^o 2 8 In I 2, ^Yekr>8 lr> \ kjer je 8% dolžina bližnjega intervala. Bližnji interval 8% smo izbrali tako, da največja absolutna vrednost argumenta Hankelove funkcije na njem ne presega 1/10, kar upravičuje uporabo aproksimacije te funkcije pri majhnih argumentih. Integral v preostalem delu segmenta S£p izračunamo numerično. Lotimo se sedaj četrtega integrala v (2.42). Njegov integrand je enak nič (glej (2.37a)), dg0(T,Tp) n Pp n P --------------- = --------- = ------- dn 2nP£ 2n(2 ------- ------- U, [zAlj saj sta n in Pp pravokotna na integracijskem intervalu e < £ < 8/p/2, e -»■ 0 (glej sliko 2.2), zato je Sžp/2 /W^ih /^ = 0. (248) 8£p e Enako velja tudi za drugi integral v (2.42), 8£p kajti tudi njegov integrand vsebuje skalami produkt n • Pp (glej (2.37b)). Preostal je še zadnji integral v (2.42) (glej [47, razdelek 3.1.2]). Tudi pri njem razdelimo integracijsko območje na bližnji interval 8% točke Tp, kjer je upravičena apro-ksimacija (2.44) Hankelove funkcije pri majhnih argumentih, in preostali del S£p \ S2£p. 18 2. Koncept ekvivalentnih ploskovnih virov Pri integralu na preostalem delu integrand ni singularen. Ker je točka Tp+ infinitezimalno blizu točke T, velja f dg(T,Tp+) f dg(T,Tp)__ dn dl= fo*1- (2-50) 8£p\82£p 8£p\82£p Pri izračunu integrala na bližnjem intervalu 8% upoštevajmo, da je normalni odvod d/dn = d/dr] (glej sliko 2.2), in dobimo f dg(T,Tp+).^ 1 f d\n(TTp+) dn d/ - ~2n dn d/ 82£p 82£p 2 - lim lim / lim------v w lp+J-----d£ (2.51) S%/2 = - lim lim / ^ d£ = - lim arctan^ = -. Velja tudi / 9g(T,Tp+) f 3<7(r,rp) 1 On dn + 2' l j saj je desni integral enak nič (glej enačbo (2.49)). Če združimo rezultate (2.50) in (2.52) ter upoštevamo še (2.49), dobimo f dg(T,Tp+)A] . fdg(T,Tp) 1 1 9n d/ = ^n dl + 2 = 2 • (2-53) 2.4 Izračun magnetnega polja Ko enkrat rešimo sistem sklopljenih robnih integralskih enačb in izračunamo ekvivalentne ploskovne vire ter potencialne razlike, se lahko lotimo izračuna polja kjerkoli v prostoru. 19 2.4. Izračun magnetnega polja 2.4.1 Izračun vektorskega magnetnega potenciala Izraza za vektorski magnetni potencial znotraj (T' e S) ali izven valjev (T' e SQ) sledita iz enačb (2.13a) in (2.13b): (2.54a) 4,(T' eS) = -f (a2(T)d9^T>) + /xr(T_)ai(T)g(T,T')) dž 4, (T' eSQ)=f (a2{T)d9o^T,) + ai(r)5o(T,T7))d/ + X>(TW*,T'^Iwi- (2.54b) Ploskovni integral v zgornji enačbi preoblikujemo v krivuljnega, podobno kot smo to storili v razdelku 2.1, in upoštevamo, da je T1 e S. Na preseku npr. j-tega valja dobimo /Verhörter, T')ds =ßcj7cjr^ / g(T,T')ds '0 "O 'o 'o 5Cj 5Cj =^,■70,— i (-V • V5(r,r') - 52D(TT7))ds 5Cj ^C7 =ßcj7cj^ 2 ggn dž + / ^d(TT') ds (2.55) za celoten presek S valjev pa velja ß{THTms{T, T')ds=^-[ i^iZl^(J'T)d^+ f E^£lß2D(TT') ds ] 1 U{T)dg(T,T') 1 [/(T') dg(T,T') d, wk 9^d/+]^ + —V. (2.56) 20 2. Koncept ekvivalentnih ploskovnih virov Če to upoštevamo v enačbi (2.54), dobimo8 AZ(T G S) c + /xr(T_)a {T) g {T,T')) dl + 1 U(T) dg(T,T') 1 U(T') (2.57a) c ur €«,) = (a,m^ dn + a1(T)g0(T,r)dl + Ytg0(TWi,r)(i0IWi. (2.57b) Ti enačbi lahko jedrnato zapišemo v matrični obliki podobno, kot smo to že počeli v razdelku 2.1.1. V ta namen vpeljimo naslednje okrajšave: K(T, T') -liT(T_)g(T,T) -dg(T,T')/dn go(T,T>) dg0(T,T>)/dn Ku(T,T') 1 dg{T,T) Wo 1 [1 1 ^ W °J dn 0 ^ Gw(T') Matrična oblika enačb (2.57) je raw 0 E 9o(TWi,T')iioIwi ^ (2.58a) (2.58b) (2.58c) (2.58d) AZ(V g 5) 4* (T' g So) K(T, T') • u(T)dl + K^T, T')f/(T)d/ + CuU(T') + G W (T'). (2.59) Če krivuljne integrale diskretiziramo na podoben način, kot smo to storili v razdelku 2.3, dobimo sledeče: AZ(T g Š) AZ{T' g So) ™E / f ( nc \ \ Y^ / K(T, T')d/ • am + ^ UGj(cj,m / Et, (T, T')d/ ~~ 8£m •?— 8£m rac + Qu J2 Ucj ^ ( f dK(T,T>)A1 - f^ \ f dKuJT,!*)^ \ /_j I f) / —m * / j UCjSCj,m f) f I ' v-^-ö4J kjer je spremenljivka w' ali x' ali y'. Odvod predzadnjega člena v enačbi (2.60) je enak nič, saj so potencialne razlike (2.28) konstante na območju valjev. Integrale v enačbi (2.64) izračunamo numerično tudi tokrat. 22 2. Koncept ekvivalentnih ploskovnih virov Odvodi Greenovih funkcij V izrazu (2.64) nastopajo odvodi elementov matrik K, K^ in Gw, v njih (glej enačbo (2.58)) pa Greenove funkcije (2.7); določiti je torej potrebno parcialne odvode prvega oziroma drugega reda teh funkcij: dg(T,T') ------------ dw' dg0(T,T') dw' 1 dH0 (2), kTTT>) dw' 1 dlnTT ~2*lh^ rßg {T, T1) dw'dn d2go(T,T') -------------------------------- zzz dw'dn 1 d2H{2\kTTT') 4j dw'dn 1 d2\nTT 2% dw'dn ' (2.65) kjer je w' ali x' ali y'. Argument funkcij je razdalja P = TT' = ^{x' - x)2 + (y' - y)2 med točkama T(x,y) in T'(x',y'). Določimo najprej njen odvod po koordinati w': dP d w' 9w, WVMHW p w (2.66) kjer je w = x, če je w' = x', inw = y, če je w' = y'. Tretji odvod v (2.65) je dg0(T,T') _ 1 d In P _ 1 dP dw' ~2*ThT ~2^Pdw' Nadaljujmo s prvim (glej enačbo (2.38)): dg{T,T') 1 dH{2\kTP) 1 u{2)MkTP) w'-w (2.67) ajHi {kTP) dw' dw' 4j dw' Lotimo se odvoda normalnega odvoda (2.37a): 1 / 1 9(n • P) d(l/P2 - kT%pW) H[2\kTP). (2.68) d2g0(T,T') dw'dn ) 1 1 d{nx{x' - x) + ny{y' - y)) n ■ P dP 2% yP2 dw' 2 P3 dw' U^-2{w'-w)^\=^(r^-{w'-w)X^f] 2% P2 v P4 %P2 2 v P2 (2.69) kjer je P = TT' = (x' -x,y' - y) distančni vektor od točke T k T'. Za konec ostaja še odvod normalnega odvoda (2.37b): d2g(T,T') kT (d{-n-V)H\2)(kTP) d(l/P) (2) n-PM,(2)feP) ^šdn = Tj\dw~'----------P------+ n'P^^^ (kTP) + —--------^ kT (nwH[2\kTP) n-PdP ,2)n = ^ --------P-------------P^d^ž^ {kTP) + ^(^>-^iW))^ kTP 1 dw' kTP) + nw p^ kT(w' - w)X^H{2){kTP) + Lw-2(w'-w)^H[2\kTP) (2.70) _ _ _ _ _ _ 23 2.4. Izračun magnetnega polja kjer smo upoštevali [53, str. 463], daje (2) dH\ ^ = H^\u) - -H\2\u). (2.71) du u 2.4.3 Izračun joulske moči v cilindru Pri računanju joulske oziroma delovne moči v posameznem cilindru izhajamo iz Poyntin-govega stavka v kompleksnem prostoru [45, razdelek 2.20]. Iskano moč na dolžinsko enoto v j-tem cilindru določa integracija Poyntingovega vektorja po njegovem obodu, p,,c3 = - Re[Sn]dl, (2.72) kjer je Re[Sn] realni del normalne komponente Poyntingovega vektorja9 S = E x H*. Normalno komponento Sn = -EZH? na robu CCj moremo po enačbi (2.12) izraziti z ekvivalentnimi viri in potencialno razliko UCj na tale način: dz ßcj on 'o ^o --------------------------Ctl h ßo ßo Če upoštevamo še Amperov zakon (podobno kot v enačbi (2.21)), dobimo Vi a = Re fe / ^dl - i— i a\a2dl = ,Cj Gj + — i Im[aja2]d/. (2.74) lo ßo ßo h ßo Z upoštevanjem stopničaste aproksimacije (2.27) in koeficientov pripadnosti (2.31) se integral prevede v vsoto, Re[[/Cj/cj] nE V i,G j ~-------------- H-----2^ Ccj,m5/mlm[ä1)mä2,m], (2.75) l0 ßo m=\ tok ICj v prvem sumandu pa izrazimo z enačbo (2.35c). Kadar se kompleksorji nanašajo na amplitude harmoničnih količin, je S = (E x H*)/2. 24 Poglavje 3 Verifikacija predlagane metode Rezultate predlagane metode smo primerjali z: • rezultati analitično rešljivih primerov, • rezultati meritev in • rezultati dveh drugih numeričnih metod. Iz primerjav sta razvidna natančnost in obseg uporabnosti predlagane metode. 3.1 Cevasti zaslon kroˇznega preseka 3.1.1 Vodnik v koaksialnem zaslonu Na sliki 3.1 je prikazan primer žice znotraj zaslona, ki geometrijsko ustreza koaksialnemu kablu. Vir primarnega polja je harmoničen tok I0 v žici oziroma žili. Oklop ima vlogo zaslona ter ima odprta konca. Rezultančno magnetno polje je v tem primeru analitično določljivo. Za funkcijo radialne porazdelitve inducirane tokovne gostote J = ezJz v zaslonu velja Besselova diferencialna enačba ničtega reda [31, razdelek 28] AJz(p) + k 2cJz(p) = ^ + -^ + k 2cJz(p) = 0, pi p2- (3.4) Če v enačbo (3.3) vstavimo nastavek (3.2) za tokovno gostoto in izračunamo njen odvod, dobimo nastavek za magnetno poljsko jakost oziroma komponento Hv v zaslonu: H^p) = (CiJi(fccp) + C2Ni{kGp)) /he, Piz potenciala primarnega polja razvijemo v Fourierevo vrsto, /J0/o^ 1 A0,z(p, P2> P < P\ (3.8a) Pl < P < P2 (3.8b) (3.8c) ra=l kjer so tf^Li in H^_x Hankelovi funkciji prve oziroma druge vrste reda (2n - 1) ter kc valovno število zaslona. Konstante Cx,n, C2,n, C3,n in C4,„, n = 1,2,..., se določi iz mejnih pogojev (2.12) na površinah zaslona. Vrste v enačbi (3.8) konvergirajo zelo hitro (razen v bližnji okolici žic, ko je p « 2d ali p < 2d), zato je dovolj upoštevati le prvih nekaj členov. Pri izračunih smo upoštevali prvih deset členov. Primerjava radialnih odvisnosti (vzdolž polpremice tp = ti/4) absolutne vrednosti magnetne poljske jakosti, izračunanih po metodi robnih elementov (MRE) in analitično, je prikazana na sliki 3.6. Podatki za ta primer so: /0 = 100 A, / = 50Hz, 2d = 30 mm, 29 3.1. Cevasti zaslon krožnega preseka pi = 50mm, p2 = 55 mm, 7c = 10MS/m in pc = Wp0. Omembe vredna odstopanja je opaziti le ob površinah zaslona. 300 250 200 150 100 50 • analit. MRE « \ » v ^nJ 45 50 55 p j mm 60 (a) Magnetna poljska jakost 10 7,5 5 2,5 0 -2,5 -5 -7,5 m m 45 50 55 6 p j mm (b) Relativna napaka Slika 3.6: Primerjava numerične in analitične rešitve za primer dvovoda v cevastem zaslonu. Diskrepanco med numerično in analitično izračunano absolutno vrednostjo magnetnega polja v odvisnosti od razmerja 81/5C smo analizirali na enakem zaslonu in pri enaki frekvenci kot v prejšnjem primeru (glej razdelek 3.1.1). Primarni vir polja je tokrat tok v dvovodu medosne razdalje 2d = 10 mm. Na sliki 3.7 so prikazani rezultati te analize v točki (p = 15mm,^ = ti/4). Absolutna vrednost relativnega odstopanja narašča spet približno potenčno z razmerjem 5Z/<5C, pri 81/5C ~ 1 je okrog 1 %. Hoburg in ostali [60] so merili magnetno polje dvovoda v cevastem zaslonu krožnega preseka. Primerjava vrednosti faktorjev zastiranja, ki so določene analitično, s predlagano metodo in z meritvijo, je za aluminijast zaslon prikazana na sliki 3.8a, za jeklen zaslon pa na sliki 3.8b. Številski podatki so /0 = 500 A, / = 60Hz, 2d = 6cm, p1)A1 = 25,435 cm, p1)Fe = 25,43 cm, p2 = 25,5 cm, 7c,ai = 30,5 MS/m, 7C)Fe = 6,76 MS/m in pc,Fe = 190/xo2. Ostali podatki meritev so podrobno opisani v [60]. Vrednosti faktorjev zastiranja so izračunane na oddaljenosti p = 0,62 m in kotih 0 <

yc+t)= f /(e)e-^cos(^)de. (3.9) o Funkcija /(£) sledi iz mejnih pogojev (2.12) na površinah zaslona, 2 _________ , _________ _________ /(£) = oßCßo \k2 - k2c e*? ( 2/iGßoS\k2 - k% coshft J(2 - k2^) + UU2 + ßl (C2 - fcc)) sinhU^£2 -k%\\ . (3.10) Z odvajanjem potenciala (3.9) pridobimo izraze za komponenti vektorja gostote magnetnega pretoka nad zaslonom (glej enačbo (2.62)): oo Bx(x,y> yc+t) = dMx,y>yc+t) = _ f ^me -(y cos^x)d^ (311a) o oo By(x,y> yc+t) = _dA^,y>Vc+t) = f ^me -ivsin^x)d^ (3.llb) o Integrale (3.9) in (3.11) računamo numerično, neskončno zgornjo mejo pa nadomestimo z ustrezno končno (v primerih, ki sta opisana v nadaljevanju, smo za zgornjo integracijsko mejo izbrali vrednost 104). Absolutne vrednosti gostote magnetnega pretoka, izračunane po metodi robnih elementov in z analitično aproksimacijo, smo primerjali v nekaj točkah nad zaslonom (na 32 3. Verifikacija predlagane metode višini yR = 15 cm nad žico). Za nemagnetni zaslon (7c = 30,5MS/m in /xc = ßo) je primerjava prikazana na sliki 3.10a, za magnetnega (7c = 10 MS/m in /xc = 10/x0) pa na sliki 3.10b. Številski podatki so: /0 = 100 A, / = 50 Hz, t = 2 mm, yc = 5 cm in l = 4 m. Ob sredini zaslona je ujemanje solidno, proti robovom zaslona (x = ±1/2 = ±2 m) pa odstopanje pričakovano narašča. 80 60 40 20 aildllL. ■ MRE / \ / \ • • / \ • • •••• ..•J, r ^ «... • ••• -1 0 1 x j m 120 100 80 60 40 20 0 a malit. lR A / \ / \ ••.. •••5 > • a K. ,••• ..•* -1 0 x j m (a) Nemagnetni zaslon (b) Magnetni zaslon Slika 3.10: Primerjava numeričnih rezultatov z analitično aproksimacijo v primeru zelo širokega ravnega zaslona. 3.3 Raven zaslon Rezultate predlagane metode smo primerjali z rezultati nekaterih drugih metod in meritev tudi na primeru ravnega zaslona nad ˇzicami (slika 3.11). l/2 l/2 it •yc, ßc yc Iwi IW2 IW3 -©-^-----Ö^1-----©- yR M+ _ ^El = E}^y /w. inTT^. (3.13) ^f(T) 2% Iq 2% *■—' s %=l Tudi potencialne razlike U med konci posameznih valjev so v splošnem neznanke, podobno kot v razdelku 2.2, zato je potrebno tej pridružiti še enačbe tipa (2.21), (2.22), (2.23) in/ali (2.24), v katerih izrazimo tok z integralom tokovne gostote po preseku cilindra, Icj= jJz(T')ds', j = 1,2,...,ne. (3.14) SCj Dobimo sistem sklopljenih integralskih enačb. Tudi tega smo numerično reševali z mo-mentno metodo, le da je bilo tokrat potrebno diskretizirati površine presekov cilindrov. Frix in Karady [3] sta merila magnetno polje dvovoda, ki je zaslanjano z ravnim aluminijastim zaslonom. Struktura, ki sta jo uporabila, ustreza primeru na sliki 3.11, pri kateri je srednja žica odstranjena in sta toka v krajnih dveh nasprotna, /Wi = -iw3 = ^o-Uporabila sta dva različna aluminijasta zaslona. V prvem primeru so številski podatki: /0 = 100 A, / = 60 Hz, yc = 6,35 cm, 2d = 7,62 cm, l = 30,48 cm, t = 3,175 mm, 34 3. Verifikacija predlagane metode 7c = (103/36,99)MS/m in /xc = ßo- Zaslon ima odprta konca in gostota magnetnega pretoka je merjena na višini yR = 14,13 cm nad žicama dvovoda. V drugem primeru so podatki spremenjeni pri debelini zaslona, t = 6,35 mm, specifični prevodnosti, 7c = (103/45)MS/m, in višini merilnih točk, yR = 14,45 cm. Preostali podatki meritev so podrobno opisani v [3]. Navajamo primerjavo rezultatov meritev, rezultatov metode robnih elementov (MRE), rezultatov metode tokovnih niti (MTN) in rezultatov hibridne metode (HM). Iz slike 3.12 je razvidno, da je ujemanje rezultatov MRE in MTN izvrstno (poprečje absolutnih vrednosti relativnega odstopanja je |er|sr. = 0,03%) in daje ujemanje z rezultati eksperimenta pri teh dveh metodah (|er|sr. = 5,1%) veliko boljše kot pri HM (|er|sr. = 16,2%). 45 40 35 25 20 15 10 5 k, X meritev N k 0 MRE \ \ HM ,\ «v A s< ^ sk ^ k < X fr- ^g 0 5 10 15 20 25 30 35 x j cm (a) Tanjši zaslon (t = 3,175 mm) 30 25 20 15 10 5 x *> s MRE "> \ 0 MTN N ^ \ k g tK. S-o 0 5 10 15 20 25 30 35 x j cm (b) Debelejši zaslon (t = 6,35 mm) Slika 3.12: Primerjava rezultatov meritev, metode robnih elementov, metode tokovnih niti in hibridne metode pri ravnem aluminijastem zaslonu. 3.3.3 Metoda ploskovnih tokov (MPT) Metoda ploskovnih tokov sloni na magnetostatiˇcni aproksimaciji in je primerna le za analizo magnetnih cilindrov. Na povrˇsinah magnetnih cilindrov je gostota ploskovnih tokov magnetizacije K sorazmerna tangencialni komponenti (n × B) gostote magnetnega pretoka [36], [31, razdelek 22]: K(T) + —ß(T)n(T) x B(T) 0, (3.15) kjer je ß(T) = (/x(T_) - /x0) /(KT-) + ^o)- Če izrazimo vektor B z vsoto primarnega polja (B0) in prispevka gostote K = ezKz ploskovnih tokov, dobimo integralsko enačbo _ 35 3.3. Raven zaslon za porazdelitev le-te: v ' li ' P2 ez ■ (n(T) x B0(T)), (3.16) kjer je P = TT' = (x' -x,y'- y) distančni vektor med opazovano točko T in integracijsko T1. Rezultate treh metod, metode robnih elementov, metode ploskovnih tokov in hibridne metode, smo podali za primer ravnega magnetnega zaslona, ki je idealno dvostransko oze-mljen in se nahaja nad žicami trifaznega voda (slika 3.11): /Wi = /0exp(j27t/3), JW2 = h in /w3 = ioexp(-j27i/3). Na sliki 3.13 so prikazane izračunane vrednosti gostote magnetnega pretoka na višini yR = 40 cm nad žicami. Številski podatki primera na sliki 3.13a so: J0 = 400 A, / = 50Hz, yG = 10cm, d = 8cm, l = 40cm, t = 2mm, jG = lOMS/m in fic = 1000/io- Ujemanje rezultatov HM z rezultati MRE ni prav dobro (|er|sr. = 9,5%). Ujemanje rezultatov MPT z rezultati MRE je še slabše (|er|sr. = 13,4%), saj magneto-statična aproksimacija ne upošteva induciranih vrtinčnih tokov v cilindrih. Opravili smo še izračun (glej sliko 3.13b) pri / = 1Hz. Tokrat je ujemanje z rezultati MRE pri MPT izvrstno (|er|sr. = 0,09%) in tudi boljše kot pri HM (|er|sr. = 1,1%). Magnetostatična aproksimacija je očitno uporabljiva le, če so inducirani vrtinčni toki zanemarljivi, ti pa so pri zelo nizkih frekvencah [3] in pri lameliranih ali feritnih cilindrih. 55 50 45 40 35 0 x j cm (a) / = 50 Hz 40 55 50 45 40 35 30 MR F. < o MPT ? x HM -40 -20 0 x / cm 20 (b) / = 1 Hz 40 Slika 3.13: Primerjava rezultatov predlagane metode z rezultati metode ploskovnih tokov in hibridne metode v primeru ravnega magnetnega zaslona. _ 36 3. Verifikacija predlagane metode 3.3.4 Meritev Na Elektroinštitutu Milan Vidmar smo opravili meritve na primeru zaslona, ki je prikazan na sliki 3.11. Zaslon je bil širok l = lm in dolg 8 m, dvostransko ozemljen in postavljen na višino yc = 20cm nad žicami trifaznega voda: /Wi = /0exp(j27t/3), /W2 = h, iw3 = /0exp(-j27i/3), /o = 100 A, / = 50Hz in d = 8cm. Dolžina trifaznega voda je bila 10 m. Eksperiment je bil opravljen za dva zaslona: za aluminijast zaslon (7c = 30,5 MS/m, ßc = Ho) debeline i = 4mm in jeklen zaslon (7c = 10MS/m, /xc = 100/x0) debeline t = 1,6mm. Zaslon je bil dvoslojen in iz plošč dimenzij Im x 2m. Gostota magnetnega pretoka je bila merjena z Wandel k Goltermann EFA-3 EM poljskim analizatorjem. Ostali podatki eksperimenta so opisani v [32, razdelek 3.3.3]. Vrednosti gostote magnetnega pretoka, ki so bile izmerjene na višinah yR = 30 cm in yR = 70 cm nad žicami, smo primerjali z izračunanimi po metodi robnih elementov. Rezultati primerjave za nemagnetni (aluminijast) zaslon so prikazani na sliki 3.14, za magnetni (jeklen) pa na sliki 3.15. Poprečje absolutnih vrednosti relativnega odstopanja izmerjenih vrednosti od izračunanih je na sliki 3.14a |er|sr. = 5,5%, na sliki 3.14b |er|sr. = 9,6%, na sliki 3.15a |er|sr. = 12,0% in na sliki 3.15b |er|sr. = 4,7%. Pri numeričnem izračunu je dvostranska ozemljitev zaslona privzeta kot idealna (glej razdelek 2.2). Ker sta potekala zaslon in trifazni vod vzporedno s površino zemlje, smo pri numeričnem izračunu upoštevali Greenovo funkcijo (2.6) in za specifično prevodnost zemlje smo privzeli vrednost (1/300) S/m [31, str. 291]. Izračun smo opravili tudi z Greenovo funkcijo (2.7b), vendar se rezultati niso bistveno spremenili. 37 3.3. Raven zaslon 6 4 2 0 MRE o meritev v\ 1 -0,5 0 0,5 x j m (a) yR = 30cm 1,2 1,1 1 0,9 0,8 0,7 0,6 C o ' 0 o o \ o o / \ ( ' v 0 0 O / \ o \ ° / C meritev \ | | 1 -0,5 (b) yR = 70cm 0 0,5 x j m Slika 3.14: Primerjava rezultatov predlagane metode z izmerjenimi vrednostimi v primeru ravnega nemagnetnega zaslona. 20 15 10 5 0 / >, z ,\ r °\ . / \ • «^ V4 S n^ / r 1 ^ •% — MRL 7, i«>^ C meritev ^•J 1 -0,5 (a) yR = 30cm 0 0,5 x j m 4 3,5 3 2,5 2 1,5'^ 1 -1 / '0 '°\ \ ( / \ > / \ ,/ \c C — MRJ meri tev -0,5 (b) yR = 70cm 0 0,5 x j m Slika 3.15: Primerjava rezultatov predlagane metode z izmerjenimi vrednostimi v primeru ravnega magnetnega zaslona. 8 1 1 1 1 38 Poglavje 4 Analiza zaslonskih struktur Predlagana metoda robnih elementov omogoča izračun magnetnega polja pred, v in za zasloni. Z ustreznim prikazom rezultatov tega izračuna pridobimo soliden vpogled v učinke zastiranja, kar je nepogrešljivo pri dizajniranju zaslonskih struktur. Za demonstracijo takšnega prikaza rezultatov smo izbrali splošno konfiguracijo, ki je prikazana na sliki 4.1. Nad snopom vodnikov je zaslon profila U širine l = 30 cm in višine h = 20 cm, ki naj je idealno dvostransko ozemljem Da bi prikazani rezultati ponazorili dva različna principa zastiranja magnetnega polja, ki slonita na indukciji oziroma magnetizaciji, in učinek njunega kombiniranega vpliva, smo izbrali naslednje tri zaslone: prevodnega (nemagne-tnega), 7c = 10MS/m in fic = /x0, (slabo prevodnega) magnetnega, 7c = 0,01 MS/m in ßc = 100/xo, in prevodnega magnetnega, 7c = 10 MS/m in /xc = 100/x0- Preostali številski podatki so: / = 50Hz, d = 6cm, yG = 7cm in t = lern. i i.y l/2 l/2 A yC Iwi IW2 IW3 x ' d ' ' d ' ' 7( h ßc > t Slika 4.1: Zaslon profila U nad snopom žic. h 39 4.1. Enofazno vzbujanje 4.1 Enofazno vzbujanje Na sliki 4.2 so prikazane gostotnice magnetnega pretoka v trenutku 0 s za primer le srednjega tokovodnika, IW2 = 100 A. 30 30 -15 -30 -15 0 15 30 x / cm (a) Prevodni zaslon. -15 -30 -15 0 15 30 x / cm (b) Magnetni zaslon. 30 -15 -30 -15 0 15 30 x / cm (c) Prevodni magnetni zaslon. Slika 4.2: Sive linije prikazujejo gostotnice primarnega polja, krepke pa gostotnice rezul-tanˇcnega polja. . V prevodnem nemagnetnem zaslonu se inducira tok IC = -(95,54+j10,39)A, ki je skoraj nasproten toku v ˇzici, IC ≈ -IW2 = -100A. Ta tok kompenzira veˇcji del primarnega polja za zaslonom, kot je razvidno na sliki 4.2a. V slabo prevodnem magnetnem zaslonu . se inducira tok IC = -(0,016 + j1,227)A, ki kot tak ne zmore kompenzirati primarnega polja. Ker pa se ta zaslon v polju tudi magnetizira, ”ujame“ precej magnetnega pretoka 40 4. Analiza zaslonskih struktur (‘flux shunting’) [61, razdelek B.1] [49, razdelka 3.1 in 3.5], kot je razvidno na sliki 4.2b. . V prevodnem magnetnem zaslonu se inducira precejˇsen tok, IC = -(84,71 + j9,31) A, zmanjˇsanje gostote pretoka za zaslonom pa je zaradi magnetizacije ˇse izrazitejˇse kot v obeh prejˇsnjih primerih (slika 4.2c). 4.2 Trifazno vzbujanje Slika 4.3 prikazuje porazdelitev faktorja zaslanjanja okrog zaslona nad trifaznim vodom: IW1 = I0 exp(j2π/3), IW2 = I0 in IW3 = I0 exp(-j2π/3), kjer je I0 = 100A. Inducirani 80 40 -40 -80 7 ^ r^\ \ \ s |fnfr\\ \ \ V <$ aftSrjN '2-- -a,* -80 -40 0 40 x/ cm (a) Prevodni zaslon. 80 40 0 -40 80 80 40 -40 -80 / .- 3 ~^ \ *2s VIC \ \ iw Vk»&W/^rP ^"^/yv J) x ^ 80 -40 0 40 x j cm (b) Magnetni zaslon. / / /-"—J3\. \ xl ®®® \AA ( v y \\ \ / \ 3 / \ X / \ -80 -80 -40 0 40 80 x/ cm (c) Prevodni magnetni zaslon. 80 Slika 4.3: Porazdelitev faktorja zastiranja okrog zaslona profila U pri trifaznem vzbujanju. 0 0 41 4.3. Trifazno vzbujanje in dvozaslonska struktura . tok v zaslonu je tokrat zanemarljiv: v prevodnem nemagnetnem zaslonu je IC = (0,761 + . j0,324)A, v slabo prevodnem magnetnem je IC = (0,000 + j0,015)A in v prevodnem . magnetnem zaslonu je IC = (0,973 + j0,008) A, kar je zelo blizu primera, ko bi imel zaslon odprta konca, IC = 0. Cˇe primerjamo porazdelitev faktorja zastiranja za prevodni nemagnetni zaslon (slika 4.3a) in slabo prevodni magnetni zaslon (slika 4.3b), opazimo, da je tik za zaslonom boljˇsi uˇcinek magnetnega, bolj stran pa uˇcinek prevodnega zaslona. V primeru prevodnega magnetnega zaslona, kot je razvidno na sliki 4.3c, se faktor zastiranja izdatno poveˇca povsod za zaslonom. Priˇcakovano najslabˇse zastiranje je pod odprtino. 4.3 Trifazno vzbujanje in dvozaslonska struktura Da bi izboljˇsali zastiranje tudi na spodnji strani, smo dodali ˇse raven zaslon, ki trifazni vod ˇˇ skoraj zapre. Sirina tega zaslona je 32cm, vse ostalo pa je isto kot pri zaslonu U. Sirina reˇze med zaslonoma je 2 cm. Porazdelitev faktorja zastiranja te dvozaslonske strukture je prikazana na sliki 4.4. Pri prevodnih nemagnetnih zaslonih (slika 4.4a) in slabo prevodnih magnetnih (slika 4.4b) se je zastiranje le delno izboljˇsalo, pri prevodnih magnetnih zaslonih (slika 4.4c) pa se je faktor zastiranja obˇcutno poveˇcal povsod. Zaslanjanje je ˇse najslabˇse ob reˇzi med zaslonoma, kar je tudi priˇcakovano. 4.4 Analiza povratnega uˇcinka Kot primarni vir polja smo izbrali vodnika dvovoda polmera 2 cm s tokoma ±100 A. Osi tokovodnikov sovpadata z lego prve in tretje ˇzice na sliki 4.1 (srednje ˇzice tokrat ni). Na sliki 4.5 so prikazane gostotnice magnetnega polja v trenutku 0 s (sive so gostotnice primarnega polja). Pri prevodnem in magnetnem zaslonu sta (ponovno) razvidna vpliva indukcije in magnetizacije: pri prevodnem se gostotnice “stisnejo” k vodnikoma dvovoda (slika 4.5a), magnetni zaslon pa “vsrka” veˇcino magnetnega pretoka (slika 4.5b). Pri kombiniranem zaslonu sta omenjena uˇcinka “prepletena” (slika 4.5c). Zaradi uravnoteˇzenosti primarnih tokov (±100 A) je induciran tok v kateremkoli zaslonu sicer zanemarljiv, ne pa tudi njegova gostota, zaradi katere se v cilindrih sproˇsˇca toplota. Zaslanjanje oˇcitno ni “zastonj” in ga kaˇze osvetliti ˇse s te plati. Joulska moˇc na dolˇzinskem metru j-tega cilindra je doloˇcena z enaˇcbo (2.75). V primeru prevodnega zaslona se izgubna moˇc v vodnikih dvovoda poveˇca za 2,63 %. Ker pa je celotno poveˇcanje izgub zaradi zastiranja 37,02 %, gre razlika poveˇcanj na raˇcun izgub v zaslonu. V pri-42 4. Analiza zaslonskih struktur 80 40 -40 -80 V i 3-s ^ \x J ^^ -30 -15 0 15 x / cm 30 (b) Magnetni zaslon. 30 -15 -30 -15 0 15 30 x / cm (c) Prevodni magnetni zaslon. Slika 4.5: Gostotnice magnetnega polja dvovoda, ki je zastirano z zaslonom profila U. Sive linije prikazujejo gostotnice primarnega polja. pri prevodnem pa se nekoliko zmanjˇsa, kar je povsem v skladu s priˇcakovanji. 0 44 Poglavje 5 Sklep Primerjave rezultatov predlagane metode z rezultati drugih metod in primeri analize zaslonskih struktur izpričujejo, da omogoča razvita metoda soliden izračun polja pred, v in za zasloni in ponuja vpogled v učinke zastiranja. Metodologija je za vse objekte zaslonske strukture enotna in omogoča analizo magnetnih, prevodnih, prevodno-magnetnih in sestavljenih zaslonov ter je primerno orodje za uspešno dizajniranje zaslonov ob snopu vzporednih tokovodnikov. Pomanjkljivost metode je vsekakor njena omejenost na linearne snovi, kar onemogoča upoštevanje nelinearnih lastnosti feromagnetikov. Frekvečno mejo ji teoretično določa pogoj kvazistatičnosti, praktično pa vdorna globina, kompleksnost geometrije, potrebno število robnih elementov in seveda zmogljivost računalnika. Metodo bi kazalo nadgraditi z različnimi nabori poskusnih in testnih funkcij ter analizirati natančnost in hitrost izračuna v odvisnosti od nabora funkcij in velikosti robnih elementov. Pomanjkljivost metode je tudi njena omejenost z 2D geometrijo. Pri posplošitvi na 3D geometrijo bi se kompleksnost problema zelo povečala, saj bi oblogi enojnega in dvojnega sloja ne zadoščali. Vpeljati bi morali ekvivalentne električne in magnetne ploskovne toke na površinah zaslonov [53, razdelek 3.5], ki bi zahtevali neprimerno večje računalniške zmogljivosti. 45 Prispevki k znanosti Pomembnejši izvirni prispevki k znanosti so: 1. Uporaba potencialov enojnega in dvojnega sloja pri analizi prevodno-permeabilnih zaslonskih struktur v 2D prostoru. Kvazistatično elektromagnetno polje v in ob prevodno-permeabilnih objektih je izraženo s pomočjo enojnega in dvojnega sloja oziroma oblog, ki sta porazdeljeni po površinah objektov zaslonske strukture. S pomočjo Greenove formule in lastnosti potencialov enojnega in dvojnega sloja so izpeljane robne integralske enačbe za te vire. Prednost koncepta virov je predvsem v tem, da se vprašanje polja v strukturi prevede najprej na iskanje ekvivalentnih virov po ograjah objektov, kar je numerično vsekakor manjši zalogaj, nato pa se iz njih določita še polje in učinek zastiranja. 2. Določitev matričnega sistema robnih integralskih enačb, ki vključuje različna električna povezovanja elementov zaslonske strukture, in izračun oblog enojnega in dvojnega sloja po površinah prevodnih in permeabil-nih objektov. Za določitev magnetnega polja v prostoru niso pomembni le objekti kot takšni, ampak tudi načini ozemljitve in električne povezave med njimi. Te okoliščine povzemajo dodatni pogoji v obliki integralskih enačb, ki se naslavljajo na prej omenjene ekvivalentne vire. Celoten sistem robnih integralskih enačb je zatem strnjen v matrični sistem in obravnavan z momentno metodo. Pri tem je posebna pozornost posvečena numeričnim integracijam singularnih jeder robnih integralskih enačb. 3. Analiza povratnega učinka prevodnih, magnetnih in prevodno-magnetnih zaslonov na primarne vire polja. Prevodno-magnetni zasloni imajo neželen povratni vpliv na primarno strukturo, ki se kaže v spremenjenih karakteristikah. Prednost predlagane metode je ravno v tem, da obravnava vodnike in zaslone z enako metodologijo. Kot takšna že v osnovi omogoča 46 Prispevki k znanosti analizo interakcij med elementi vodniˇsko-zaslonske strukture, zato v pogledu analize povratnega uˇcinka ne predstavlja nadgradnje neke metode, ampak njo samo. 47 Izjava Izjavljam, da sem doktorsko disertacijo izdelal samostojno pod vodstvom mentorja doc. dr. Antona R. Sinigoja. Pomoˇc ostalih sem navedel v zahvali. V Ljubljani, 18. junij 2009 mag. Edi Bulic, univ. dipl. inž. el. 48 Zahvala Zahvaljujem se mentorju doc. dr. Antonu R. Sinigoju za strokovno vodstvo pri delu, ne-pogreˇsljive nasvete pri izdelavi doktorske disertacije in pripombe ter izboljˇsave v besedilu. Za mnoge odgovore pri postavljanju dela v sistemu LATEX se zahvaljujem bratu Patriciu ter sodelavcema Aleˇsu Berkopcu in Samotu Peniˇcu. Sodelavcem pri predmetih Osnove elektrotehnike, Elektromagnetika in Elektromagnetne strukture se zahvaljujem za prijetno delovno vzduˇsje ter vzpodbudo. In nenazadnje, za veliko podporo sem hvaleˇzen starˇsema Ivanu in Ivanki. 49 Literatura [1] A. P. Wills. On the magnetic shielding effect of tri-lamellar spherical and cylindrical shells. Physical Review, vol. 9, no. 4, pp. 193–213, Oct. 1899. [2] J. F. Hoburg. A computational methodology and results for quasistatic multilayered magnetic shielding. IEEE Transactions on Electromagnetic Compatibility, vol. 38, no. 1, pp. 92–103, Feb. 1996. [3] William M. Frix and George G. Karady. A circuital approach to estimate the magnetic field reduction of nonferrous metal shields. IEEE Transactions on Electromagnetic Compatibility, vol. 39, no. 1, pp. 24–32, Feb. 1997. [4] Marko Istenicˇ. Izpostavljenost in zaˇsˇcita ljudi in naprav v okolju elektriˇcnih in magnetnih polj elektroenergetskih omreˇzij. Magistrska naloga, Fakulteta za elektrotehniko Univerze v Ljubljani, 2000. [5] Massimo Guarnieri, Federico Moro, and Roberto Turri. An integral method for extremely low frequency magnetic shielding. IEEE Transactions on Magnetics, vol. 41, no. 5, pp. 1376–1379, May 2005. [6] Richard B. Schulz, Vellar C. Plantz, and David R. Brush. Shielding theory and practice. IEEE Transactions on Electromagnetic Compatibility, vol. 30, no. 3, pp. 187–201, Aug. 1988. Originally published in Proc. 9th Tri-Service Conf. on Electromagnetic Compatibility, Oct. 1963. [7] D. A. Miller and J. E. Bridges. Geometrical effects on shielding effectiveness at low frequencies. IEEE Transactions on Electromagnetic Compatibility, vol. 8, no. 4, pp. 174–186, Dec. 1966. [8] J. R. Moser. Low-frequency shielding of a circular loop electromagnetic field source. IEEE Transactions on Electromagnetic Compatibility, vol. 9, no. 1, pp. 6–18, Mar. 1967. 50 Literatura [9] Peter R. Bannister. New theoretical expressions for predicting shielding effectiveness for the plane shield case. IEEE Transactions on Electromagnetic Compatibility, vol. 10, no. 1, pp. 2–7, Mar. 1968. [10] Richard B. Schulz, Vellar C. Plantz, and David R. Brush. Low-frequency shielding resonance. IEEE Transactions on Electromagnetic Compatibility, vol. 10, no. 1, pp. 7–15, Mar. 1968. [11] John R. Harrington and Richard B. Schulz. Design of minimum weight and maximum effectiveness of very-low-frequency shielding. IEEE Transactions on Electromagnetic Compatibility, vol. 10, no. 1, pp. 152–157, Mar. 1968. [12] Robert B. Cowdell. New dimensions in shielding. IEEE Transactions on Electromagnetic Compatibility, vol. 10, no. 1, pp. 158–167, Mar. 1968. [13] Peter R. Bannister. Further notes for predicting shielding effectiveness for the plane shield case. IEEE Transactions on Electromagnetic Compatibility, vol. 11, no. 2, pp. 50–53, May 1969. [14] R. Yang and Raj Mittra. Coupling between two arbitrarily oriented dipoles through multilayered shields. IEEE Transactions on Electromagnetic Compatibility, vol. 27, no. 3, pp. 131–136, Aug. 1985. [15] J. F. Hoburg. Principles of quasistatic magnetic shielding with cylindrical and spherical shields. IEEE Transactions on Electromagnetic Compatibility, vol. 37, no. 4, pp. 574–579, Nov. 1995. [16] Yaping Du, T. C. Cheng, and A. S. Farag. Principles of power-frequency magnetic field shielding with flat sheets in a source of long conductors. IEEE Transactions on Electromagnetic Compatibility, vol. 38, no. 3, pp. 450–459, Aug. 1996. [17] Leonardo Sandrolini, Antonio Massarini, and Ugo Reggiani. Transform method for calculating low-frequency shielding effectiveness of planar linear multilayered shields. IEEE Transactions on Magnetics, vol. 36, no. 6, pp. 3910–3919, Nov. 2000. [18] Aldo Canova, Alessandra Manzin, and Michele Tartaglia. Evaluation of different analytical and semi-analytical methods for the design of ELF magnetic field shields. IEEE Transactions on Industry Applications, vol. 38, no. 3, pp. 788–796, May/June 2002. 51 Literatura___________________________________________________________________________ [19] Robert G. Olsen, Marko Istenič, and Peter Zunko. On simple methods for calculating ELF shielding of infinite planar shields. IEEE Transactions on Electromagnetic Compatibility, vol. 45, no. 3, pp. 538-547, Aug. 2003. [20] Marko Istenič. Zaˇsˇcita pred magnetnim poljem ekstremno nizkih frekvenc z magnetnimi zasloni konˇcnih dimenzij. Doktorska disertacija, Fakulteta za elektrotehniko Univerze v Ljubljani, 2003. [21] Marko Istenič and Robert G. Olsen. A simple hybrid method for ELF shielding by imperfect finite planar shields. IEEE Transactions on Electromagnetic Compatibility, vol. 46, no. 2, pp. 199-207, May 2004. [22] P. P. Silvester and R. L. Ferrari. Finite Elements for Electrical Engineers. Cambridge University Press, Cambridge, 1983. Reprinted 1986. [23] Matthew N. O. Sadiku. Numerical Techniques in Electromagnetics. CRC Press, Boca Raton, second edition, 2001. [24] L. Hasselgren, E. Möller, and Y. Hamnerius. Calculation of magnetic shielding of a substation at power frequency using FEM. IEEE Transactions on Power Delivery, vol. 9, no. 3, pp. 1398-1405, July 1994. [25] C. Buccella, M. Feliziani, F. Maradei, and G. Manzi. Magnetic field computation in a physically large domain with thin metallic shields. IEEE Transactions on Magnetics, vol. 41, no. 5, pp. 1708-1711, May 2005. [26] Aldo Canova, Giambattista Gruosso, and Maurizio Repetto. Integral methods for analysis and design of low-frequency conductive shields. IEEE Transactions on Magnetics, vol. 39, no. 4, pp. 2009-2017, July 2003. [27] Miodrag Milutinov and Neda Pekarič-Nad. Shielding effect of non -ferrous metallic plates in vicinity of three phase conductors. Serbian Journal of Electrical Engineering, vol. 2, no. 2, pp. 147-156, Nov. 2005. [28] Chih-Ping Chang and Chang-Fa Yang. A moment method solution for the shielding properties of three-dimensional objects above a lossy half space. IEEE Transactions on Electromagnetic Compatibility, vol. 47, no. 4, pp. 723-730, Nov. 2005. 52 Literatura [29] Peter Silvester. AC resistance and reactance of isolated rectangular conductors. IEEE Transactions on Power Apparatus and Systems, vol. 86, no. 6, pp. 770-774, June 1967. [30] Pei-bai Zhou. Numerical Analysis of Electromagnetic Fields. Springer-Verlag, Berlin, 1993. [31] Anton R. Sinigoj. ELMG polje. Fakulteta za elektrotehniko, Ljubljana, 1996. [32] Breda Cestnik, Primož Hrobat, Radomir Isakovič, Igor Babic, Zoran Strubel, Anton R. Sinigoj, Edi Bulic. Možnosti za znižanje jakosti elektromagnetnega sevanja v okolju naprav in objektov distribucijskega EE-omrežja. Referat ˇst. 1605, Elektroinˇstitut Milan Vidmar, Ljubljana, 2004. [33] Breda Cestnik, Rudi Vončina, Anton R. Sinigoj, Edi Bulic. Zaslanjanje magnetnega polja s prevodnimi in feromagnetnimi zasloni. V: Marjan Zorman (urednik), 26. posvetovanje o močnostni elektrotehniki in sodobnih električnih instalacijah Kotnikovi dnevi : izobraževalni program, pages 1/VII-14/VII, Radenci, Slovenija, 31. marec in 01. april 2005. Elektrotehniˇsko druˇstvo. [34] Edi Bulic. Izraˇcun magnetnega polja okrog dolgih prevodnih nemagnetnih zaslonov. Elektrotehniški vestnik, letn. 76, ˇst. 1-2, str. 31-37, 2009. [35] Ibrahim Tekin and Edward H. Newman. Moment method analysis of the magnetic shielding factor of a conducting TM shield at ELF. IEEE Transactions on Electromagnetic Compatibility, vol. 38, no. 4, pp. 585-590, Nov. 1996. [36] Willy P. Legros and Patrick G. Scarpa. Fast computation of magnetic field in rotationally symmetric structures. IEEE Transactions on Magnetics, vol. 21, no. 6, pp. 2644-2651, Nov. 1985. [37] Piergiorgio Alotto, Massimo Guarnieri, and Federico Moro. A boundary integral formulation on unstructured dual grids for eddy-current analysis in thin shields. IEEE Transactions on Magnetics, vol. 43, no. 4, pp. 1173-1176, Apr. 2007. [38] H. Igarashi, A. Kost, and T. Honma. A three dimensional analysis of magnetic fields around a thin magnetic conductive layer using vector potential. IEEE Transactions on Magnetics, vol. 34, no. 5, pp. 2539-2542, Sep. 1998. 53 Literatura [39] O. Bottauscio, M. Chiampi, D. Chiarabaglio, F. Fiorillo, L. Rocchino, and M. Zucca. Role of magnetic materials in power frequency shielding: numerical analysis and experiments. IEE Proceedings Generation, Transmission & Distribution, vol. 148, no. 2, pp. 104–110, Mar. 2001. [40] Hermann A. Haus and James R. Melcher. Electromagnetic Fields and Energy. Prentice-Hall, Englewood Cliffs, 1989. [41] Prem Kishore Kythe. An Introduction to Boundary Element Methods. CRC Press, Boca Raton, 1995. [42] C.Pozrikidis. A Practical Guide to Boundary Element Methods with the Software Library BEMLIB. Chapman & Hall/CRC Press, London/Boca Raton, 2002. [43] Roger F. Harrington. Field Computation by Moment Methods. IEEE Press Series on Electromagnetic Waves. IEEE Press, Piscataway, reissued edition, 1993. [44] C. A. Brebbia. The Boundary Element Method for Engineers. Pentech Press, London, 1978. [45] Julius Adams Stratton. Electromagnetic Theory. IEEE Press Series on Electromagnetic Wave Theory. IEEE Press, Piscataway, reissued edition, 2007. [46] John R. Carson. Wave propagation in overhead wires with ground return. Bell System Technical Journal, vol. 5, pp. 539–554, Oct. 1926. [47] Nagayoshi Morita, Nobuaki Kumagai, and Joseph R. Mautz. Integral Equation Methods for Electromagnetics. The Artech House Antenna Library. Artech House, Boston, 1990. [48] George B. Arfken and Hans J. Weber. Mathematical Methods for Physicists. Academic Press, San Diego, fifth edition, 2001. [49] Kenneth L. Kaiser. Electromagnetic Shielding. CRC Press, Boca Raton, 2006. [50] I. N. Bronˇstejn, K. A. Semendjajev, G. Musiol, H. Mu¨hlig. Matematiˇcni priroˇcnik. Tehniˇska zaloˇzba Slovenije, Ljubljana, 2. predelana in dopolnjena izdaja, 1997. [51] J. Van Bladel. Singular Electromagnetic Fields and Sources. IEEE Press Series on Electromagnetic Wave Theory. IEEE Press, Piscataway, 1996. 54 Literatura [52] The MathWorks, Inc. MATLAB® 7, Function Reference: Volume 2 (F-O), Sep. 2007. http://www.mathworks.com/access/helpdesk/help/pdf_doc/matlab/ refbook2.pdf. [53] Roger F. Harrington. Time-Harmonic Electromagnetic Fields. IEEE Press Series on Electromagnetic Wave Theory. IEEE Press, Piscataway, reissued edition, 2001. [54] Walter Gander and Walter Gautschi. Adaptive Quadrature – Revisited. Technical Report 306, Departement Informatik, ETH Zu¨rich, Aug. 1998. [55] The MathWorks, Inc. MATLAB® 7, Function Reference: Volume 3 (P-Z), Sep. 2007. http://www.mathworks.com/access/helpdesk/help/pdf_doc/matlab/ refbook3.pdf. [56] Anton R. Sinigoj. Osnove elektromagnetike. Fakulteta za elektrotehniko, Ljubljana, 1999. [57] Heinrich Kaden. Wirbelstr¨ome und Schirmung in der Nachrichtentechnik. SpringerVerlag, Berlin, 1959. [58] Lennart Hasselgren and Jorma Luomi. Geometrical aspects of magnetic shielding at extremely low frequencies. IEEE Transactions on Electromagnetic Compatibility, vol. 37, no. 3, pp. 409–420, Aug. 1995. [59] James R. Wait and David A. Hill. Electromagnetic shielding of sources within a metal-cased bore hole, vol. 15, no. 2, pp. 108–112, Apr. 1977. [60] James F. Hoburg, Bernard A. Clairmont, David W. Fugate, and Richard J. Lordan. Comparisons of measured and calculated power frequency magnetic shielding by multilayered cylinders. IEEE Transactions on Power Delivery, vol. 12, no. 4, pp. 1704–1710, Oct. 1997. [61] Salvatore Celozzi, Rodolfo Araneo, and Giampiero Lovat. Electromagnetic Shielding. Wiley Series in Microwave and Optical Engineering. John Wiley & Sons, Hoboken, 2008. [62] George W. Hanson and Alexander B. Yakovlev. Operator Theory for Electromagnetics. Springer-Verlag, New York, 2002. 55 Literatura [63] Constantine A. Balanis. Advanced Engineering Electromagnetics. John Wiley & Sons, Hoboken, 1989. [64] M. V. K. Chari and S. J. Salon. Numerical Methods in Electromagnetism. Academic Press Series in Electromagnetism. Academic Press, San Diego, 2000. [65] Jean G. Van Bladel. Electromagnetic Fields. IEEE Press Series on Electromagnetic Wave Theory. John Wiley & Sons, Piscataway, second edition, 2007. [66] Endre Su¨li and David F. Mayers. An Introduction to Numerical Analysis. Cambridge University Press, Cambridge, 2003. [67] Charles W. Steele. Numerical Computation of Electric and Magnetic Fields. Chapman & Hall, New York, second edition, 1997. [68] Andrew F. Peterson, Scott L. Ray, and Raj Mittra. Computational Methods for Electromagnetics. IEEE/OUP Series on Electromagnetic Wave Theory. IEEE Press, Piscataway, 1998. [69] Gabrijel Tomˇsicˇ. Izbrana poglavja iz matematike – Integralske enaˇcbe. Fakulteta za elektrotehniko, Ljubljana, 1998. [70] Izrail S. Gradshteyn and Iosif M. Ryzhik. Table of Integrals, Series, and Products. Academic Press, San Diego, fifth edition, 1994. Alan Jeffrey, editor. Translated from the Russian by Scripta Technica, Inc. [71] Zoya Popovic´ and Branko D. Popovic´. Introductory Electromagnetics. Prentice-Hall, New Jersey, 2000. [72] Clayton R. Paul. Introduction to Electromagnetic Compatibility. Wiley Series in Microwave and Optical Engineering. John Wiley & Sons, New York, 1992. [73] M. H. Lean and A. Wexler. Accurate field computation with the boundary element method. IEEE Transactions on Magnetics, vol. 18, no. 2, pp. 331–335, Mar. 1982. [74] Ruey-Beei Wu and Jin-Chum Yang. Boundary integral equation formulation of skin effect problems in multiconductor transmission lines. IEEE Transactions on Magnetics, vol. 25, no. 4, pp. 3013–3015, July 1989. 56 Literatura [75] Antonije R. DjordjeviČ, Tapan K. Sarkar, and Sadasiva M. Rao. Analysis of finite conductivity cylindrical conductors excited by axially-independent TM electromagnetic field. IEEE Transactions on Microwave Theory and Techniques, vol. 33, no. 10, pp. 960-966, Oct. 1985. [76] The MathWorks, Inc. MATLAB® 7, Function Reference: Volume 1 (A-E), Sep. 2007. http://www.mathworks.com/access/helpdesk/help/pdf_doc/matlab/ refbook.pdf. 57 58 IEEE TRANSACTIONS ON ELECTROMAGNETIC COMPATIBILITY, VOL. 51, NO. 2, MAY 2009 263 An Equivalent Surface Source Method for Computation of the Magnetic Field Reduction of Metal Shields Edi Bulic´, Anton R. Sinigoj, and Breda Cestnik Abstract—A numerical method for computation of the resultant quasi-static magnetic field in the vicinity of parallel wires and metal shields is presented. The primary magnetic field source is time-harmonic currents in wires. This field is modified by conducting magnetic and/or nonmagnetic shields. The material is assumed to be linear under the applied source field. The shielding effectiveness can be estimated by a comparison between the primary and the resultant field. The reaction magnetic field is expressed by a sum of fields caused by equivalent single- and double-layer sources distributed on the shield surface. Integral equations for unknown distributions of these equivalent sources are derived from the Green’s second identity implemented inside and outside the shields. These equations are coupled integral equations, and are solved by the moment method. Numerical results of the resultant (shielded) magnetic field obtained with the proposed method are compared with the results of: 1) analytically solvable problems; 2) measurements; and 3) two different numerical methods. Index Terms—Boundary element methods, equivalent sources, integral equations, magnetic shielding, metal shields, moment methods, single and double layers. I. Introduction IN ORDER to meet nonionizing radiation requirements and/or electromagnetic compatibility criteria, magnetic field shielding is frequently needed. The disturbing magnetic field in a particular region can be minimized by shifting field sources and/or by setting metal shields between the sources and the region. Shielding effectiveness is impacted by many of the shield parameters, such as position, size, thickness, shape, and electromagnetic properties. To design shields efficiently, a method is needed for simulating shields of various shapes and electromagnetic properties, and for calculating the magnetic field in the shield vicinity. There are many methods proposed in literature. Analytical estimations of shielding efficiency [1]–[5] are based on various assumptions, such as extremely low frequency (ELF; 3 Hz–3 kHz [6]), flat shields, infinite planar shields, and similar, and are not suitable for general use. When using the finite-element method [7], which deals with various geometrical and material properties, there are relatively many elements required for calculating the resultant field in the shield vicinity [4]. By using Manuscript received January 31, 2008; revised September 4, 2008. First published March 31, 2009; current version published May 15, 2009. e. Bulic´ and a. r. Sinigoj are with the Faculty of Electrical Engineering, University of Ljubljana, Ljubljana 1000, Slovenia (e-mail: edi.bulic@fe.uni-lj.si; anton.sinigoj@fe.uni-lj.si). b. Cestnik is with Milan Vidmar Electric Power Research Institute, Ljubljana 1000, Slovenia (e-mail: breda.cestnik@eimv.si). Digital Object Identifier 10.1109/temc.2009.2014848 integral methods, difficulties caused by the numerous elements, such as the needed computational time and computer memory consumption, are reduced. One such method is the multicon-ductor method (MCM) [6], [8], [9]. It is suitable for conducting nonmagnetic shields. Another one is the surface-current method (SCM) [10] used for magnetic shields. It is drawn from a mag-netostatic approximation, and is therefore valid only for ul-tralow frequencies (ULF; below 3 Hz [6]) and conditionally for ELF. In this paper, we propose an equivalent surface source method (ESSM) enabling calculation of the magnetic field in the vicinity of parallel wires, and conducting magnetic and/or conducting nonmagnetic shields of practically any shape, thus allowing for a wide spectrum of applicability. The ESSM is based on a quasi-static description of the field. The primary field source is time-harmonic, and as all the materials are linear, the field is time-harmonic too; therefore, the field analysis is made in the frequency space, and the field quantities are considered as phasors. II. Equivalent Surface-Source Concept A. Structure Geometry and Problem Presentation The specific geometry of the problem to be examined is shown in Fig. 1. The primary magnetic field source is time-harmonic currents la, Ic2, ■ ■ ■, Icnc in long thin parallel conductors (wires), where nc is the number of the conductors. Homogeneous conducting magnetic and/or nonmagnetic shields run parallel to the wires, and modify the primary field. The shield cross-sectional areas are A = ezAz(x, y). Under quasi-static conditions, the potential at some point T(x, y) is governed by Pois-son’s equation [11] AAZ(T) = \ TT is the dis KT)JZ(T), ng TeS TeS0 (1) where TTa is the distance from point T to point Ta and S2 is the Dirac delta function in 2-D. The current density in the shields can be expressed by the electric field intensity E or by the magnetic vector potential A and electric scalar potential V Jz = lEz = -]^AZ 7^ = -J^7^+7^ (2) where j = a/—T, u> is the angular frequency of the time-harmonic field, and U is the electric potential difference between the beginning and end of the shield (the direction of the z-axis defines which end of the shield is the beginning one and which is the ending one). The potential difference U is constant over the cross-sectional area of one shield, but, in general, it varies from shield to shield. Applying (2) into (1) yields U(T) AAz(T)+kz(T)Az(T) = -M(T)7(T) /o AAZ(T)= -J2^oIciÖ2 (TTa), T e S TeS (3a) (3b) 8=1 where k = a/— jw/x7 is the wavenumber of the (metal) shield. Equation (3a) is the Helmholtz’s equation [12]. Equation (3a) is solved with the boundary element method [13]-[15]. The boundary integral equations are derived as follows. The Green’s functions inside and outside the shields for the 2-D Helmholtz’s and Poisson’s equations, respectively, are [16] 1 (2) g(T,r)= -H0 3o(T,T')= -ln= (fc(T)TT7) 1 (4a) (4b) where Hq is the zero-order Hankel function of the second kind. Inserting the Green’s functions (4) and magnetic potential z-component Az into the Green’s second theorem, we obtain f (g(T, T')AAZ(T) - Az(T)Ag(T,T')ds = I (V,r')^(T-) ^(r_}9g(r,ro\ / (do (T, T')AAZ (T) - Az (T)Ag0 (T, T')) *(T+) ds dgo(T,V) + Az(T+y" d/ dn dr (5b) where T' is any point and T is the integration point. Taking (3) into consideration and the boundary conditions for tangential components of the magnetic and electric field intensities on the shield boundaries 9AJT+) 1 dAz(T) /r () -------d^ = -^(T-) dn =aiiT) m AZ{T+)= AZ (T_)= a2 (T) (6b) where yur is the relative permeability, we have f (az(T)ö2 (TT) -M(T)7(T)^^ff(T,T'))ds = I (-^(T)ai(T)g(T,T')-a2(T)^g(J'V) \d/ (7a) ng -J2go(Tci,T')i^Ici+ AZ(T)S2(TT .--1 S ds fL(T)go(T,r) + Q2(T)%» (T,r) dn d/. (7b) Authorized licensed use limited to: UNIVERSITY OF LJUBLJANA. Downloaded on June 15, 2009 at 05:52 from IEEE Xplore. Restrictions apply. BULIC et al.: EQUIVALENT SURFACE SOURCE METHOD 265 To find the limit of (7a) and (7b) as the point T" approaches the shield boundaries C from outward (T" —> T'+) and from inward (T" —> Ti), respectively, we take into account the results 1 from [17, Sec. 3.1]. The resulting equations are coupled integral equations for functions a.\ and a2, which are density distributions of the single layer and double layer, respectively, on the shield boundaries a2{T') + 2 I ( md9(.T>T') + fr(T)ai(T)g(T,T')dl )ds =0, T' e C (8a) s 'o a2 (Tr) -2f (a2 (T)9go(T,r) + ai (T)go (T, T')) d; rac =2^g0(Tci,T')Mo-fci, T'e £. (8b) These two layers are assumed to be equivalent surface sources of the electromagnetic field. Surface integrals over the shield cross-sectional areas can be transformed into curve integrals over their boundaries by the Gauss’s theorem / MT)7(T)^W ;)ds ^Sj '0 cSg, MSj 7Sj ^ 1 / dg{T,T'+) S j Csj dn dl (9) where U$j is the electric potential difference between the beginning and end of the jth shield and k$j is the wavenumber of the jth shield. Explicit expressions for the vector magnetic potential at a point inside or outside the shields are derived from (7a) or (7b), respectively + fr(T)ai(T)g(T,T')\dl + f KTh(Tm9(T,r)ds s k (10a) + )+ a2(T)-------^ dl — L V + J Lr(T)ai(T)5(T,T')+ a2(:r) / (a1(T)gQ(T,r)+ a2(T)d9°^'T-) dl L V dn dg(T,T') dn (ai(T)5o(T,r') + a2(T) L ^ dg0(T,T') dn dl, dl, T' eC a2(T') T' eC AZ{T' eS0 £a2{T)d90^r) + ai(T)go(T,T')\ + J29o(,Tci,T')^IGi. dl (10b) The Green’s function outside the shields (4b) must be modified when the wires and shields run parallel to the ground [18] 1 r\ P-W4 g0(T,T') = — In —- (11) where dc — 1.85/\kg\ is the Carson’s distance and kg is the wavenumber of the ground. If the Carson’s distance is much larger than the distances (TT') from the wires and shields to the observed points, the impact of the ground effect can be neglected. This condition is usually fulfiled at low frequencies when the Carson’s distance is quite large. Shields can have open ends or can be grounded for safety reasons. If the jth shield is grounded at its beginning and its end (two-sidedly grounded), an additional current path is provided for current circulation and the net-induced current I$j in that shield is not equal to zero [5]. The current I$j is proportional to the integral of the single-layer density over the boundary Cgj, according to Ampere’s law, and also to Usj, according to Ohm’s law (for an ac circuit) *y i Mo £sj ai(T)d/ (12) where Zgisj is the sum of the grounding impedances of the jth shield. This equation forms a system of equations with (8) because Usj is also an unknown variable. If Zg$j is very small compared to the shield impedance, grounding can be considered ideal. In this case, Us j —> 0 and (12) does not need to join (8) because Us j is no longer unknown. If the jth shield is grounded only at one point (one-sidedly grounded) or has open ends, no additional current path is provided and the net-induced current Is j is equal to zero hi = — ai(T)d/ =0. Mo JcSj (13) Shields can also be electrically connected. If, for example, the jth and the ith shield are connected at least at their beginnings and ends, their electrical potential differences are the same, i.e., Us j = Usi. When they are also one-sidedly grounded, their common net-induced current is equal to zero, and the equation hj+Isi = — Mo Mo cSjucSi ai(T)d/ =0 (14) joins (8). C. Numerical Method In our paper, the system of coupled integral equations [(8) and equations like (12), (13), and/or (14)] is solved by using 2 Authorized licensed use limited to: UNIVERSITY OF LJUBLJANA. Downloaded on June 15, 2009 at 05:52 from IEEE Xplore. Restrictions apply. 266 IEEE TRANSACTIONS ON ELECTROMAGNETIC COMPATIBILITY, VOL. 51, NO. 2, MAY 2009 TABLE I CPU Time ™BE TEss (min) Tb (s) 27 0.14 1.0 54 0.53 1.9 99 1.74 3.5 201 7.09 6.9 401 29.2 14.2 801 95.9 27.6 Fig. 2. Coaxial shield. the moment method [19]. The boundaries of the shield cross-sectional areas are discretizated into short straight segments (boundary elements) of lengths AZi, AZ2, • • •, A/nBE, where nBE is the number of the segments. The pulse functions and Dirac’s delta functions are used as expansion and testing functions, respectively. The integrals of the Green’s functions and their normal derivatives over the boundary elements are computed by using Simpson’s rule. Solution to the system of coupled integral equations are single-layer and double-layer densities. From this solution, we first compute the vector magnetic potential in the nodes of a rectangular grid [by using (10)], and then the magnetic flux density from B = V x A by using finite differences. To give an idea of the computational effort for various numbers of the boundary elements, the CPU times Tess and T# required to calculate: 1) density distributions of the single layer and double layer, and 2) magnetic flux density in one point, respectively, are shown in Table I. Computations were performed by a nonoptimized MATLAB code on a personal computer equipped with an Intel 1.83-GHz processor and 2 GB of RAM. III. Verification of Results The results obtained with the ESSM are compared with: 1) results of analytically solvable problems; 2) experimental results; 3) results obtained with two different numerical methods. A. Cylindrical Shields 1) Coaxial Shield: The case in Fig. 2 corresponds to the geometry of the coaxial cable. The time-harmonic current Iq in the inner conductor is the source of the primary magnetic field. The outer conductor plays the role of a shield and has open ends. The function of radial distribution of the induced current density J = ezJz (r) in the shield, which is source-free (J = 7s E) with the wavenumber k$, is governed by the Bessel’s differential equation [12] AJz(r) + klJz(r) d2Jz H-----------V k$Jz{r)=0 (15) and is expressed with the zero-order Bessel functions of the first and second kind Jz(r) = d Jo(ksr) + C2N0(ksr). (16) The constants C\ and C2 result from the boundary conditions on the surfaces of the outer conductor, iJ^(ri) = Io/(2irri) and Htp (r2) = Io/(27rr2), and have the following values: Ci c2 ksIo N1(ksr2)/r1 ^ Ji(fcsri)M(fcsr2) fcs/o Ji(ksr1)/r2 27T J1(ksr1)N1(ksr2) JVi(fcsri)/r2 Ji(fcsr2)M(fcsri) Ji(ksr2)/r1 Ji(fcsr2)iVi(ifesri) (17a) (17b) where J\ and N\ are the first-order Bessel functions of the first and second kind, respectively. The net-induced current in the shield is equal to zero due to its open ends (but this does not imply that the induced current density is equal to zero). Beacuse of this and since the current density is axial-symmetrically distributed in the shield, such a shield modifies the primary field only in its own interior, where T\ 50 mm or r —> 55 mm, i.e., at the points of the single layer and double layer. The analytical solution in the case of the coaxial shield can be used in estimation of accuracy of the proposed numerical method. The discrepancy between the numerically and analytically computed value of the magnetic field is dependent primarily on the ratio between the length Al of the boundary elements and the skin depth S = y/2/(cj/xs7s) of the shield (to allow for straightforward dependence, all the elements are of the same length). This dependence is analyzed. The numerical data were f = 250 Hz, v\ = 10 mm, r2 = 20 mm, 7s = 10 MS/m, and /xs = 100 /xo. Fig. 5 illustrates absolute values of the relative discrepancy in the shield at r = 15 mm. The discrepancy increases progressively with the ratio Al/5, and at Al/5 = 1, it does not yet reach 1%. 2) Shielded Pair of Conductors: Consider a bifilar pair of parallel conductors, which is centered inside a cylindrical shield Authorized licensed use limited to: UNIVERSITY OF LJUBLJANA. Downloaded on June 15, 2009 at 05:52 from IEEE Xplore. Restrictions apply. BULIC et al.: EQUIVALENT SURFACE SOURCE METHOD 267 Fig. 3. Comparison with the analytical solution for the coaxial nonmagnetic-metal shield. (a) Magnetic field intensity. (b) Discrepancy. Fig. 4. Comparison with the analytical solution for the coaxial magnetic-metal shield. (a) Magnetic field intensity. (b) Discrepancy. Fig. 5. Dependence of accuracy of the magnetic field intensity calculation on the ratio between the length of the boundary elements and the skin depth for the coaxial shield. (see Fig. 6). The field in this case is analytically solvable. The analytical solution for a shield with a small thickness compared to the shield skin depth {r r2 (19c) where f£ 2n x and i?2 ri- are the Hankel functions of the first and second kind, respectively, of order (2n — 1) and fcg is the wavenumber of the shield. The constants Ci?n, C2,n, C^, and C4?n result from the boundary conditions on the shield surfaces, which state continuity of tangential components Az and H^ of the magnetic vector potential and magnetic field intensity, respectively. As series in (19) converge very rapidly, only the first few summands need to be considered. We considered the first ten summands. A comparison of the radial dependence of the magnetic field intensity magnitudes, along the straight line at cp = 7r/4, obtained with the ESSM and with the analytical approach, is shown in Fig. 7. The only significant deviation from the analytical solution is again (see Figs. 3 and 4) detected at the shield surfaces. The analytical solution is once more used in estimation of accuracy of the proposed numerical method. The analysis of dependence of the discrepancy between the numerically and analytically computed values of the magnetic field on the ratio Al/5 is made for the same shield and same frequency as in the case of the coaxial shield. The primary field source is a pair of conductors whose spacing is 2d = 10 mm. Fig. 8 shows the results of this analysis in the point (r = 15 mm, cp = 7r/4). The discrepancy again increases progressively with the ratio Al/5, and at Al/5 ~ 1, it is near 1%. Hoburg et al. [24] measured the magnetic field from a pair of conductors surrounded by a cylindrical shield. Fig. 9 provides a comparison between the shielding factor values obtained: 1) analytically; 2) with the ESSM; and 3) with the experiment (in [24]) for a magnetically nonpermeable (Al) shield. Fig. 10 provides the same comparison for a magnetically permeable (low carbon steel) shield. The numerical data were: 1) for the primary field source Iq = 500 A, / = 60 Hz, and 2d = 6 cm; 2) for Fig. 8. Dependence of accuracy of the magnetic field intensity calculation on the ratio between the length of the boundary elements and the skin depth for a cylindrical shield surrounding a pair of conductors. Fig. 9. Comparison with the analytical solution and experimental results for a pair of conductors inside a cylindrical nonmagnetic-metal shield. the Al shield v\ = 25.435 cm, r^ =25.5 cm, 7s =30.5 MS/m, and /is = Mo; and 3) for the steel shield v\ = 25.43 cm, r^ = 25.5 cm, 7s =6.76 MS/m, and /is = 190 /io2. The experiment 2In [24] it reads: “The measured initial relative permeability of the steel was /ir = 182.5, and the relative permeability remained below /ir = 190 for flux densities in the steel below 10-3 T.” But the flux density in the steel shield is about 15 mT, so in calculations, we used the value of 190 for the relative permeability. The relative permeability is likely larger than 190, and this probably also causes the discrepancy between the measured and calculated results. The proposed method is applicable only to linear materials, but there are methods that can deal with the nonlinearity of ferromagnetic materials [25]. Authorized licensed use limited to: UNIVERSITY OF LJUBLJANA. Downloaded on June 15, 2009 at 05:52 from IEEE Xplore. Restrictions apply. BULIC et al.: EQUIVALENT SURFACE SOURCE METHOD 269 Fig. 10. Comparison with the analytical solution and experimental results for Fig. 12. Comparison with the analytical approximation for an exceedingly a pair of conductors inside a cylindrical magnetic-metal shield. wide flat nonmagnetic-metal shield. Fig. 11. Exceedingly wide flat shield above the wire. setup is described in detail in [24]. The shielding factor is calculated at a radial distance of r = 0.62 m for angles 0 < cp < ir/2. It was measured at the same distance for two angles, cp = 0 and if = 7r/2. Figs. 9 and 10 show that the ESSM yields an excellent agreement (within 0.02 %) with the analytical solution. Experimental results for the Al shield deviate from the analytical one by —9.0 % and 0.0 % for angles cp = 0 and cp = 7r/2, respectively. These deviations for the steel shield are about 4.7 % and 7.8 %. B. Exceedingly Wide Flat Shield We deal with the case of an exceedingly wide and two-sidedly, ideally grounded shield above the wire shown in Fig. 11. The numerical data are I0 = 100 A, / = 50 Hz, t = 2 mm, ys = 5 cm, and / = 4 m. To estimate the reaction field of such shield, we use the analytical solution for the infinite planar shield [1], [18]. The analytical solution for the vector magnetic potential above the shield is expressed as the Fourier cosine transform OO __________ Az (x,y>ys+t)= / /(0e_ v^-^ž * o cos(£c)d£ (20) where fco is the wavenumber of free space (or air) and the function /(£) is obtained from the boundary conditions of the electromagnetic field on the lower and upper surface of the shield, and is given by TV V x (2 /is^oyS2 £2etV^I C kl cosh t ^ 0 + Ms (s^-fco) +f4 (C2-*!) smht^/e-kl Fig. 13. Comparison with the analytical approximation for an exceedingly wide flat magnetic-metal shield. (21) Fig. 14. Finite-width flat shield above the wires. The integral in (20) is solved numerically. Comparisons of the magnetic flux density magnitudes above the shield (at the height of y^ = 15 cm above the wire) obtained with the ESSM and with analytical approximation are shown in Figs. 12 and 13 for the nonmagnetic (7s =30.5 MS/m and Ms = Mo) and magnetic shield (7s =10 MS/m and /is = 10 /xq ), respectively. The wide shield resembles the infinite planar shield in the vicinity of its center. The discrepancy between the wide shield and infinte planar shield increases rapidly with the distance from the center. C. Finite-Width Flat Shields Our comparison between the results of several methods is made for a finite-width flat shield placed above the wires (Fig. 14). Authorized licensed use limited to: UNIVERSITY OF LJUBLJANA. Downloaded on June 15, 2009 at 05:52 from IEEE Xplore. Restrictions apply. 270 IEEE TRANSACTIONS ON ELECTROMAGNETIC COMPATIBILITY, VOL. 51, NO. 2, MAY 2009 1) Hybrid Method (HM): In [3], a method for calculating the shielding performance of a 2-D finite-width metal shield is described. It is based on an analytical approximation for the magnetic field behind the shield. This approximation combines the expressions for penetration of the field through the infinite planar shield and for leakage around the finite-width shield, and is therefore called HM. 2) Multiconductor Method (MCM): This method [6], [8], [9] is suitable only for analysis of the nonmagnetic shields. It is based on the solution of the integral equation for distribution of the induced current density J = ez Jz in shields. The equation is derived as follows. The sources of the magnetic field are currents in wires and shields MT) Mo 2tt £/ciln 1 TT a + ^ j Jz(T')h n JL= ds'. (22) If the shields and wires run parallel to the ground, the Green’s functions in this equation have to be modified by means of (11). By inserting (22) in (2), the following integral equation is obtained: JZ(T) jwno [ Jz{T')\ s n JL= ds' U(T) ^T i=l ( (23) a The electric potential differences U between the beginning and end of the shields are also unknowns in the integral equation. Thus, to obtain its unique solution, (12) or (13) for the two-sidedly or one-sidedly grounded shields, respectively, have to be added hj Jz{T')ds' = Sj 0 3 = 1,2, (24) In solving numerically the acquired system of coupled integral equations (23) and (24) by using the moment method [19], the cross-sectional areas of the shields are discretizated. It is advantageous to use the MCM rather than the ESSM for thin shields (requiring less unknowns, i.e., discretization elements, for the same accuracy of results) and, inversely, for thick shields. Frix and Karady [6] measured the magnetic field from a pair of wires (2c i = ^o and 2c3 = — 2q, the central wire in Fig. 14 is absent) near a flat Al shield. The measurement was performed for two different Al shields. In the first case, the numerical data were 20 = 100 A, / = 60 Hz, ys = 6.35 cm, 2d = 7.62 cm, I = 30.48 cm, t = 3.175 mm, 7s = (103/36.99) MS/m, and Ms = Mo,the shield had open ends, and the magnetic flux density was measured at the height of ^r = 14.13 cm above the wires. In the second case, the only changed parameters were thickness (t = 6.35 mm) and conductivity (7s = (103/45) MS/m) of the shield and the height of measurement points (y*& = 14.45 cm). The experiment setup is described in detail in [6]. We made a comparison between the results obtained with the experiment, ESSM, MCM, and HM. Figs. 15 and 16 show: 1) that the ESSM yields an excellent agreement with the MCM (the average absolute value of the relative difference is emean abs. = 0.03 %) and Fig. 15. Comparison with the results of the experiment, MCM, and HM for the finite-width flat Al shield of the thickness oft = 3.175 mm. Fig. 16. Comparison with the results of the experiment, MCM, and HM for the finite-width flat Al shield of the thickness oft = 6.35 mm. 2) that the agreement of these two methods with the experiment (čmean abs. = 5.1 %) is much better than that of the HM with the experiment (emean abs. = 16.2 %). 3) Surface-Current Method (SCM): This method is suitable only for analysis of magnetic shields. It is based on the magne-tostatic approximation. On surfaces of the magnetic shields, the magnetization surface-current density K is proportional to the tangential component of the magnetic flux density [10] K(T) + —ß(T)n(T) x B(T) = 0 Mo (25) where ß(T) = (fi(T) - /i0) /(/x(T) + /x0). By expressing the magnetic flux density B in this equation as a sum of the primary field (Bo) and contribution of the magnetization surface-current density K = ez Kz, an integral equation for the distribution of the latter on the shield surfaces is obtained Kz{T) + ^Uz{Ti)^,l> IT£ PZ 2ß(T) Mo ez.(n(T) xB0(T)) (26) where P is the distance vector between points T' and T. A comparison between the results of the three methods, namely ESSM, SCM, and HM, is made for a two-sidedly (ideally) grounded flat magnetic-metal shield placed above wires of a three-phase system (see Fig. 14): 2ci = Iq exp(j27r/3), 2c2 = 2o, and 2c3 = Iq exp(— j2tt/3). Results for the magnetic flux density at the height of ?/r = 40 cm above the Authorized licensed use limited to: UNIVERSITY OF LJUBLJANA. Downloaded on June 15, 2009 at 05:52 from IEEE Xplore. Restrictions apply. BULIC et al.: EQUIVALENT SURFACE SOURCE METHOD 271 Fig. 17. Comparison with the SCM and HM for the finite-width flat magnetic- ld Fig. 19. Comparison with the experimental results for the finite-width flat nonmagnetic-metal shield at y^ =30 cm. Fig. 18. Comparison with the SCM and HM for the finite-width flat magnetic- f Fig. 20. Comparison with the experimental results or the finite-width flat metal shield at an ULF. nonmagnetic-metal shield at ?/r =70 cm. wires are given in Figs. 17 and 18. For the case in Fig. 17, the numerical data are Iq = 400 A, / = 50 Hz, y$ = 10 cm, d = 8 cm, I = 40 cm, t = 2 mm, 7s = 10 MS/m, and /iS = 1000 /xq. The HM agreement with the ESSM is not quite good (čmean abs. = 9-5%). The SCM agreement with the ESSM is even worse (emean abs. = 13.4 %) because the magnetostatic approximation does not consider the induced eddy currents. This is why a second test (see Fig. 18) is made where only the frequency is changed (/ = 1 Hz). This time, the SCM agreement is excellent (emean abs. =0.09 %) and better than the HM one (čmean abs. = 1-1 %). This proves that the magnetostatic approximation is justifiable only when the impact of the induced eddy currents is negligible, this being at a very low frequency [6] (or for laminated magnetic shields and ferrite shields). 4) Experimental Results: Our experiment is made for the configuration shown in Fig. 14. The two-sidedly grounded flat shield is I = 1 m wide and 8 m long, and is placed ys = 20 cm above the wires of a three-phase system: 2c 1 — ^0 exP(j27r/3), ic2 = IoJcz = 20exp(-j27r/3),20 = 100A,/ = 50 Hz,and d = 8 cm. Grounding is considered to be ideal in the numerical calculation. Two types of shields are used: 1) t = 4 mm thick Al shield (7s =30.5 MS/m and /ig = /xq) and 2) t = 1.6 mm thick steel shield (7s =10 MS/m and /ig =100 /xq). The magnetic flux density is measured by the Wandel and Goltermann EFA-3 electromagnetic (EM) field analyzer. A comparison between the measured and calculated magnetic flux density magnitudes above the shield is made at the heights of Vr — 30 cm and ?/r = 70 cm above the wires. Figs. 19 and Fig. 21. Comparison with the experimental results for the finite-width flat magnetic-metal shield at 2/r = 30 cm. 20 provide a comparison for the nonmagnetic (Al) shield, and Figs. 21 and 22 for the magnetic (Fe) shield. The average absolute values of the relative deviation of the measured results from the calculated one are 5.5 %, 9.6 %, 12.0 %, and 4.7 % for the data on Figs. 19, 20, 21, and 22, respectively. IV. Presentation of the Results The ESSM enables a prominent presentation of its results. How the calculated field or how the shielding factor can be presented is demonstrated on the configuration shown in Fig. 23. The shield is two-sidedly (ideally) grounded. The numerical data are /0 = 400 A, / = 50 Hz, ys = 13 cm, d = 8 cm, / = 70 cm, h = 40 cm, t = 1 cm, 7s = 10 MS/m, and /is = 100 /xq (steel). Authorized licensed use limited to: UNIVERSITY OF LJUBLJANA. Downloaded on June 15, 2009 at 05:52 from IEEE Xplore. Restrictions apply. 272 IEEE TRANSACTIONS ON ELECTROMAGNETIC COMPATIBILITY, VOL. 51, NO. 2, MAY 2009 Fig. 22. Comparison with the experimental results for the finite-width flat magnetic-metal shield at yR = 70 cm. Fig. 23. Shield of a U-shape. The wires are a three-phase system: 7c i 7/0exp(j27r/3),7/C2 =/0,and/C3 =/0 exp(-j27r/3). Fig. 24. Magnetic field around a U-shaped shield in the case of one-phase excitation. Fig. 24 illustrates the magnetic flux lines in the presence of only the middle wire (one-phase excitation). It is evident that the ferromagnetic shield “catches” a lot of the magnetic flux (flux shunting). Fig. 25 shows contours of the constant shielding factor values for the case of three-phase excitation. The reason why we here draw the contours of the shielding factor values from 1 to 20 alone is because they are sufficient to show the shielding tendency and because the contours of the higher shielding factor values are getting thicker. Attenuation of the magnetic field is evidently greater than ten times almost everywhere above the shield (y > ys -\-1). Let us now add another shield. It is flat and 80 cm wide (thickness and material parameters are the same as for the first Fig. 25. Shielding factor around a U-shaped shield in the case of three-phase excitation. Fig. 26. Shielding factor around a two-shield structure in the case of three-phase excitation. shield). When placing this additional shield horizontally below the first one, so that the wires are surrounded by the shields, the shielding factor is considerably increased everywhere and not only below the shields, as illustrated in Fig. 26 (showing only the contours of the shielding factor values from 5 to 100). The worst shielding effect is in the vicinity of the slots between the shields. Presentation of the calculated field, like in Figs. 24, 25, and 26, offers an excellent insight into the shielding effect. Such analysis is found to be a very useful tool in designing shielding structures. V. C ONCLUSION By comparing the ESSM with several other methods and by presenting its results, it is clear that it allows for a highly accurate calculation of the magnetic field either before, inside, or behind shields, thus providing good insight into the effectiveness of a shield we intend to choose. The ESSM is suitable for lossy nonferrous as well as ferrous shields provided the material media are linear. Its methodology is the same for the nonmagnetic and magnetic shields meaning that structures that contain either of the two can be analyzed. Authorized licensed use limited to: UNIVERSITY OF LJUBLJANA. Downloaded on June 15, 2009 at 05:52 from IEEE Xplore. Restrictions apply. BULIC et al.: EQUIVALENT SURFACE SOURCE METHOD 273 The proposed method does not impose any restrictions on the geometry of the longitudinally symmetric structures or the number of shields. Unlike many other methods, our method is not limited to the ELF fields. Though the frequency is theoretically limited by the quasi-static field assumption, it is in practice limited by computer capability, i.e., by the number of the boundary elements the computer can deal with. This number depends on the admissible boundary element length that is related to the frequency-dependent skin depth. Until now, none of the methods presented in literature have considered the impact of the ground vicinity, different types of grounding, and electric connections between shields. Judging from the facts presented before, the proposed method can be considered as one of the most general and efficient 2-D calculation method for estimation of the magnetic field shielding effectiveness. [19] R. F. Harrington, Field Computation by Moment Methods. Piscataway, NJ: IEEE Press, 1993. [20] H. Kaden, Wirbelstro¨me und Schirmung in der Nachrichtentechnik. Berlin, Germany: Springer-Verlag, 1959. [21] J. R. Wait and D. A. Hill, “Electromagnetic shielding of sources within a metal-cased bore hole,” IEEE Trans. Geosci. Electron., vol. 15, no. GE-2, pp. 108–112, Apr. 1977. [22] J. F. Hoburg, “Principles of quasistatic magnetic shielding with cylindrical and spherical shields,” IEEE Trans. Electromagn. Compat., vol. 37, no. 4, pp. 574–579, Nov. 1995. [23] J. F. Hoburg, “A computational methodology and results for quasistatic multilayered magnetic shielding,” IEEE Trans. Electromagn. Compat., vol. 38, no. 1, pp. 92–103, Feb. 1996. [24] J. F. Hoburg, B. A. Clairmont, D. W. Fugate, and R. J. Lordan, “Comparisons of measured and calculated power frequency magnetic shielding by multilayered cylinders,” IEEE Trans. Power Del., vol. 12, no. 4, pp. 1704– 1710, Oct. 1997. [25] O. Bottauscio, M. Chiampi, D. Chiarabaglio, F. Fiorillo, L. Rocchino, and M. Zucca, “Role of magnetic materials in power frequency shielding: Numerical analysis and experiments,” in Proc. Inst. Elect. Eng. Gener. Transmiss. Distrib., Mar. 2001, vol. 148, pp. 104–110. References [1] R. G. Olsen, M. Istenicˇ, and P. Zˇunko, “On simple methods for calculating ELF shielding of infinite planar shields,” IEEE Trans. Electromagn. Compat., vol. 45, no. 3, pp. 538–547, Aug. 2003. [2] R. G. Olsen and P. Moreno, “Some observations about shielding extremely low-frequency magnetic fields by finite width shields,” IEEE Trans. Elec-tromagn. Compat., vol. 38, no. 3, pp. 460–468, Aug. 1996. [3] M. Istenicˇ and R. G. Olsen, “A simple hybrid method for ELF shielding by imperfect finite planar shields,” IEEE Trans. Electromagn. Compat., vol. 46, no. 2, pp. 199–207, May 2004. [4] L. Sandrolini, A. Massarini, and U. Reggiani, “Transform method for calculating low-frequency shielding effectiveness of planar linear multi-layered shields,” IEEE Trans. Magn., vol. 36, no. 6, pp. 3910–3919, Nov. 2000. [5] Y. Du, T. C. Cheng, and A. S. Farag, “Principles of power-frequency magnetic field shielding with flat sheets in a source of long conductors,” IEEE Trans. Electromagn. Compat., vol. 38, no. 3, pp. 450–459, Aug. 1996. [6] W. M. Frix and G. G. Karady, “A circuital approach to estimate the magnetic field reduction of nonferrous metal shields,” IEEE Trans. Elec-tromagn. Compat., vol. 39, no. 1, pp. 24–32, Feb. 1997. [7] L. Hasselgren and J. Luomi, “Geometrical aspects of magnetic shielding at extremely low frequencies,” IEEE Trans. Electromagn. Compat., vol. 37, no. 3, pp. 409–420, Aug. 1995. [8] P. Silvester, “AC resistance and reactance of isolated rectangular conductors,” IEEE Trans. Power App. Syst., vol. PAS-86, no. 6, pp. 770–774, Jun. 1967. [9] A. Canova, G. Gruosso, and M. Repetto, “Integral methods for analysis and design of low-frequency conductive shields,” IEEE Trans. Magn., vol. 39, no. 4, pp. 2009–2017, Jul. 2003. [10] W. P. Legros and P. G. Scarpa, “Fast computation of magnetic field in rotationally symmetric structures,” IEEE Trans. Magn., vol. MAG-21, no. 6, pp. 2644–2651, Nov. 1985. [11] H. A. Haus and J. R. Melcher, Electromagnetic Fields and Energy. En-glewood Cliffs, NJ: Prentice-Hall, 1989. [12] R. F. Harrington, Time-Harmonic Electromagnetic Fields. New York: Wiley Interscience, 2001. [13] P. K. Kythe, An Introduction to Boundary Element Methods. Boca Raton, FL: CRC Press, 1995. [14] C. Pozrikidis, A Practical Guide to Boundary Element Methods with the Software Library BEMLIB. London, U.K./Boca Raton, FL: Chapman & Hall/CRC Press, 2002. [15] P. Zhou, Numerical Analysis of Electromagnetic Fields. Berlin, Germany: Springer-Verlag, 1993. [16] M. V. K. Chari and S. J. Salon, Numerical Methods in Electromagnetism. San Diego, CA: Academic Press, 2000. [17] N. Morita, N. Kumagai, and J. R. Mautz, Integral Equation Methods for Electromagnetics. Boston, MA: Artech House, 1990. [18] J. R. Carson, “Wave propagation in overhead wires with ground return,” Bell Syst. Tech. J., vol. 5, pp. 539–554, Oct. 1926. Edi Bulic was born in Pula, Croatia, in 1973. He received the B.Sc. and M.Sc. degrees in electrical engineering from the University of Ljubljana, Ljubljana, Slovenia, in 1997 and 2000, respectively. Since 1997, he has been an Assistant Professor of electromagnetics and basics in electrical engineering in the Faculty of Electrical Engineering, University of Ljubljana. His current research interests include numerical computation of the electromagnetic field. Anton R. Sinigoj was born in Ljubljana, Slovenia, in 1955. He received the B.Sc., M.Sc., and Ph.D. degrees in electrical engineering from the University of Ljubljana, Ljubljana, in 1978, 1982, and 1989, respectively. Since 1989, he has been an Assistant Professor in the Faculty of Electrical Engineering, University of Ljubljana. His current research interests include electromagnetic fields, numerical methods in electromagnetics, antennas and propagations in optical fibers, electric and magnetic shielding, and induction warming. He has authored or coauthored the books Fundamentals of Electrical Engineering and Electromagnetic Fields (Ljubljana, Slovenia: FE and FRI Publishing, 1994 and 1996, respectively). Breda Cestnik was born in Slovenia in 1967. She received the B.Sc. and M.Sc. degrees in electrical engineering from the Faculty of Electrical Engineering, University of Ljubljana, Ljubljana, Slovenia, in 1990 and 1999, respectively. Since 1990, she has been a Researcher at Milan Vidmar Electric Power Research Institute, Ljubljana. Her current research interests include effects of electric devices on the environment. Authorized licensed use limited to: UNIVERSITY OF LJUBLJANA. Downloaded on June 15, 2009 at 05:52 from IEEE Xplore. Restrictions apply.