УНИВЕРЗИТЕТ У БЕОГРАДУ ЕЛЕКТРОТЕХНИЧКИ ФАКУЛТЕТ Миодраг С. Тасић ИТЕРАТИВНО РЕШАВАЊЕ ИНТЕГРАЛНИХ ЈЕДНАЧИНА ЕЛЕКТРОМАГНЕТСКОГ ПОЉА КОРИШЋЕЊЕМ МЕТОДА ФИЗИЧКЕ ОПТИКЕ докторска дисертација Београд, 2012 UNIVERSITY OF BELGRADE SCHOOL OF ELECTRICAL ENGINEERING Miodrag S. Tasić ITERATIVE SOLUTION OF ELECTROMAGNETIC FIELD INTEGRAL EQUATIONS USING PHYSICAL OPTICS METHODS Doctoral Dissertation Belgrade, 2012 Ментор др Бранко Колунџија, редовни професор Електротехнички факултет Универзитета у Београду Чланови комисије др Антоније Ђорђевић, редовни професор Електротехнички факултет Универзитета у Београду др Братислав Миловановић, редовни професор Електронски факултет Универзитета у Нишу др Братислав Иричанин, доцент Електротехнички факултет Универзитета у Београду Датум одбране ______________ Захваљујем се: Свом ментору, професору др Бранку Колунџији Члановима комисије, професору др Антонију Ђорђевићу, професору др Братиславу Миловановићу и доценту др Братиславу Иричанину Својим колегама са електротехничког факултета у Београду Својим колегама запосленим у фирми WIPL-D у Београду Подаци о докторској дисертацији Наслов докторске дисертације: Итеративно решавање интегралних једначина електромагнетског поља коришћењем метода физичке оптике Резиме: Нумерички егзактно решавање интегралних једначина електромагнетског поља могуће је методом момената (МоМ). Када се примени на електрично велике проблеме (нпр. летелице), МоМ резултује великим системима линеарних једначина, чије решавање представља уско грло методе. Проблем је двојак: рачунарска меморија потребна за смештање матрице система је велика (сразмерна квадрату броју једначина), а рачунарско време решавања дуго (сразмерно трећем степену броја једначина, у случају решавања LU декомпозицијом). Један начин (делимичног) превазилажења овог проблема јесте паралелизација решавања система линеарних једначина. Корисницима обичних („персоналних“) рачунара на располагању су решења заснована на коришћењу вишејезгарних централних процесора или графичких процесора. На овај начин може се значајно скратити време решавања система линеарних једначина, али не и меморијско заузеће. У сваком случају, примена МоМ-а на персоналним рачунарима има јасне границе. Развијене су и бројне апроксимативне методе, између осталог засноване и на МоМ-у, за решавање електрично великих проблема. Свака од њих има своје предности и недостатке, али се ни за једну не може рећи да је универзална у погледу решавања произвољних проблема. Циљ нашег истраживања био је развијање нове апроксимативне методе, засноване на МоМ-у, за решавање интегралних једначина електромагнетског поља. Нова метода требало је да, на персоналном рачунару, омогући решавање проблема за један ред веће електричне величине од оних који се могу решити применом МоМ-а. Нова метода реализована је као итеративни поступак у којем се правац корекције одређује модификованом PO техником, а резултујући систем линеарних једначина је, за електрично велике проблеме, за више редова величина мањи од одговарајућег МоМ система. Нумерички експерименти су показали да је време потребно за добијање решења (задовољавајуће тачности) значајно краће од оног потребног за добијање МоМ решења, док је меморијско заузеће занемарљиво у односу на оно код примене МоМ-а. Електрична величина проблема који се могу анализирати персоналним рачунаром повећана је за најмање ред величине у односу на МоМ. Кључне речи: електромагнетско поље, интегралне једначине, метода момената, итеративне технике, апроксимативне методе, физичка оптика, расејачи, антене Научна област: Електротехника Ужа научна област: Електромагнетика, антене и микроталаси УДК број: 621.3 About Doctoral Dissertation Doctoral Dissertation Title: Iterative solution of electromagnetic field integral equations using physical optics method Abstract: Numerically exact solution of electromagnetic field integral equations is possible using Method of Moments (MoM). When applied to electrically large problems (e.g. airplanes), MoM results in large system of linear equations. Solving it is a bottleneck of the method, for two reasons: occupied computer memory is large (proportional to number of equations squared) and computer solution time is long (proportional to third power of number of equations, when using LU decomposition solver). A way to (partially) circumvent this problem is parallelization of solution process. Personal computer (PC) users have at their disposal software solutions based on multi-core CPUs and GPUs, which are able to significantly decrease solution time, but memory occupation stays intact. Anyway, application of MoM on PCs has clear limitations. Different approximate methods, many of them based on MoM, are developed for treating electrically large problems. Each of them has its own advantages and shortcomings, but neither is universal in respect of treating an arbitrary electromagnetic problem. The goal of our research is development of new approximate method, based on MoM, for solution of electromagnetic field integral equations. The new method should enable solution of the order of magnitude electrically larger problems than MoM. The new method is formulated as iterative procedure in which the direction of correction is determined using modified PO technique. For electrically large problems the resulting system of linear equations is much smaller than the corresponding MoM system. Numerical experiments show that (accurate enough) results can be obtained in much shorter time than with MoM, whereas memory occupancy can be neglected when compared to those of MoM. Electrical size of problems that can be analyzed using PC is increased by (at least) an order of magnitude. Keywords: electromagnetic field, integral equations, method of moments, iterative techniques, approximate methods, physical optics, scatterers, antennas Scientific area: Electrical engineering Specific scientific area: Electromagnetics, antennas and microwave UDK number: 621.3 Садржај 1. Увод…………………………………………………………………………… 1 2. Решавање интегралних једначина електромагнетског поља методом момената………………………………………………………………………. 18 2.1 Електромагнетско поље у хомогеној изотропној линеарној средини…...... 18 2.2 Електромагнетско поље у случају две хомогене средине – површински принцип еквиваленције………………………………………………….. 20 2.3 Интегралне једначине електромагнетског поља за две хомогене средине………………………………………………………………………… 27 2.4 Савршено проводан расејач у електромагнетском пољу………………....... 29 2.5 Интегралне једначине електромагнетског поља у део по део хомогеној средини……………………………………………………………………... 31 2.6 Електромагнетско поље антене на платформи……………………………... 34 2.7 Решавање (интегралних) операторских једначина методом момената…...... 34 3. Метода момената вођена физичком оптиком – нова метода за итеративно решавање интегралних једначина електромагнетског поља……………………. 38 3.1 Примена физичке оптике на савршено проводан расејач………………….. 40 3.2 Примена PDM методе на савршено проводни расејач…………………….. 42 3.3 Примена PDM методе на антену на платформи…………………………… 49 4. Детаљи реализације методе момената вођене физичком оптиком……………... 55 4.1 Геометријско моделовање и апроксимација струја МоМ функцијама базиса…………………………………………………………………………. 56 4.2 Одређивање корекционих струја…………………………………………... 62 4.3 Креирање макро функција базиса…………………………………………... 65 4.3.1 Издвајање групе плочица за МоМ зону………………………….. 68 4.3.2 Запреминско груписање плочица из PO зоне……………………. 68 4.3.3 Површинско груписање плочица из PO зоне……………………. 70 4.3.4 Креирање група МоМ функција базиса………………………….. 71 4.3.5 Конструкција Макро функција базиса коришћењем група МоМ функција базиса……………………………………………………. 72 4.4 Рачунање матрице система……………………………………………… 73 5. Примери примене методе момената вођене физичком оптиком на расејаче………………………………………………………………………... 78 5.1 Савршено проводна сфера у побудном електромагнетском пољу…… 79 5.2 Савршено проводна коцка у побудном електромагнетском пољу…… 86 5.3 Авион у побудном електромагнетском пољу………………………….. 93 5.4 Хеликоптер у побудном електромагнетском пољу……………………. 101 6. Примери примене методе момената вођене физичком оптиком на антене на великим платформама…………………………………………………….. 106 6.1 Микрострип печ антена монтирана на хеликоптер……………………. 106 6.2 Полуталасна дипол антена на авиону…………………………………... 117 6.3 Планарни низ микрострип печ антена монтираних на хеликоптер…... 121 7. Закључак………………………………………………………………………. 123 Литература………………………………………………………………………… 128 Прилог 1 – Изјава о ауторству Прилог 2 – Изјава о истоветности штампане и електронске верзије докторског рада Прилог 3 – Изјава о коришћењу Биографија аутора 1 1. Увод Тема ове дисертације суштински припада великој научној области познатој под називом електромагнетика. Назив је кованица две речи – електрицитет и магнетизам. Дисциплине које су се бавиле проучавањем особина наелектрисаних, односно намагнетисаних тела прилично независно су се развијале до средине деветнаестог века. Затим се, после неких неуспешних покушаја, 1864. појавила теорија која је трајно објединила ове две области [1]. Премда је за све време од свог настанка била изложена критикама колико и прихватањем, она је, гледајући из времена садашњег, оставила неизбрисиво наслеђе, које се пре свега огледа кроз тзв. систем Максвелових једначина (према имену свог творца), којима успешно описујемо макроскопске појаве у непокретним срединама. Врло исцрпну историју електромагнетизма написао је Витакер [2], а у новије време и Даригол (Darrigol) [3], додатно расветливши рад појединих француских истраживача у овој области, пре свега Ампера. Илустративан је и наслов његове књиге, „Електродинамика од Ампера до Ајнштајна“ – први је, поред значајних открића, заслужан и за термин „електродинамика“, док је рад другог означио прекретницу у развоју електромагнетизма – од класичне електромагнетике (описане Максвеловим једначинама) ка квантној електродинамици. У оквиру ове дисертације бавићемо се само класичном електромагнетиком, чија је теорија описана у бројним књигама: од пионирских радова Максвела [4], Хевисајда [5] (који је дао савремени облик „Максвеловим“ једначинама) и Херца [6] (који је, између осталог, први експериментално показао истинитост Максвелове претпоставке о „струјама помераја“, што је био одговор на чувено наградно питање Берлинске академије из 1879. године), преко класичног Стратоновог дела [7], до свеобухватног Џексоновог приступа [8]. Домаћим су читаоцима на располагању и универзитетски уџбеници наших професора и академика: Сурутке [9], Поповића [10] и Ђорђевића [11]. Према класичној теорији, електромагнетско понашање потиче од присуства наелектрисаних честица које се крећу или мирују. Електромагнетика описује 2 промене које овакво понашање узрокује у својој околини (при чему последњу реч треба условно схватити; електромагнетске промене се у вакуму простиру брзином светлости, тј. могу стићи веома далеко за кратко време). Пошто је „очити“ механизам њиховог деловања сила (на друге наелектрисане честице, у стању мировања или кретања), уведен је концепт електричног и магнетског поља, као дела простора у којем делују силе на наелектрисања. (Концепт поља обично се приписује Фарадеју и Максвелу, премда је полемика о томе увек актуелна [12].) Временом је електрична сила на јединично наелектрисање названа вектором електричног поља, а подужна магнетска сила на јединичну струју (под струјом подразумевамо брзину протока наелектрисања) вектором магнетске индукције. Узрочници (извори) ових вектора су наелектрисања, како она која су у „празном“ простору, тј. вакууму (слободна), тако и она која чине нераскидиви део материјала (везана). Ако се не бавимо простором између наелектрисаних честица (у самим атомима), дискретну природу наелектрисања можемо да апроксимирамо њиховом континуалном расподелом, а утицај везаних наелектрисања може се укључити у још два вектора који описују електромагнетско поље: вектор диелектричног помераја и вектор магнетског поља. Релације између вектора електричног поља и вектора диелектричног помераја, односно вектора магнетске индукције и вектора магнетског поља, називају се конститутивним релацијама. Електромагнетски систем део је природног система, па мора постојати начин да с њим размењује енергију. Још је Хевисајд у [5] стране силе које делују на наелектрисања електромагнетског система назвао побудним (impressed forces). Уобичајено је да се ове побудне силе искажу преко побудног поља (стране силе по јединичном наелектрисању) или побудних струја (које стране силе обезбеђују у делу простора). Ако поље електромагнетског таласа посматрамо довољно далеко од извора таласа, онда можемо сматрати да је повратни ефекат на примарне изворе занемарљив, па поље самог таласа, на посматраном месту, можемо прогласити побудним пољем. (Херц је 1888. године експериментално показао да се електромагнетске промене у вакууму простиру коначном брзином (светлости). Простирање тих промена названо је електромагнетским таласом.) 3 Под електромагнетском анализом подразумевамо налажење електромагнетског поља (тј. вектора поља), или неких „сродних“ величина. Максвелове једначине могу се решавати у диференцијалном или у интегралном облику. У првом случају, величина коју директно добијамо је електромагнетско поље (или помоћне функције при налажењу поља, тзв. потенцијали). Полазећи од Максвелових једначина у интегралном облику могуће је, у зависности од проблема, добити различите интегралне једначине електромагнетског поља. Њиховим решавањем обично добијамо расподелу струја. У линеарној средини у којој делују временски простопериодични извори Максвелове једначине могу се решавати у комплексном домену. Тиме се решавање у извесној мери упрошћава, а често је баш оваква анализа од интереса. У наставку текста бавићемо се искључиво питањем решавања интегралних једначина електромагнетског поља у комплексном домену. (Теорија временски простопериодичног електромагнетског поља може се пронаћи у многим књигама, нпр. [13].) Интегралне једначине електромагнетског поља спадају у линеарне операторске једначине, чији је општи облик gf =L . Познати су нам оператор L и побуда g , а треба да одредимо одзив f . Најзаступљенија општа метода за нумеричко одређивање решења f (бар у електромагнетици) је метода момената (МоМ) [14]. Творац ове методе, Харингтон, био је подстакнут методом руског инжењера Галеркина из 1915. и варијационим нумеричким методама које су се почетком шездесетих година прошлог века користиле [15]. МоМ представља уопштење поменутих техника. Апроксимативни облик решења усваја се у виду линеарне комбинације N познатих функција (базиса). Непознати тежински коефицијенти (непознате) ове линеарне комбинације одређују се из услова да „пројекција“ грешке усвојеног апроксимативног решења на сваку од N познатих (тест) функција буде једнака нули (под одређеним условима овакав поступак води ка минимизирању грешке у одређеном смислу). На тај начин добијамо систем од N линеарних једначина по непознатим тежинским коефицијентима (чиме смо „дискретизовали“ операторску једначину). Овај се систем уобичајено представља 4 у матричном облику, па говоримо о матрици система – МоМ матрици и колони слободних чланова – „члановима с десна“. Кључна питања код примене МоМ-а су: избор функција базиса и тест функција, начин израчунавања чланова МоМ матрице (и чланова с десна - у наставку се нећемо њима бавити, јер их има много мање него чланова МоМ матрице, а рачунају се на сличан начин) и начин решавања резултујућег система линеарних једначина („инверзија“ матрице). Избор функција базиса утиче на тачност решења, као и на једноставност израчунавања чланова МоМ матрице, њену величину и „стабилност“. Јасна правила за избор функција базиса не постоје, већ се може говорити о препорукама и запажањима у вези са тим [14], [16], [17]. Показује се, рецимо, да метода функционише и у појединим ситуацијама када функције базиса не испуњавају, стриктно, потребне математичке услове, нпр. када нису линеарно независне. Нумерички експерименти остају најбоља провера „оптималности“ употребљених функција базиса. Ако су функције базиса и тест функције исте, МоМ се уобичајено назива Галеркиновом методом. Ако пак операторску једначину задовољимо у одређеном броју дискретних тачака (што одговара тест функцијама које су делта импулси), говоримо о „подешавању у тачкама“ (Point Matching). Формално математички, чланови МоМ матрице одређују се пројекцијом грешке усвојеног апроксимативног решења у простор тест функција. Те се пројекције, у поступку познатом као тестирање, одређују коришћењем унутрашњег производа – операције која је, у случају решавања интегралних једначина (које су нама од интереса), интеграција. Пошто се интеграција у реалним проблемима спроводи нумерички, потребно је обратити пажњу на поједине детаље, нпр. број тачака интеграције [18]. Такође, постоји могућност апроксимације подинтегралне функције, о чему ће бити речи у наставку. Системи линеарних једначина у принципу се решавају на два начина: директно и итеративно. Код директних метода се до решења, ако постоји, долази након тачно одређеног броја рачунских операција (због матричне природе операција које се врше, обично се под једном операцијом подразумева једно множење и једно сабирање). У директне методе спадају (чувена) Гаусова 5 елиминација, која је, уједно, релативно погодна и у смислу броја потребних операција. За решавање система линеарних једначина на рачунару најчешће се користи „матрична“ верзија Гаусове елиминације – LU декомпозиција (свака Гаусова елиминација има „своју“ LU декомпозицију, али обратно не важи). За систем са 1>>N линеарних једначина израз за број потребних операција код LU декомпозиције пропорционалан је са 3N . Тада кажемо да је број потребних операција ( )3NO . Код итеративног решавања система линеарних једначина усвајамо почетно решење, које затим покушавамо да „поправимо“ понављањем одређених операција у циклусима које називамо итерацијама. За решавање интегралних једначина електромагнетског поља најчешће се користе методе засноване на тзв. Криловљевим просторима. Поједностављено гледано, код тих метода се „правац“ поправке усваја тако да се „минимизира“ грешка текућег решења. Број потребних операција по итерацији је у општем случају ( )2NO , али може бити значајно мањи када је матрица система „растресита“ (има много нула). Брзина конвергенције зависи од „стабилности“ матрице система – ову особину изражавамо преко једног реалног броја, чији је енглески назив condition number. Тај број је количник највеће и најмање сингуларне вредности матрице. Што је он већи, операције с матрицом подложније су пропагацији нумеричке грешке, а за матрицу кажемо да је мање стабилна. За стабилне матрице, са добро груписаним сингуларним вредностима, довољно тачно решење постиже се у броју итерација много мањем од броја једначина система. При томе, код симетричних матрица система, чије су све сопствене вредности позитивне, могуће је итеративни поступак организовати тако да у изразима за решење у текућој итерацији не фигуришу решења из свих осталих итерација (као у општем случају), већ из много мање (нпр. две), што скраћује израчунавања у свакој итерацији. Аспекти примене МоМ-а на које смо се осврнули су добро познати у математици, а књиге у којима су описане теме од интереса: апроксимација функција, специјалне функције, нумеричка интеграција, решавање система линеарних једначина, факторизација матрица, анализа путем сопствених 6 вредности и друге, које ћемо користити у наставку, наведене су абецедним редом [19]-[27]. Већина ових књига покрива више различитих тема (од наведених). У наставку уводимо још једно ограничење: посматрамо само решавање површинских интегралних једначина електромагнетског поља. Код примене МоМ- а за њихово решавање уобичајено је да се површ од интереса подели на више делова, па да се функције базиса дефинишу над сваким од делова. Седамдесетих и осамдесетих година прошлог века, са напредовањем рачунарске технологије (а на темељима дискретизације површинских интегралних једначина МоМ-ом), развијена су бројна решења за одређивање расподеле струја на површи металних објеката произвољног облика [28]-[36]. У овим се радовима површ обично моделује троугаоним или четвороугаоним елементима, а присутан је и модел жице, као „једнодимензионалног“ елемента. Коришћено је Галеркиново тестирање или point matching. (Уобичајени енглески назив за правац развоја електромагнетике подстакнут значајним упливом рачунара у ту област је Computational electromagnetics – CEM.) За усвојени поступак дискретизације површинских интегралних једначина и један посматрани модел, са порастом учестаности на којој се анализа врши (тј. повећањем „електричне величине проблема“, исказане у таласним дужинама), расте и број функција базиса ( N ) потребних за добијање довољно тачног решења. Рачунарски програм који треба да одреди решење има два очита задатка: да израчуна чланове МоМ система линеарних једначина („пуњење матрице“) и да, потом, тај систем реши („инверзија матрице“). Рачунарско време потребно за сваки од ова два задатка пропорционално је броју потребних операцијa. За пуњење матрице, то време је ( )2NO . Време инверзије матрице је ( )3NO код LU- декомпозиције, односно ( )2NO по итерацији (у општем случају; за „растресите“ матрице система може бити знатно краће) код неке од итеративних метода решавања система линеарних једначина (са теоријски максимално N потребних итерација). Број потребних меморијских локација је ( )2NO . Избор функција базиса утиче и на њихов потребан број, N (број непознатих). У радовима [28]-[36] ради се са тзв. функцијама базиса нижег реда – константним или линеарним. Димензије површинских елемената тада најчешће морају бити 7 реда десетине таласне дужине λ (или мање). Нпр. потребан број функција базиса дефинисан над троугаоним површинским елементима у [33] је око 300 по 2λ . Ако се ради са полиномским функцијама базиса вишег реда, тада се могу користити и површински елементи чије су димензије између једне или две таласне дужине. Неки од примера примене функција базиса вишег реда наведени су у [37]-[43], а ретроспектива оваквих решења у [44]. Поређења ради, потребан број функција базиса дефинисан над четвороугаоним површинским елементима у [37] је између 20 и 35 по 2λ , што је за ред величине боље него у [33]. На основу изложеног у претходна два пасуса, јасно је да електрична величина проблема којег можемо да решимо применом МоМ-а (а са усвојеним функцијама базиса) зависи од тога колики систем линеарних једначина, са расположивим рачунарским средствима и у задатом временском оквиру, можемо да решимо. Увешћемо претходно појам „стандардног МоМ решења“, под којим ћемо подразумевати: 1. израчунавање интеграла (при пуњењу матрице), применом неке од нумеричких формула, у довољном броју тачака (тако да интегрална сума конвергира) и без апроксимације подинтегралне функције, и 2. инверзију матрице неком од директних метода (у пракси, LU декомпозицијом), уз евентуално коришћење паралелизације (итеративне методе су, нумерички, много осетљивије на стабилност матрице система, те се, у општем случају, не могу узети као мера за тачност добијеног решења). На овај начин добијамо нумерички егзактно решење (изузимајући „грешке при заокруживању“). Одличан показатељ домета стандардног МоМ решења је рад [45]. Гледајући тамо приказане резултате, закључујемо да је, на „кућном рачунару“ са почетка друге декаде двадесет првог века, стандарно МоМ решење временски прихватљиво (под тиме овде мислимо да траје краће од једног дана; ово је, наравно, релативно) за проблеме са 100.000 непознатих. Преведено у електричну величину проблема (за функције базиса из [37]), то је око 23103 λ⋅ , или сфера полупречника око λ15 . (Из истог рада произилази да се, коришћењем паралелизације на графичким процесорима, граница подиже на читавих 300.000 непознатих.) 8 За сваки појединачни проблем, решен применом стандардног МоМ-а, можемо да поставимо питање: да ли се тај проблем могао решити (уз занемарљиву разлику у решењу) са мање непознатих? Одговор на то питање се, теоријски, може добити разлагањем матрице система применом поступка скраћеног назива SVD (Singular Value Decomposition). SVD омогућује запис матрице у виду линеарне комбинације спољашњих производа унитарних вектора [26]. Тежински коефицијенти ове суме су сингуларне вредности матрице. Одбацивањем одређеног броја чланова са најмањим сингуларним вредностима добијамо апроксимативни (компримовани) запис те матрице. Број чланова који се може одбацити без значајног губитка у тачности приказа (матрице) је мера редундантности те матрице. Очито, проблеми описани редундантним МоМ матрицама могли су бити решени и са мањим бројем непознатих (мањим бројем степени слободе). При томе, (велику) матрицу могуће је апроксимирати скраћеним записом са готово произвољном тачношћу, а коришћењем апроксимативне матрице добићемо апроксимативно решење. Тако, бар теоријски, можемо контролисати одступање апроксимативног решења од стандардног МоМ решења, односно балансирати између тачности и ефикасности решења (у смислу потребних рачунарских средстава, укључујући ту и време анализе). Пошто је SVD разлагање читаве матрице нумерички захтевнији поступак од стандардног МоМ-а, његова примена у решавању није исплатива (премда се практично користи на блоковима матрице). Али се апроксимација може вршити у подинтегралним функцијама, на добијеној МоМ матрици, при њеној инверзији, или њиховом комбинацијом. Апроксимација може сезати и до нивоа самих интегралних једначина, у потпуности или делимично заобилазећи МоМ. У наставку излагања усредсредићемо се на изналажење одговарајућег апроксимативног решења (стандардном МоМ решењу). Најочигледнији начин добијања апроксимативног МоМ решења јесте итеративно решавање МоМ система једначина, при чему се углавном користе две методе засноване на ортогонализацији Криловљевих простора: метода коњугованих градијената (Conjugate Gradients, CG) и метода бикоњугованих градијената (Biconjugate Gradients, BCG) [37]-[43]. Обе су методе са тзв. скраћеном рекурзијом: при добијању решења у текућој итерацији користе се 9 (суштински) решења из (само) претходне две. CG има смисла код симетричних матрица система чије су сопствене вредности позитивне, док је BCG модификација која ради и код несиметричних матрица система, а за симетричне се своди на CG (цена коју плаћамо за ову општост су додатне рачунске операције, намењене трансформацији проблема у симетричан). Мера блискости са стандардним МоМ решењем јесте одступање које настаје уношењем апроксимативног решења у МоМ систем једначина – грешка остатка или резидијум. Обе методе (и CG и BCG) имају две битне карактеристике: 1. брзина конвергенције је већа за стабилније матрице система, и 2. број рачунских операција је мањи код растреситијих матрица система. Дакле, пожељна је стабилна и растресита матрица система. Зато су код итеративног решавања важни поступци (preconditioner-и) који треба да трансформишу систем линеарних једначина тако да решење остане исто, али да нови систем буде погоднији за итеративно решавање [24]. Под preconditioning-ом се обично подразумева дејство на већ формирану матрицу система, али идентичан ефекат има и избор функција базиса такав да обезбеди растреситу и стабилну (тј. добро условљену – well conditioned) матрицу. Тако је, на пример, у [52] и [53] понуђено двојако решење: које може да се користи и као нова функција базиса или као preconditioner за постојеће МоМ решење. Тамо коришћене функције базиса давале су растреситу матрицу. Растресите и стабилне матрице могу се, такође, добити коришћењем Лежандрових полинома [54] („обичне“ степене функције виших редова дају нестабилне МоМ матрице). Особина која им то омогућава, међусобна ортогоналност, уједно их онемогућава да обезбеде једнакост нормалних компоненти (струја које представљају) на спојевима површинских елемената, те се морају додатно модификовати. Вероватно најпопуларнија апроксимација која се користи у последњих двадесет година позната је под скраћеницом FMM (Fast Multipole Method) [55]. FMM се заснива на апроксимацији подинтегралне функције при рачунању чланова МоМ матрице (тј. дела познатог под називом Гринова функција), на 10 начин да се један део тог интеграла учини заједничким за више различитих чланова, тј. групу (грешка која се тако чини може се контролисати). Ово за последицу има убрзано рачунање матрично-векторског производа (који је најважнија операција при итеративном решавању система линеарних једначина). Хијерархијским груписањем у више нивоа добијамо MLFMM (Multi Level FMM) [56] - [58]. Ова метода је добар пример „споја математике, физике, инжењерства и рачунарства“. Пријављена су решења проблема са више милиона непознатих, али постоје и проблеми где је овакав приступ мање ефикасан (видети нпр. [59]). Апроксимацију је могуће извршити и на израчунатој матрици система. Наиме, ако је матрица велика, а функције базиса „добро распоређене“ (тако да чланови матрице који су далеко од главне дијагонале потичу од физички удаљених функција базиса), могуће је блокове такве матрице компримовати (коришћењем SVD-а, или неког деривата тог поступка), а то затим искористити при LU декомпозицији. Типични примери су радови [60] и [61]. И FMM и технике компресије матрице садрже у себи примесе приступа који се обично назива „разлагање на области“ (domain decomposition). У наставку ћемо описати нека решења која се формално могу свести под domain decomposition. Приступ заснован на тзв. површинској еквиваленцији примењен је у радовима [62] и [63]. Структура се дели на области и за сваку се област уочавају граничне површи са другим областима. Затим се свака од области анализира независно, а утицај остатка структуре на посматрану област еквивалентира се површинским изворима на њеним граничним површима. Од великог значаја је избор граничних површи и начин одређивања еквивалентних извора (што се може учинити апроксимативно, као што је напоменуто у [62], или се могу поставити егзактни гранични услови, као у [63]). У [64] и [65] се, након геометријске декомпозиције на области, за сваку област креирају нове функције базиса. Полази се од следећег запажања: расподела струја у посматраној области биће различита за различите комбинације расподела струја у преосталим областима. Ради симулације тих комбинација, расподела струје у преосталим областима еквивалентира се фиктивним изворима изван (или на граници) посматране области. Затим се одређује одзив посматране области за 11 различите комбинације спољашњих извора, чиме се формира простор свих решења. Затим се, редукцијом тог простора, уз услов да се задржи задата тачност, формира минималан скуп нових функција базиса у посматраној области потребан да се решење представи са задатом тачношћу (чији је број мањи од оригиналног предвиђеног МоМ поставком). Нове функције базиса називају се синтетичким, и оне се даље користе у процесу МоМ дискретизације, чиме се смањује број непознатих. Може се контролисати врста и број фиктивних извори, као и начин компресије простора свих решења (нпр. у [65] примењен је SVD). У [66] се такође врши геометријска декомпозиција на области, али се користи искључиво оригинална МоМ матрица. У првом пролазу се за сваку област тражи расподела струја као да је усамљена (ради се са исечком матрице који одговара тој области, као и одговарајућим исечком из колоне с десне стране), затим се у другом пролазу за сваку област тражи решење када побуда „стиже“ из других области (за те побуде се користе струје одређене у првом пролазу), затим се решење (које је редундантно, пошто сопствене струје сваке од области залазе мало и у суседне) усвојено у виду линеарне комбинације струја (карактеристичне функције базиса) добијених у два пролаза увршћује у МоМ систем једначина, да би се затим извршило тестирање карактеристичним функцијама базиса. Резултат је систем са 2M линеарних једначина, где је M број области. Ако је 2M много мање од броја полазних једначина, додатни рад уложен у креирање карактеристичних функција базиса има смисла. Суштински, у радовима [64] и [65] (и у [66], премда то можда није очигледно) смањивање броја непознатих постигнуто је конструкцијом функција базиса прилагођених посматраној геометрији проблема. Такав приступ могао би да се сматра посебном техником, па су због тога, а и због одређених сличности са нашим истраживањем, ове технике нешто подробније описане. Преглед апроксимативних метода завршавамо тзв. асимптотским методама. Њиховом применом добијамо апроксимативно решење Максвелових једначина, које тежи тачном када учестаност анализе тежи бесконачности (тј. таласна дужина тежи нули). Уобичајено се асимптотске технике деле на оне засноване на анализи зракова (ray based) и оне засноване на анализи расподеле струја (current based). 12 Првој техници припадају геометријска оптика и геометријска теорија дифракције, а другој физичка оптика (PO) и физичка теорија дифракције. Код физичке оптике сматрамо да се поље уз површ савршено проводног објекта састоји из поља директног и рефлектованог таласа [67]. Струја на тој површи се стога одређује као двоструки векторски производ спољашњег орта на површ и вектора магнетског поља инцидентног таласа, ако „зрак“ таласа доспева до те тачке, а у супротном је та струја нула (ова апроксимација струје позната је као PO струја). Асимптотске технике одликује велика брзина рада, али и недовољна тачност. Зато су учињени бројни покушаји комбиновања (хибридизације) асимптотских метода са МоМ-ом [68] - [73]. Идеја је да се асимптотске технике употребе у делу структуре где њиховом употребом не долази до драстичног опадања тачности, а да се МоМ употреби у преосталом делу структуре. Ефикасност је боља код анализе антенских проблема, док код анализе расејача број непознатих обично може да се (само) преполови. У овој дисертацији описаћемо наше истраживање чији је циљ изналажење нове апроксимативне методе за решавање површинских интегралних једначина електромагнетског поља. Та нова метода треба да омогући решавање електрично великих проблема употребом једног кућног (персоналног) рачунара. У конкретним бројевима, желимо да проблем описан са милион МоМ непознатих решимо у року од (највише) једног дана. Полазна основа је постојеће МоМ решење описано у [74] и реализовано у софтверском пакету WIPL-D Pro [75]. Геометријско моделовање се врши билинеарним (четвороугаоним) површима - плочицама (и зарубљеним конусима - жицама), користе се полиномске функције базиса, дискретизација површинских интегралних једначина електромагнетског поља врши се Галеркиновим тестирањем, а примарни начин решавања МоМ система линеарних једначина је LU декомпозиција. Могуће је анализирати расејаче (побуђене равним електромагнетским таласом) или антене (побуђене напонским делта генератором). За посматрани модел на располагању су нам израчунати елементи оригиналне МоМ матрице, као и могућност рачунања електромагнетског поља за задату 13 расподелу струја (другим речима, те операције нису нам од интереса, јер их користимо готове). У првој фази развоја нове методе бавили смо се анализом савршено проводних расејача у пољу побудног електромагнетског таласа. Циљ анализе је одређивање расподеле струја на површи расејача, итеративним путем. За почетно решење усвојамо PO струју. У свакој наредној итерацији струје из текућег решења покушавамо да поправимо увођењем корекционих струја. Корекционе струје израчунавамо коришћењем „модификованог“ ПО-а (у свакој тачки израз за PO струје примењујемо као да зрак инцидентног таласа стиже до ње). Локално, ове струје треба до пониште магнетско поље непосредно испод површи расејача (ово поље за тачну расподелу струја треба да је нула). Глобално, то се неће десити, јер ће и корекционе струје на другим деловима површи стварати поље у свим тачкама у расејачу. Практично, очекујемо да макар релативан однос између корекционих струја на мањим (глатким) деловима површи буде „погођен“, тј. да одговара односу тачних струја. Корекционе струје (као и почетне PO струје) одређујемо у дискретним тачкама на површи, а затим их представљамо преко линеарне комбинације оригиналних (МоМ) функција базиса, чије тежинске коефицијенте израчунавамо током овог „пресликавања“. Пресликавање се обавља посебно за сваку плочицу, па није рачун(ар)ски захтевно. Корекционе струје смо, при поправци решења, покушали да употребимо на више начина. Најпре смо решење за струје у текућој итерацији одређивали додавањем корекционих струја, помножених (једним) тежинским коефицијентом, решењу за струје из претходне итерације. Тај један тежински коефицијент рачунали смо из услова да магнетско поље испод површи расејача буде минимално [76]. Затим смо решење за струје у текућој итерацији одређивали као линеарну комбинацију корекционих струја из сваке итерације (текуће и свих претходних) и почетне (PO) струје. Тежинске коефицијенте линеарне комбинације (којих је сада било за један више од броја итерација – тај „један више“ је тежински коефицијент почетне струје) и овде смо рачунали из услова да магнетско поље испод површи расејача буде минимално [77]. 14 Заправо, корекционе струје из сваке итерације, представљене у виду линеарне комбинације оригиналних функција базиса, можемо сматрати новом (великом, или макро) функцијом базиса. То се односи и на почетне (PO) струје (и оне су једна макро функција базиса). Задржавајући апроксимативно решење (за струје у текућој итерацији) у виду линеарне комбинације макро функција базиса (као у [77]), пробали смо да тежинске коефицијената те линеарне комбинације одредимо на неки другачији начин. У [78] смо апроксимативно решење (схваћено као линеарна комбинација макро функција базиса) уврстили у полазну операторску једначину, а тежинске коефицијенте одређивали смо решавањем система линеарних једначина добијеног Галеркиновим тестирањем (макро функцијама базиса). У [79] смо апроксимативно решење формално свели на линеарну комбинацију оригиналних (МоМ) функција базиса. Затим смо тежинске коефицијенте такве линеарне комбинације уврстили у оригинални (МоМ) систем једначина. Пошто уврштени тежински коефицијенти нису добијени решавањем тог система једначина, лева и десна страна у општем случају неће бити исте. Разлика десне и леве стране представља грешку „остатка“ (резидијум) посматране једначине. Укупним резидијумом (или само: резидијумом) назваћемо суму квадрата модула тих грешака (за сваку од једначина). Тежинске коефицијенте макро функција базиса добијамо из услова да резидијум буде минималан. Приступи у [77], [78] и [79] разликују се само по начини одређивања тежинских коефицијената апроксимативног решења. Они се одређују: мимимизирањем магнетског поља испод површи расејача у [77], Галеркиновим тестирањем полазне операторске једначине макро функцијама базиса у [78], и минимизирањем резидијума оригиналног (МоМ) система једначина (на претходно описани начин) у [79]. Посматрањем резулата из [77], [78] и [79], закључујемо да „минимизирање резидијума“ обезбеђује најбржу конвергенцију и највећу тачност блиског и далеког поља, па у наставку користимо само ову технику. Ради даљег убрзања конвергенције, почетне (PO) струје и корекционе струје (добијене у једној итерацији) делимо на више делова – више макро функција 15 базиса. Практично то радимо геометријском поделом структуре на области, након чега оригиналне (МоМ) функције базиса сврставамо у области по физичкој припадности (групишемо), па на основу група МоМ функција базиса формирамо макро функције базиса. Оваквом модификацијом методе из [79] добијамо текућу верзију нове итеративне апроксимативне методе за електромагнетску анализу расејача, која је описана у [80]. Метода смо назвали метода момената вођена физичком оптиком (скраћено PDM, према енглеском називу PO Driven Method of Moments). Да би се PDM метода могла употребити и за анализу антенских проблема (што сматрамо другом фазом развоја), потребне су одређене измене у односу на верзију изложену у [80]. Пре свега, корекција модификованим PO струјама има смисла само за савршено проводне затворене објекте. Стога се део структуре који не припада затвореним савршено проводним објектима (жице, генератори, диелектрични делови, отворене површи) смешта у тзв. МоМ зону, а остатак структуре у тзв. PO зону. За почетно решење усваја се расподела струја добијена МоМ анализом МоМ зоне (издвојене из остатка структуре) са оригиналном побудом. Корекционе струје у PO зони одређују се као у [80], модификованим ПО-ом, док се корекционе струје у МоМ зони одређују њеном МоМ анализом, али са побудом у виду резидијума из претходне итерације (о чему ће више речи бити при опису PDM методе у трећем поглављу). Остатак поступка сличан је оном примењеном на расејаче. Ако је број макро функција базиса много мањи од броја оригиналних (МоМ) функција базиса (тј. PDM систем има много мање једначина од одговарајућег МоМ система), отклања се потреба за решавањем великих система једначина (које даје примена МоМ-а при решавању на електрично великих проблема). Питањем ефикасности одређивања макро функција базиса бавићемо се у четвртом поглављу. Наредна три пасуса садрже математички формализам PDM-а. Мање амбициозан (или мање нестрпљив) читалац може их прескочити, јер ће у споријем ритму и са више детаља PDM бити описан у наставку. 16 У свакој итерацији PDM анализе израчунавамо корекционе струје на читавој површи структуре, које затим исказујемо линеарном комбинацијом N оригиналних (МоМ) функција базиса (са познатим, из корекционих струја израчунатим, тежинским коефицијентима). Издвајањем из ове суме NM << парцијалних сума, добијамо M макро функција базиса. Усвојену почетну расподелу струја на исти начин прерасподељујемо на макро функције базиса. Апроксимативно решење (за расподелу струја) у n -тој итерацији усвајамо у виду линеарне комбинације ( )Mn 1+ макро функција базиса из текуће и свих претходних итерација (укључујући и почетну, нулту итерацију). Апроксимативно решење у n -тој итерацији може се, очито, формално изразити линеарном комбинацијом N оригиналних (МоМ) функција базиса. У изразима за (МоМ) тежинске коефицијенте овакве суме фигуришу ( )Mn 1+ тежинских коефицијената макро функција базиса. Уврштавањем израза за N МоМ тежинских коефицијената у N оригиналних МоМ једначина, пребацивањем те суме на другу страну знака једнакости и узимањем модула те разлике, добијамо N „грешки остатка“. Уводимо резидијум – суму квадрата ових грешки – као меру одступања нашег апроксимативног решења од нумерички егзактног МоМ решења. Из услова да резидијум буде најмањи могући добијамо систем од ( )Mn 1+ линеарних једначина по ( )Mn 1+ непознатих тежинских коефицијената. Остаје нам да, решавањем овог система једначина, израчунамо те непознате тежинске коефицијенте, односно апроксимативно решење (за расподелу струја) у n -тој итерацији. PDM се може тумачити на различите начине. С једне стране, макро функције базиса подсећају на синтетичке функције базиса (из [64] и [65]) и карактеристичне функције базиса (из [66]) – креирају се на већим деловима структуре (па их је знатно мање од оригиналних функција базиса) и зависе од геометрије проблема. Али, креирају се на другачији начин од синтетичких и карактеристичних функција базиса. Ако пак PDM поредимо са CG методом (за итеративно решавање система линеарних једначина), можемо рећи да се суштински разликују по начину одређивања “правца“ поправке решења: код CG-а 17 то се ради на „математички“ начин (условљавањем ортогоналности), док је код PDM-а одређивање правца поправке „електромагнетско“ (у највећем делу структуре модификованим PO струјама, ради локалне корекције магнетског поља). Неколико речи о организацији текста у наставку. У другом поглављу су најпре приказане површинске интегралне једначине електромагнетског поља на којима се заснива електромагнетски модел који користимо у овој дисертацији, а затим и општа метода за њихово решавање - МоМ. Реч је о добро познатој теорији, чији се детаљи могу наћи нпр. у [74] и [14]. У трећем поглављу изложен је кључни допринос ове дисертације: нова итеративна метода за решавање површинских интегралних једначина електромагнетског поља - PDM. Одређивање расподеле површинских струја применом PDM-а приказано је за два случаја: савршено проводних (физички затворених) расејача и антена на (електрично великим) платформама. Описана математичка формулација независна је од коришћених површинских интегралних једначина, начина геометријског моделовања и употребљених функција базиса. У четвртом поглављу приказана је једна конкретна PDM реализација, која се ослања на МоМ поставку уграђену у нумеричко језгро софтверског пакета WIPL-D Pro [75]. Изведена је и анализа потребног рачунарског времена и меморијског заузећа. У петом поглављу приказани су примери примене PDM-а на расејачке проблеме. Одабрана су два каноничка и два реална примера. Број непознатих димензионисан је на 100.000 (број који смо у претходном излагању означили као горњу границу употребљивости „стандардног“ МоМ-а). У шестом поглављу приказани су примери примене PDM-а на антенске проблеме. Поред реалних примера са око 100.000 непознатих, приказан је и, за кућне рачунаре, екстремно велики пример од милион и по непознатих. На крају су дати закључак и списак литературе. 18 2. Решавање интегралних једначина електромагнетског поља методом момената У овом поглављу бавимо се одређивањем електромагнетског поља у део по део хомогеној (изотропној и линеарној) средини. Објаснићемо како се долази до погодних интегралних једначина који описују електромагнетско поље и како се те интегралне једначине нумерички решавају применом методе момената. 2.1 Електромагнетско поље у хомогеној изотропној линеарној средини Електромагнетско поље зависи од средине у којој постоји. За неку средину кажемо да је изотропна ако њене особине не зависе од правца вектора поља, а да је линеарна ако јој особине не зависе од интензитета поља. Средина је хомогена ако су јој особине свуда исте. За средину која поседује сва три наведена својства, Максвелове једначине, које описују макроскопско електромагнетско поље, могу се у написати у комплексном облику [74] HME ωµ−−= jrot (2.1а) EJH ωε+= jrot (2.1б) ερ=Ediv (2.1в) µρ= mdivH . (2.1г) У једначинама (2.1а-г) су подебљаним словима означени комплексни представници простопериодичних вектора: електричног ( E ) и магнетског ( H ) поља, и густина побудних запреминских струја, електричне ( J ) и магнетске ( M ). Обичним словима означени су комплексни представници простопериодичних скалара: запреминске густине електричног (ρ ) и магнетског ( mρ ) оптерећења, као и комплексна пермитивност ( ε ) и комплексна пермеабилност ( µ ) средине. 19 Кружна учестаност побудних струја означена је са ω . У изразе за ε и µ на уобичајен начин укључени су и губици, тј. ( )ωσ+ε ′′−ε′=ε j и µ′′−µ′=µ j , где је σ специфична проводност средине, а ε ′′ и µ ′′ ненегативни реални бројеви којима се описују електрични и магнетски хистерезисни губици. С обзиром да се бавимо искључиво простопериодичним побудама у хомогеној изотропној линеарној средини и сви одзиви биће простопериодични, па ћемо и све изразе писати у комплексном облику. При томе, комплексне представнике нећемо посебно обележавати (подвлачењем и сл.) јер не може доћи до забуне. Ваља напоменути да су магнетске струје и оптерећења формално уведени јер у оваквом математичком моделу омогућавају еквивалентирања појединих сложених проблема једноставнијим. Њиховом физичком остварљивошћу нећемо се бавити. Налажењем дивергенције израза (2.1а) и (2.1б) добијамо mjdiv ωρ−=M (2.2а) ωρ−= jdiv J . (2.2б) Једначина (2.2б) је диференцијални облик познате једначине континуитета, а (2.2а) је одговарајућа веза између магнетских струја и оптерећења. Пошто сходно једначинама (2.2а-б) струје и оптерећења нису независни, у изразима који следе најчешће ће фигурисати само струје - електричне и магнетске. Решавањем једначина (2.1а-г) добијају се следећи изрази за E и H [74] ( ) ( ) ( )MJrE KLZ +−= (2.3а) ( ) ( ) ( )MJrH LYK −−= , (2.3б) где је r вектор положаја тачке (у односу на усвојени координатни систем) у којој се рачуна поље, εµ=Z комплексна импеданса средине, ZY 1= комплексна адмитанса средине, а интегрални оператори L и K , уведени ради краћег записа, имају следеће изразе 20 ( ) ( ) ( ) ( ) vRRRL R v R de 111 4 2 2 γ− ∫               γ + γγ ′∇′ + γ ′ pi γ = irXrXX (2.4а) ( ) ( ) ( ) vRRK R v R de 11 4 2 2 γ− ∫               γ + γ ×′ pi γ −= irXX . (2.4б) У једначинама (2.4а-б) је εµω=γ j комплексни коефицијент простирања, v ознака области у којој се налазе извори поља ( J и M , уместо којих у изразима стоји X ), r′ вектор положаја произвољне тачке области v , rr ′−=R , а ∇′ је познати Хамилтонов оператор диференцирања који делује у тачки са вектором положаја r′ . Закључујемо да, уколико су познате (побудне) струје J и M , електромагнетско поље ( E и H ) у хомогеној изотропној линеарној средини може да се израчуна коришћењем израза (2.3а-б) и (2.4а-б). 2.2 Електромагнетско поље у случају две хомогене средине – површински принцип еквиваленције На слици 2.1 приказан је случај електромагнетског поља у две физички спрегнуте хомогене изотропне линеарне средине. Средина комплексних пермитивности 2ε и пермеабилности 2µ заузима област 2v коју ограничава површ S , док остатак простора, означеног са 1v и ограниченог површима S (унутрашња) и ∞ S (спољашња, сфера у бесконачности; у самом разматрању није од значаја пошто сматрамо да на њој нема поља, јер је бесконачно далеко од извора) заузима средина комплексних пермитивности 1ε и пермеабилности 1µ . У области 1v делују побудне струје 1J и 1M , а у области 2v побудне струје 2J и 2M . Потребно је одредити изразе за електромагнетско поље, дато векторима E и H . Једначине (2.3а-б) и (2.4а-б) у овом случају се не могу директно искористити за рачунање електромагнетског поља, с обзиром да су изведене под условом да су побудне струје у јединственој хомогеној средини. У овом потпоглављу 21 одредићемо изразе за електромагнетско поље у којима фигуришу, математички уведене, непознате површинске електричне и магнетске струје. Слика 2.1. Хомогена средина комплексних параметара 2ε , 2µ у области 2v , у којој делују побудне струје . 2J и 2M , окружена је хомогеном средином комплексних параметара 1ε , 1µ у области 1v , у којој делују побудне струје 1J и 1M , тако да их разграничава затворена површ S ; електромагнетско поље означено је векторима E и H . Електромагнетско поље у области 1v може се исказати у следећем облику [81] ( ) ( ) ( ) ( ) ( )EnHnMJrE ×+×−−+−= 1s1s11111 KLZKLZ (2.5а) ( ) ( ) ( ) ( ) ( )EnHnMJrH ×−×−−−−= 1s11s1111 LYKLYK , (2.5б) где је 111 εµ=Z , 11 1 ZY = , 1n је спољашњи јединични вектор за област 1v на површи S , а интегрални оператори sL и sK имају следеће изразе 22 ( ) ( ) ( ) ( ) SRRRL R S R de 111 4 2 s 2 s γ− ∫               γ + γγ ′∇′ + γ ′ pi γ = irXrXX (2.6а) ( ) ( ) ( ) SRRK R S R de 11 4 2 2 s γ− ∫               γ + γ ×′ pi γ −= irXX . (2.6б) У изразима за операторе sL и sK (2.6а-б) и за операторе L и K (2.4а-б), када се користе у изразима за поље (2.5а-б), комплексни коефицијент простирања је 111 j µεω=γ=γ . Интеграциона површ у (2.6а-б) је раздвојна површ S на слици 2.1. Ако изразе за операторе sL и sK упоредимо са изразима за операторе L и K видимо формалну сличност. Изрази (2.6а-б) могу се добити из израза (2.4а-б) заменом запреминске интеграције површинском ( Sv → , Sv dd → ) и заменом оператора просторног диференцирања оператором површинског диференцирања ( s∇′→∇′ ). Управо бисмо ове замене урадили када бисмо желели да помоћу израза (2.3а-б) у непрекидној хомогеној средини израчунамо електромагнетско поље површинских електричних и магнетских струја, sJ и sM . Другим речима, крајња два члана са десних страна једначина (2.5а-б) се формално могу тумачити као допринос површинских струја електромагнетском пољу у средини 1. Приметимо сада да на десним странама једнакости (2.5а-б) фигуришу искључиво побудне струје у средини 1 ( 1J и 1M ) и особине средине 1 ( 1Z и 1γ ), док се присуство средине 2 огледа само посредно, преко чланова Hn ×− 1 и En ×1 на раздвојној површи S , на које (чланове ) утиче присуство средине 2. Ако бисмо на раздвојну површ S унели фиктивне електричне и магнетске површинске струје дате са S∈×−= rHnJ 11s (2.7а) S∈×= rEnM 11s (2.7б) 23 а у области 2v средину 2 заменили средином 1 и укинули побудне струје 2J и 2M , ништа се не би променило у изразима (2.5а-б), па би и електромагнетско поље у области 1v остало исто. Слика 2.2. Електромагнетско поље у области 1v остаће исто као у случају са слике 2.1 ако се побудне струје 2J и 2M укину, у области 2v средина 2 замени средином 1, а на раздвојној површи S две области уведу фиктивне електричне и магнетске струје, 1sJ и 1sM , према изразу (2.7а-б); показује се да је поље у области 2v за овај проблем једнако нули. Ситуација након акције описане у претходној реченици приказана је на слици 2.2. Имамо хомогену средину параметара 1ε , 1µ са побудним запреминским и површинским струјама. Електромагнетско поље у овом еквивалентном проблему, које ћемо означити са ( )1E и ( )1H , може се израчунати на основу израза ( )( ) ( ) ( ) ( ) ( )1ss11ss11111111 MJMJrE KLZKLZ +−+−= (2.8а) 24 ( )( ) ( ) ( ) ( ) ( )1ss111ss1111111 MJMJrH LYKLYK −−−−= , (2.8б) формално добијених уврштавањем површинских струја дефинисаних изразима (2.7а-б) у (2.5а-б). Индекс 1 у ознакама оператора 1L , 1K , s1L и s1K указују да се у одговарајућим изразима за ове операторе користи комплексни коефицијент простирања средине 1 111 j µεω=γ=γ . Већ смо закључили да ће електромагнетско поље у области 1v бити непромењено у односу на проблем приказан на слици 2.1, тј. ( )( ) ( ) ( )( ) ( ) 11111111 ,, v∈== rrHrHrErE , (2.9а) а, као занимљив резултат, може се показати да ће електромагнетско поље у области 2v у проблему приказаном на слици 2.2 бити нула [81] ( )( ) ( )( ) 222121 ,0,0 v∈== rrHrE . (2.9б) Замена проблема са слике 2.1 проблемом са слике 2.2 уз припадајуће једначине и закључке назива се површински принцип еквиваленције. Приликом коришћења израза (2.8а-б) за одређивање електромагнетског поља у непосредној близини раздвојне површи S , нпр. за тачку са вектором положаја + 1r на слици 2.2, интеграли у изразима (2.6а-б) постају сингуларни. Да бисмо заобишли сингуларитет, интеграциону површ S делимо на део у непосредној околини сингуларитета, δS , и остатак површи, 0S , као што је приказано на слици 2.3. За површ 0S спроводимо интеграцију као (2.6а-б) сматрајући да 0→δS (резултат се обично назива Кошијева главна вредност интеграла). Овакву процедуру при рачунању (2.6а-б) ћемо назначити коришћењем ознака 0sL и 0sK . За површ δS се могу увести одређене апроксимације које воде ка коначном изразу за операторе sL и sK у тачкама у непосредној близини сингуларитета [74] 25 ( ) ( ) ( ) ( ) ( )       =⋅∇ γ − =⋅∇ γ + = − + 0 2 1 0 2 1 s 0 s s 0 s s zL zL L nXX nXX X (2.10а) ( ) ( ) ( )    =×+ =×− = − + 0 2 1 0 2 1 0 s 0 s s zK zK K nXX nXX X . (2.10б) У (2.10а-б) је n јединични вектор површи у посматраној тачки, += 0z указује на тачку мало померену од сингуларне тачке у смеру n , −= 0z указује на тачку мало померену од сингуларне тачке у смеру n− , а X замењује sJ или sM . Слика 2.3. Да бисмо избегли сингуларитет интеграциону површ S са слике 2.2 делимо на део у непосредној околини сингуларитета, δS , и остатак површи, 0S . 26 Да не би дошло до забуне, поновимо још једном: оператори sL и sK рачунају се према изразима (2.10а-б) у близини сингуларитета (тј. површинских струја), а према изразима (2.6а-б) у свим осталим тачкама. У оба случаја ознаке оператора су исте ( sL и sK ), а у зависности од ситуације користе се или изрази (2.6а-б) или изрази (2.10а-б). Слика 2.4. Електромагнетско поље у области 2v остаће исто као у случају са слике 2.1 ако се побудне струје 1J и 1M укину, у области 1v средина 1 замени средином 2, а на раздвојној површи S две области уведу фиктивне електричне и магнетске струје, 2sJ и 2sM , према изразу (2.12а-б); показује се да је поље у области 1v за овај проблем једнако нули. Поступак површинске еквиваленције можемо спровести и за средину 2, као што је илустровано на слици 2.4. Изрази за фиктивне површинске струје су 1s12s2 JHnHnJ rr −=×=×−= ∈∈ SS (2.11а) 1s122s MEnEnM rr −=×−=×= ∈∈ SS . (2.11б) 27 Видимо да се површинске струје s2J и s2M од одговарајућих струја s1J и s1M разликују само по знаку. Електромагнетско поље у овом еквивалентном проблему, које ћемо означити са ( )2E и ( )2H , може се израчунати на основу израза ( )( ) ( ) ( ) ( ) ( )1ss21ss22222222 MJMJrE KLZKLZ −++−= (2.12а) ( )( ) ( ) ( ) ( ) ( )1ss221ss2222222 MJMJrH LYKLYK ++−−= . (2.12б) У (2.12а-б) је 222 εµ=Z , 22 1 ZY = , индекс 2 у ознакама оператора 2L , 2K , s2L и s2K указују да се у одговарајућим изразима за ове операторе користи комплексни коефицијент простирања средине 2, 222 j µεω=γ=γ . Слично као у случају површинске еквиваленције за средину 1, за изразе (2.12а-б) важи ( )( ) ( ) ( )( ) ( ) 22222222 ,, v∈== rrHrHrErE (2.13а) ( )( ) ( )( ) 111212 ,0,0 v∈== rrHrE . (2.13б) 2.3 Интегралне једначине електромагнетског поља за две хомогене средине Потребно је да одредимо електромагнетско поље у области коју чине две хомогене изотропне линеарне средине приказане на слици 2.1. У претходном потпоглављу, користећи се површинским принципом еквиваленције, одредили смо изразе за електромагнетско поље у свакој од области, дате изразима (2.8а-б, 2.9а-б) и (2.12а-б, 2.13а-б). Међутим, у изразима фигуришу електричне и магнетске површинске струје, s1J и s1M , које нису познате. Да бисмо их одредили, морамо да истакнемо неки услов који електромагнетско поље мора да испуњава. У оригиналном проблему, приказаном на слици 2.1, електромагнетско поље у области 1v , дато изразом (2.5а-б), исказано је преко параметара те области и 28 тангенцијалних компоненти електромагнетског поља на граничној површи области 1v . У еквивалентном проблему, приказаном на слици 2.2 електромагнетско поље у области 1v , дато изразом (2.8а-б), исказано је преко параметара те области и фиктивних електричних и магнетских површинских струја на граничној површи области 1v . Пошто је електромагнетско поље у области 1v исто у оба случаја, очито да уведене површинске струје морају да обезбеде исте тангенцијалне компоненте вектора поља на граничној површи области 1v као у оригиналном проблему. Сличан закључак важи и за област 2v . Из израза (2.10а-б) јасно је да увођење површинских струја ствара дисконтинуитет на раздвојној површи S , али су те површинске струје резултат математичког моделовања, а у оригиналном проблему приказаном на слици 2.1 постоје само запреминске побудне струје. Стога из интегралног облика једначина (2.1а-б) следе тзв. гранични услови ( )( ) ( )( ) 0221111 =×−× ++ rr EnEn (2.14а) ( )( ) ( )( ) 0221111 =×−× ++ rr HnHn , (2.14б) који су овде написани имајући у виду слике 2.2 и 2.4 ( +1r , односно +2r , је вектор положаја произвољне тачке приљубљене уз граничну површ S са стране области 1v , односно 2v ). Изрази (2.14а-б) кажу да су тангенцијалне компоненте електричног и магнетског поља са две стране раздвојне површи S исте. Непознате површинске струје s1J и s1M могу се добити управо користећи овај услов. Уврштавањем израза за поље (2.8а-б) и (2.12а-б) у (2.14а-б), и коришћењем оператора sL и sK датих изразима (2.10а-б) добијамо ( ) ( ) ( ) ( )[ ] ( ) ( )[ ]2inc1inc11ss21ss11ss221ss111 EEnMMJJn −×=−−+× KKLZLZ (2.15а) ( ) ( ) ( ) ( )[ ] ( ) ( )[ ]2inc1inc11ss21ss11ss221ss111 HHnJJMMn −×=+++× KKLYLY , (2.15б) где су 29 ( ) ( ) ( ) 2,1,inc =+−= iKLZ iiiiii MJE (2.16а) ( ) ( ) ( ) 2,1,inc =−−= iKLY 1iiiii JMH (2.16б) познате величине (јер су познате побудне запреминске струје iJ и iM ). У (2.16а- б) индекс i означава средину, док ћемо ( )iincE , односно ( )iincH , називати побудним електричним пољем, односно побудним магнетним пољем. Изрази за електромагнетско поље у граничним условима (2.14а-б) проистекли су из Максвелових једначина, а према Теореми јединствености решења Максвелових једначина [74] гранични услови (2.14а-б) обезбеђују јединственост решења проблема приказаног на слици 2.1. Закључујемо да решавањем интегралних једначина (2.15а-б) добијамо електричне и магнетске површинске струје, s1J и s1M , на основу којих коришћењем израза (2.8а-б, 2.9а-б) и (2.12а-б, 2.13а-б) добијамо јединствено решење за електромагнетско поље које постоји у области сачињеној од две хомогене изотропне линеарне средине са познатим побудним запреминским електричним и магнетским струјама. Једначине (2.15а-б) називамо интегралним једначинама електромагнетског поља. 2.4 Савршено проводан расејач у електромагнетском пољу Реч је о специјалном случају проблема две хомогене средине приказаног на слици 2.4 и описаног интегралним једначинама поља (2.15а-б). Расејач заузима област 2v , у којој нема побудних струја ( 02 =J , 02 =M ). Област 2v испуњена је хомогеном средином врло велике проводности, која у математичком моделу тежи бесконачности, ∞→σ2 (савршено проводан материјал, тј. савршено проводан расејач). Оваква поставка значајно упрошћава једначине (2.15а-б), те ћемо таксативно навести све последице уведених претпоставки: 1. за 02 =J и 02 =M из (2.16а) следи ( ) 02inc =E и из (2.16б) следи ( ) 02inc =H 30 2. ( ) ∞−→ε⇒ωσ+ε′′−ε′=ε jj 222 3. 02222 →⇒εµ= ZZ 4. ∞→⇒= 222 1 YZY 5. ( )∞+→γ⇒µεω=γ j1j 2222 6. за ( )∞+→γ j12 из (2.6а-б) следи да ( ) 0s2 →XL и ( ) 0s2 →XK за произвољно мало али коначно R , што је испуњено за једначине (2.15а-б) 7. због последица 4 и 6 члан ( )1ss22 MLY није одређен, па се решење једначине (2.15б) може добити само ако тај члан не постоји, тј. ако је 01s =M . Коришћењем закључака 1-7, на основу (2.8а-б) и (2.16а-б) добијамо изразе за електромагнетско поље у средини 1, ( )( ) ( ) ( )1ss111inc1 JErE LZ−= (2.17а) ( )( ) ( ) ( )1ss11inc1 JHrH K−= , (2.17б) на основу (2.12а-б) добијамо изразе за електромагнетско поље у средини 2, ( )( ) 02 =rE (2.18а) ( )( ) 02 =rH . (2.18б) док се интегралне једначине електромагнетског поља (2.15а-б) своде на: ( ) ( )1inc11s0s111 EnJn ×=× LZ (2.19а) ( ) ( )1inc11s1s0s11 2 1 HnJJn ×=−× K . (2.19б) Једначина (2.19а) позната је као интегрална једначина електричног поља (EFIE - Electric Field Integral Equation), а једначина (2.19б) као интегрална једначина магнетског поља (MFIE-Magnetic Field Integral Equation). Непознате струје 1sJ 31 могу се одредити из било које од једначина (2.19а-б), али и коришћењем њихове комбинације, познате под називом комбиноване интегралне једначине поља (CFIE - Combined Field Integral Equation). Ако су побудне струје у области 1v веома далеко од области 2v , можемо сматрати да оне ни на који начин не зависе од области 2v . Тада поље ових побудних струја, ( )1incE и ( )1incH , можемо посматрати као познато побудно поље. Решавањем једначина (2.19а-б) добијамо непознате површинске струје 1sJ , а на основу тих струја и све друге величине од интереса, пре свега расејано електромагнетско поље. 2.5 Интегралне једначине електромагнетског поља у део по део хомогеној средини Део по део хомогена средина, приказана на слици 2.5, сачињена је од N области iv ( Ni ,..,1= ). У области iv , у којој делују побудне електричне и магнетске струје iJ и iM , постоји хомогена изотропна линеарна средина комплексних параметара iε , iµ . Потребно је одредити електромагнетско поље E , H . Уопштавањем површинског принципа еквиваленције закључујемо да се поље у области iv неће променити ако у свим другим областима укинемо побудне струје и испунимо их средином која је у области iv , а на све граничне површи области iv , ikS , уведемо фиктивне електричне и магнетске струје, iksJ и iksM , према изразима ikS ikik ∈×−= rHnJs (2.20а) ikS ikik ∈×= rEnMs , (2.20б) где је ikn јединични вектор нормале површи ikS усмерен од области iv . 32 Слика 2.5. Део по део хомогена средина сачињена од N области; у области iv ( Ni ,..,1= ), у којој делују побудне струје iJ и iM , постоји средина комплексних параметара iε , iµ ; гранична површ области iv и kv означена је са ikS ; област 1v није ограничена са спољашње стране; електромагнетско поље означено је векторима E и H . Еквивалентни проблем у односу на област iv приказан је на слици 2.6. Уопштавањем једначина (2.8а-б) електромагнетско поље у области iv може се изразити као ( )( ) ( ) ( ) ( ) ( )[ ]∑ = −−= N j ijiijii i KLZ 1 ssss i inc MJrErE (2.21а) ( )( ) ( ) ( ) ( ) ( )[ ]∑ = +−= N j ijiiiji i LYK 1 ssss i inc MJrHrH , (2.21б), 33 где су ( )iincE и ( )iincH рачунају према изразима (2.16а-б), оператори iLs и iKs дати су изразима (2.6а-б) или (2.10а-б), а површинске струје ijsJ и ijsM , дате изразима (2.20а-б), постоје само ако се области iv и jv граниче. Слика 2.6. Електромагнетско поље у области iv остаће исто као у случају са слике 2.5 ако се у свим областима, осим iv , укину побудне струје и средина у тој области замени средином из области iv , а на све граничне површи области iv , ikS , уведу фиктивне електричне и магнетске струје, iksJ и iksM , према изразу (2.20а-б); показује се да је поље изван области iv за овај проблем једнако нули. За проблеме приказане на сликама 2.5 и 2.6 важи ( )( ) ( )    ∉ ∈ = i ii v v r rrE rE ,0 , (2.22а) ( )( ) ( )    ∉ ∈ = i ii v v r rrH rH ,0 , . (2.22б) 34 Сада се за сваку граничну површ између две области ikS може написати по пар једначина сличних са (2.15а-б) ( ) ( )[ ] ( ) ( )[ ] ( ) ( )[ ] ikSkiik N j kjkkjkk N j ijiijiiik KLZKLZ incinc 1 ssss 1 ssss EEn MJMJn −× =         −−−× ∑∑ == (2.23а) ( ) ( )[ ] ( ) ( )[ ] ( ) ( )[ ] ikSkiik N j kjkkkjk N j ijiiijiik LYKLYK incinc 1 ssss 1 ssss HHn MJMJn −× =         +−+× ∑∑ == . (2.23б) Решавањем једначина (2.23а-б) добијамо непознате површинске струје ijsJ и ijsM , а помоћу њих и електромагнетско поље у произвољној области iv коришћењем израза (2.21а-б). 2.6 Електромагнетско поље антене на платформи Реч је о специјалном случају проблема изложеног у потпоглављу 2.5. Област 1v испуњена је вакуумом и у њој нема побудних струја. Преостале области чине платформу на којој се налази антена. Антена и платформа сачињени су од диелектричних и металних делова. У областима које чине антену постоје побудне струје. Решавањем једначина (2.23а-б) добијамо непознате површинске струје, а одатле и остале величине од интереса, пре свега електромагнетско поље зрачења антене. 2.7 Решавање (интегралних) операторских једначина методом момената Линеарна операторска једначина може се написати у компактном облику 35 gf =L , (2.24) где је L линеарни оператор, g позната векторска функција (побуда), а f непозната векторска функција коју треба одредити (одзив). Оператор врши придруживање векторске функције f векторској функцији g . За оператор се каже да је линеаран ако за сваки пар векторских функција 1f и 2f и константни скалар a важи ( ) 2121 ffff LLL +=+ , ( ) 11 ff LaaL = . (2.25) Решење за f у (2.24) се ретко може исказати преко коначног броја познатих функција и оператора (тј. добити аналитички, у тзв. затвореном облику). Нумеричко решавање (2.24) помоћу методе момената [14] почиње усвајањем апроксимативног израза за f , af , у виду линеарне комбинације N познатих, међусобно независних, векторских функција kf (познатих као функције базиса) ∑ = = N k kka 1 a ff . (2.26) Начин усвајања функција базиса је посебна тематика, али није од интереса при опису општег поступка методе момената. Непознате скаларне коефицијенте ka треба одредити тако да апроксимативно решење што је могуће мање одступа од тачног, тј. да грешка решења ffE −= a (2.27) буде најмања могућа. У реалној ситуацији када апроксимативно решење af тражимо јер немамо тачно решење f израз (2.27) не може нам користити за налажње коефицијената ka . Резидијум линеарне операторске једначине (2.24) gfR −= aL (2.28) 36 индиректни је показатељ одступања апроксимативног решења af од тачног решења f . За ff =a је резидијум једнак нули, а за ff ≠a модуо резидијума је већи од нуле. Стога очекујемо да са минимизацијом неке погодно одабране функције пропорционалне резидијуму (2.28) апроксимативно решење af тежи тачном решењу f . Једначина (2.28) важи у континуалном простору дефинисаности функције g , док је за апроксимативно решење (2.26) потребно одредити N коефицијената ka . Зато је, да бисмо минимизацијом резидијума добили коефицијенте ka , потребно претходно извршити дискретизацију једначине (2.28). Дискретизација се врши поступком познатим под називом тестирање при којем се формира N унутрашњих производа између познатих векторских функција jw и резидијума R Rw ,j , Nj ,..,1= . (2.29) Унутрашњи производ, који је симболички представљен угластим заградама , и чији је резултат обавезно скалар, део је ширег математичког концепта, а за две комплексне векторске функције, какве су у општем случају jw и R , важи jjj j Ω⋅= ∫ Ω d, * RwRw , (2.30) где * означава коњуговано-комплексну вредност, а jΩ је област дефинисаности jw у којој је 0≠jw . Заменом израза за af из (2.26) у (2.28), затим заменом тако добијеног израза за R из (2.28) у (2.29), те изједначавањем сваког од унутрашњих производа (2.29) са нулом (чиме се истиче услов ортогоналности резидијума R са сваком од функција jw ) и коришћењем особина линерних оператора датих са (2.25), добијамо систем линеарних једначина 37 gwfw ,, 1 j N k kjk La =∑ = , Nj ,..,1= . (2.31) Увођењем смена kjjk Lz fw ,= (2.32) и gw ,jjv = (2.33) систем (2.31) може се записати као j N k kjk vaz =∑ =1 , Nj ,..,1= . (2.34) Решавањем система (2.34) добијамо непознате коефицијенте ka , тиме из (2.26) апроксимативно решење af , које представља нумеричко решење једначине (2.24) коришћењем методе момената. 38 3. Метода момената вођена физичком оптиком – нова метода за итеративно решавање интегралних једначина електромагнетског поља Решавање интегралних једначина електромагнетског поља методом момената (Method of Moments, МоМ), приказано у поглављу 2, укључује неколико корака: 1. избор интегралних једначина, 2. избор функција базиса (геометријско моделовање и апроксимација струја), 3. одређивање система линеарних једначина („пуњење“ матрице), 4. решавање система линеарних једначина („инверзија“ матрице). Непознате величине које добијамо решавањем система линеарних једначина су комплексни тежински коефицијенти. Њима множимо познате функције базиса да бисмо добили тражену величину – у нашем случају површинске електричне и магнетске струје. Под решењем зато равноправно подразумевамо и тежинске коефицијенте (непознате) и расподелу струја (добијену на основу тежинских коефицијената). У погледу потребних рачунарских средстава најзахтевнији су смештање матрице система (МоМ матрице) и решавање система линеарних једначина. Матрица система састоји се од 2N комплексних бројева, а сваки од њих обично се представља са 4 бајта или 8 бајтова. Ако се определимо за 8 бајтова, за 510=N то је GB160 потребне меморије, а за 610=N сто пута више - TB16 меморијског простора заузеће ова матрица система. На почетку друге декаде двадесетпрвог века ово је на персоналним рачунарима достижно само смештањем на хард дискове, што за последицу има знатно спорији одзив него да су подаци у оперативној меморији. Методе за решавање система линеарних једначина се обично деле на директне и итеративне. Код директних метода се до решења, ако постоји, долази након тачно одређеног броја операција (једна операција укључује једно множење 39 комплексних бројева и једно сабирање комплексних бројева). Типичан представник директних метода је LU -декомпозиција, која се може посматрати као матрична верзија чувене Гаусове елиминације. Показује се да је за велики број непознатих N (шта је велики број је релативна ствар, али ћемо у овом раду сматрати да је велико 510~N , или веће) број потребних операција (самим тим и потребно рачунарско време) за решавање система линеарних једначина применом LU -декомпозицијe пропорционалан са 3N , што ћемо обележавати са ( )3NO . (Пријављене су и методе код којих је број операција ( )8,2NO и ( )376,2NO - видети нпр. [26] или [27] - али се њима овде нећемо бавити.) Итеративне методе полазе од неког усвојеног решења и затим настоје да га поправе кроз циклус истоветних корака које називамо итерацијама. За итеративно решавање система једначина насталих дискретизацијом интегралних једначина електромагнетског поља обично се користе методе коњугованих градијената (Conjugate Gradient, CG) и бикоњугованих градијената (Biconjugate Gradient, BCG). Прва метода користи се код симетричних матрица система (тј. Хермитових, када су чланови матрице комплексни бројеви), а друга нема то ограничење (видети нпр. [24] и [27]). Број операција по једној итерацији је ( )2NO , док се, под условом да поступак конвергира, до тачног решења долази у највише N итерација. Брзина конвергенције зависи од реалног броја који је количника највеће и најмање сингуларне вредности матрице система (condition number). Што је овај број већи (а он расте са порастом N ), конвергенција је у принципу спорија. Како број N расте са порастом учестаности, за електрички велике проблеме потребно је много меморије, директно решавање постаје временски (пре)дуго а код итеративног решавања брзина конвергенције опада. Другим речима, МоМ, овако реализован, постаје неефикасан. Метода момената вођена физичком оптиком (за коју ћемо, према енглеском називу Physical optics Driven Method of moments користити скраћеницу PDM) направљена је тако да заобилази коришћење великих матрица, а да ипак тежи ка МоМ решењу. Зато она задржава поједине кораке из МоМ-а, али и уводи неке нове. Тако избор интегралних једначина и функција базиса остаје непромењен. 40 Систем једначина (2.34) решавамо итеративно, тако што правац побољшања текућег решења (тј. непознатих тежинских коефицијената) процењујемо на „електромагнетски“ начин (уместо на „математички“, као код CG-а). Процена се врши на основу текућег решења. Функције базиса са процењеним вредностима тежинских коефицијената групишу се у веће целине, које ћемо звати макро функцијама базиса. Апроксимација решења се представља као линеарна комбинација макро функција базиса, а њихови тежински коефицијенти се одређују тако да решење минимално одступа од МоМ решења. Итеративни поступак додавања нових макро функција базиса у апроксимацију решења наставља се до добијања задовољавајућег решења. Техника процене вредности тежинских коефицијената слична је тзв. апроксимацији физичке оптике. Стога ћемо најпре објаснити ову апроксимацију, а затим изложити PDM методу, у две варијанте: примењену на расејаче и примењену на антене са платформом. 3.1 Примена физичке оптике на савршено проводан расејач Тело од савршено проводног материјала налази се, у вакууму, у пољу инцидентног таласа познатих вектора електричног и магнетског поља, incE и incH . Према анализи у потпоглављу 2.4 на телу ће постојати површинске електричне струје које стварају тзв. расејано поље, па такво тело називамо расејачем. Површинске струје расејача могу се добити решавањем једначина (2.19а-б) применом МоМ-а. У случају да је расејач вишеструко већи од таласне дужине оваква процедура постаје рачунарски неефикасна. Тада се може применити тзв. апроксимација физичке оптике (Physical Optics, PO). Под зраком ћемо подразумевати произвољни правац паралелан правцу простирања таласа. Површ расејача делимо на два дела: део 1S , до чије сваке тачке допире један зрак, и део 2S , до које не допиру зраци (тј. чија је свака тачка заклоњена). За прву зону кажемо да је осветљена, а за другу да је у сенци. Ово је 41 приказано на слици 3.1. Површинске електричне струје расејача сада одређујемо коришћењем израза ( ) ( )    ∈ ∈× = 2 1incPO s 0, ,2 S S r rrHn rJ . (3.1) Струје дате изразом (3.1) се уобичајено називају струјама физичке оптике. У осветљеној зони користи се израз за струје индуковане на бесконачној савршено проводној плочи изложеној дејству инцидентног таласа, а у сенци се сматра да нема струја. Слика 3.1. Савршено проводан расејач у побудном електромагнетском пољу incE , incH инцидентног таласа који се простире у правцу и смеру орта incn ; у остатку простора је ваккуум, гранична површ је S , а n је нормални спољашњи орт расејача. Према томе да ли до ње допире зрак таласа, површ расејача делимо на осветљену део 1S (светло) и део у сенци 2S (сенка); тзв. површинске струје физичке оптике добијају се помоћу израза (3.1). 42 3.2 Примена PDM методе на савршено проводни расејач Проблем расејача приказан је на слици 3.1 и објашњен у претходном потпоглављу. Да бисмо добили површинске струје расејача решавамо једну од једначина (2.19а-б). Та једначина се може написати у операторском облику (2.24), где се ознака f односи на површинске струје расејача. Нумеричка апроксимација за f дата је изразом (2.26), где су kf усвојене познате векторске функције базиса, којих има N . Циљ је да одредимо тежинске коефицијенте ka у изразу (2.26), а тиме и решење проблема. Струје дате изразом (3.1) представљају једно апроксимативно решење проблема. Оне се могу приказати преко функција базиса, у облику (2.26) ( ) )0( a )0( 1 0PO s fFfJ ==≅ ∑ = N k kka . (3.2) Овде индекс нула означава да је реч о почетној (или нултој) итерацији, ( )0F ћемо називати почетном расподелом струја, а ( )0ka су почетне вредности тежинских коефицијената функција базиса. Истовремено, узећемо да је ова процена, добијена физичком оптиком, почетно решење нашег проблема, ( )0af . Ово почетно решење код реалних проблема увек одступа од (нумерички) тачног решења, добијеног МоМ-ом. Да бисмо закључили како да побољшамо почетно решење посматраћемо проблем са слике 3.2, еквивалентан ономе са слике 3.1 у односу на област ван расејача. У еквивалентном проблему поље у области расејача не постоји (изузев близу резонантних учестаности, у случају када се користи само једна од једначина (2.19а-б), а не њихова комбинација; овим се проблемом нећемо бавити, и сматраћемо да смо далеко од резонантних учестаности). 43 Слика 3.2. Проблем еквивалентан ономе са слике 3.1 у односу на област v . Утицај расејача замењен је увођењем површинских струја sJ добијених решавањем (2.19а-б). Електромагнетско поље у области v исто је као у проблему са слике 3.1, док је електромагнетско поље у области расејача, која у еквивалентном проблему има средину као и област v , једнако нули. За апроксимативну расподелу струја у расејачу ће постојати неко поље грешке. Да бисмо у текућој итерацији, n , поправили решење из претходне итерације, ( )1-anf , покушаћемо да на граничну површ додамо корекционе струје које ће умањити поље грешке и тиме нас евентуално приближити тачном решењу. Посматрајмо зато магнетско поље које у расејачу постоји у претходној итерацији, а које је приказано на слици 3.3. Употребићемо магнетско поље на површи −S , непосредно испод граничне површи расејача. Оно је једнако суми инцидентног магнетског поља incH и расејаног магнетског поља ( )( )1-anfH , ( ) ( ) inc)1(a1 HfHH += −− nn . (3.3) 44 Слика 3.3. Уколико се површинске струје sJ на слици 3.2 замене њиховом апроксимативном вредношћу ( )1- a nf , поље у расејачу биће различито од нуле. Поправку решења вршимо на основу вредности магнетског поља на површи −S , непосредно испод граничне површи расејача S . Расејано магнетско поље рачуна се коришћењем израза ( ) )1(aj-2)1(a)1(a 21dej1)( −β−− ×+β+×= ∫ nS R R nn S R R ' fnirffH (3.4) где је 002 µεpi=β f , r је вектор положаја тачке на површи −S (где рачунамо магнетско поље), r′ је вектор положаја тачке на површи S (где су површински извори магнетског поља), rr ′−=R и ( ) rrrri ′−′−=R . Локално, тангенцијална компонента поља ( )1−nH може се поништити додавањем површинских струја ( ) ( ) −−×=∆ Snn 2 1s HnJ . (3.5) 45 Претходно следи из израза за магнетско поље које се у непосредној близини струјног елемента може апроксимирати магнетским пољем бесконачног струјног листа. Ове струје могу се представити као линеарна комбинација функција базиса, ( ) ( )∑ = =≅∆ N k k n k nn a 1 )( s fFJ , (3.6) где ћемо )(nF називати корекционим струјама у n -тој итерацији ( 0>n ), а ( )nka су корекционе вредности тежинских коефицијената функција базиса у n -тој итерацији. Струје дате изразом (3.5) поништавају тангенцијално магнетско поље непосредно испод тачке на површи где су додате, али ову поправку кваре додатне струје постављене у свим осталим тачкама површи S . Нумерички експерименти показују да је при простом додавању корекционих струја )(nF постојећем решењу ( )1- a nf поље грешке у расејачу веће него пре њиховом додавања. Другим речима, овако организован итеративни поступак дивергира. Ако се, пак, корекционе струја )(nF пре додавања постојећем решењу ( )1-anf помноже тежинским коефицијентом погодно одабране вредности, итеративни поступак конвергира, али веома споро. Да бисмо убрзали конвергенцију, корекционе струје разложићемо на M делова, где сваки део представља једну макро функцију базиса. На тај начин ћемо корекционе струје представити као суму макро функција базиса ∑ = = M l n l n 1 )()( FF . (3.7) Макро функције базиса могу се представити као линеарна комбинација свих МоМ функција базиса ( ) Mla N k k n kl n l ,...,1, 1 )( == ∑ = fF , (3.8) 46 где ( )nlF означава l -ту макро функцију базиса у n -тој итерацији. Заменом израза за )(nF из (3.6) и израза за ( )nlF из (3.8) у (3.7), након измене редоследа сумирања са десне стране, директним поређењем са левом страном добијамо Nkaa M l n kl n k ,..,1, 1 )()( == ∑ = . (3.9) На основу претходног израза закључујемо да се тежински коефицијенти макро функција базиса, ( )nkla , добијају поделом тежинских коефицијената МоМ функција базиса, ( )nka . Ова подела зависи од коришћених МоМ функција базиса, о чему ће бити речи у наредном поглављу. Апроксимативно решење у n -тој итерацији, ( )naf , усвојићемо као линеарну комбинацију свих макро функција базиса, уведених у претходним и у текућој итерацији ( )∑∑ = = = n i M l i l n il n c 0 1 )()( a Ff , (3.10) где су )(nilc непознати тежински коефицијенти које треба одредити. Након уврштавања израза (3.8) за макро функцију базиса, (3.10) постаје ∑∑ ∑ = = =       = n i M l N k k i kl n il n ac 0 1 1 )()()( a ff , (3.11) да би се након измене редоследа сума у (3.11) ( )naf могло записати као ∑ ∑∑ = = =       = N k k n i M l i kl n il n ac 1 0 1 )()()( a ff . (3.12) Увођењем смене ∑∑ = = = n i M l i kl n il n k acA 0 1 )()()( (3.13) 47 израз (3.12) постаје ( ) 1 )( a ∑ = = N k k n k n A ff . (3.14) Израз (3.14) је формално сличан апроксимативном МоМ решењу (2.26). Међутим, док су у изразу (2.26) непознате које треба одредити N тежинских коефицијената ka , у (3.14) непознате су ( )Mn 1+ тежинских коефицијената )(nilc , који су у коефицијенте ( )nkA уграђене путем израза (3.13). Да бисмо апроксимативно решење (3.14) што више приближили МоМ решењу (2.26) потребно је да коефицијенти ( )nkA буду одабрани тако да се минимално разликују од коефицијената ka . Пошто коефицијенте ka , који се добијају се решавањем система линеарних једначина (2.34), не знамо, преостаје нам да коефицијентима ( )nkA што боље задовољимо систем (2.34). Зато за сваку од једначина система (2.34) конструишемо меру одступања, ( ) ( ) NjvAzR j N k n kjk n j ,...,1 , 1 =−= ∑ = (3.15) и дефинишемо средњи квадратни остатак (резидијум) ∑ = = N j n j n R N R 1 2)()( 1 (3.16) као меру у којој апроксимативно решење ( )naf , дато изразом (3.14), не задовољава систем (2.34). Заменом израза за ( )njR из (3.15) у (3.16) добијамо ( )∑∑ = = −= N j j N k n kjk n vAz N R 1 2 1 )( 1 , (3.17) па након уврштавања израза (3.13) за ( )nkA у (3.17) )(nR постаје 48 ∑∑ ∑∑ = = = = −= N j j N k n i M l i kl n iljk n vacz N R 1 2 1 0 1 )()()( 1 . (3.18) Након измене у редоследа сумирања, израз (3.18) се може написати у облику ∑∑ ∑∑ = = == −= N j j n i N k i kljk M l n il n vazc N R 1 2 0 1 )( 1 )()( 1 . (3.19) Увођењем нове ознаке за унутрашњу суму у (3.19), ( ) ∑ = = N k i kljk i jl azZ 1 )( , (3.20) добијамо израз за резидијум )(nR у функцији непознатих тежинских коефицијената )(nilc ( ) ∑∑∑ = = = −= N j n i j M l i jl n il n vZc N R 1 2 0 1 )()(1 . (3.21) Очекујемо да се смањивањем резидијума )(nR апроксимативно решење ( )naf приближава МоМ решењу. Зато као услов налажења тежинских коефицијената )(n ilc истичемо минимизирање резидијума )(nR , тј. Mmnk c R n km n ,...,1,,...,0,0)( )( === ∂ ∂ . (3.22) Резултат примене услова (3.22) је систем једначина ( ) ( ) ( ) ( ) Mmnk ZvZZc N j k jmj n i M l N j k jm i jl n il ,...,1,,...,0 1 * 0 1 1 * == =         ∑∑∑ ∑ == = = . (3.23) Решавањем система линеарних једначина (3.23) добијамо тежинске коефицијенате )(nilc , самим тим и апроксимативно решење у n -тој итерацији, ( )n af , 49 дато изразом (3.14). Затим можемо израчунати резидијум на основу израза (3.21). Ако резидијум није довољно мали, можемо поновити поступак додавањем нових макро функција базиса и одређивањем новог апроксимативног решења. Уколико је ( ) NMn <<+1 систем (3.23) је далеко економичнији за решавање од система (2.34). Питањем укупног броја потребних операција у оба случаја бавућемо се у наредним поглављима. Пошто за почетно решење усвајамо оно добијено физичком оптиком, дато изразом (3.2), корекционе струје у свакој итерацији одређујемо коришћењем израза (3.5), који је апроксимација физичке оптике у осветљеној зони, и како коначно тежимо МоМ решењу, описана метода добила је назив метода момената вођена физичком оптиком. 3.3 Примена PDM методе на антену на платформи Проблем антене на платформи укратко је изложен у потпоглављу 2.6. Антенски проблем подразумева постојање побудних струја, а дозволићемо и постојање диелектричних делова. Решавањем интегралних једначина (2.23а-б) добијамо еквивалентне површинске струје, а одатле остале величине од интереса, нпр. електромагнетско поље у зони зрачења. Једначине (2.23а-б) приказујемо у операторском облику (2.24), а коначни циљ је „заобилазним“ путем – уместо решавањем система једначина (2.34) - доћи до коефицијената ka у изразу (2.26). Корекционе струје дате изразом (3.5) и (3.6) имају смисла само на раздвојној површи савршено проводног материјала. Зато ћемо површ проблема поделити на два дела: зону напајања и границе диелектричних области доделићемо тзв. МоМ зони, а остатак тзв. PO зони, као што је приказано на слици 3.4. Практично, у МоМ зону улази антена са ближом околином, а у PO зону (најчешће драматично већа) платформа. 50 Слика 3.4. Проблем антене на платформи делимо на МоМ зону, која обухвата зону напајања и диелектричне делове, и PO зону, која обухвата делове сачињене од савршеног проводника, ван зоне напајања. Околна средина је вакуум. Поправку решења у PO зони вршимо на основу вредности магнетског поља на површи −S , непосредно испод граничне површи PO зоне и вакуума. У складу са овом поделом, и корекционе струје у n -тој итерацији ћемо поделити на МоМ део и PO део, ( ) ( ) ( ) ( )nnN Nk k n k N k k n k n aa POMoM 11 )( MoM MoM FFffF +=+= ∑∑ +== . (3.24) Изразом (3.24) исказујемо да првих MoMN функција базиса припадају МоМ зони, а да преостале функције базиса припадају PO зони. У свакој итерацији корекционе струје у МоМ зони чине једну функцију базиса, коју ћемо нумерисати као нулту у тој итерацији, тј. ( ) ( )nn 0MoM FF = , (3.25), 51 док ћемо корекционе струје ( )nPOF , као и код расејача, разлажити на M макро функција базиса ( ) ∑ = = M l n l n 1 )( PO FF , (3.26) које се могу се представити као линеарна комбинација свих функција базиса из PO зоне ( ) Mla N Nk k n kl n l ,...,1, 1 )( MoM == ∑ += fF , (3.27) где ( )nlF означава l -ту макро функцију базиса у n -тој итерацији и где, као и код расејача, важи NNkaa M l n kl n k ,..,1, MoM 1 )()( +== ∑ = . (3.28) тј. тежински коефицијенти макро функција базиса, ( )nkla , добијају се поделом тежинских коефицијената МоМ функција базиса, ( )nka . Формалним издвајањем МоМ зоне из система једначина (2.34) добијамо делимични систем једначина j N k kjk vaz =∑ = MoM 1 , MoM,..,1 Nj = , (3.29) где су чланови jkz и jv познати, дати изразима (2.32) и (2.33). Коефицијенти ka добијени решавањем делимичног система једначина (3.29) не представљају тачно решење за првих MoMN коефицијената у систему једначина (2.34), јер је занемарен утицај функција базиса из PO зоне. Ипак, они могу послужити као добро почетно решење ( ) ( ) ( ) ( )0 0 0 MoM 1 00 a MoM FFff === ∑ = N k kka , (3.30) 52 где су коефицијенти ka добијени решавањем делимичног система једначина (3.20) овде погодно означени са ( )0ka . Апроксимативно решење у n -тој итерацији, ( )naf , усвојићемо као линеарну комбинацију свих макро функција базиса, уведених у итерацијама почевши од прве до текуће ( )∑∑ = = = n i M l i l n il n c 1 0 )()( a Ff , (3.31) где су )(nilc непознати тежински коефицијенти које треба одредити. Израз за апроксимативно решење код антена (3.31) се од одговарајућег израза код расејача (3.10) разликује само у почетним вредностима индекса i и l . Наиме, у нултој итерацији имамо само једну макро функцију базиса ( )00F , дату изразом (3.30), која се, као што ћемо видети, користи неизмењена у првој итерацији, па бројач итерација i сада креће од 1. С друге стране, да бисмо потцртали увођење додатне макро функције базиса у МоМ зони, нумерисали смо је као нулту, те бројач функција базиса l сада креће од 0. Како су промене израза (3.31) у односу на израз (3.10) формалне природе, тежинске коефицијената )(nilc у (3.31) добијамо истоветно као оне у (3.10), тј. минимизирањем резидијума ( ) ∑∑∑ = = = −= N j n i j M l i jl n il n vZc N R 1 2 1 0 )()(1 , (3.32) а што се своди на решавање система једначина ( ) ( ) ( ) ( ) Mmnk ZvZZc N j k jmj n i M l N j k jm i jl n il ,...,0,,...,1 1 * 1 0 1 * == =         ∑∑∑ ∑ == = = . (3.33) Видимо сличност систему (3.23) – уз већ изречену опаску о индексима. Остаје нам да прецизирамо добијање корекционих струја. 53 Корекционе струје у n -тој итерацији у PO зони одређујемо као код расејача: најпре коришћењем израза (3.4) рачунамо магнетско поље које на површи −S са слике 3.4 постоји у претходној итерацији, затим на основу израза (3.5) рачунамо површинске струје ( )nsJ∆ , које коначно представљамо преко функција базиса у PO зони, ( ) ( ) ( )∑ += =≅∆ N Nk k n k nn a 1 POs MoM fFJ . (3.34) Корекционе струје у n -тој ( 1>n ) итерацији у МоМ зони одређујемо полазећи од резидијума за првих MoMN једначина из система (2.34), у 1−n -ој итерацији ( ) MoM 1 1 0 )()1(1 ,...,1 , NjvZcR n i j M l i jl n il n j =−= ∑∑ − = = −− . (3.35) Уврштавањем ових резидијума на десне стране једначина (3.29), уместо слободних чланова jv , добијамо систем једначина ( ) ( ) MoM 1 1 ,..,1, MoM NjRaz nj N k n kjk == − = ∑ . (3.36) Да резидијуми дати изразом (3.35) потичу само од функција базиса из МоМ зоне, додавањем коефицијената )(nka , добијених решавањем система (3.36), одговарајућим коефицијентима добијеним у претходној итерацији остварили бисмо потпуну поправку решења. То, наравно, није случај, јер резидијуми зависе и од функција базиса из PO зоне. Ипак, овако добијени коефицијенти )(nka могу послужити за креирање нових корекционих струја, односно нове макро функције базиса у МоМ зони, ( ) ( ) ( ) 1, MoM 1 0MoM >== ∑ = na N k k n k nn fFF . (3.37) У првој итерацији као макро функцију базиса користимо почетно решење (3.30), добијено решавањем делимичног проблема, 54 ( ) ( ) ( ) ( )∑ = === MoM 1 00 a 1 0 1 MoM N k kka ffFF . (3.38) Поновимо укратко описани итеративни поступак. У свакој итерацији додајемо нове функције базиса, у МоМ зони дате изразима (3.37) или (3.38), а у PO зони дате изразом (3.34). Апроксимативно решење представљамо у облику (3.31), а непознате тежинске коефицијената )(nilc добијамо решавањем система једначина (3.33). Итеративни поступак је завршен када резидијум дат изразом (3.32) постане мањи од задатог. 55 4. Детаљи реализације методе момената вођене физичком оптиком Mетода момената вођена физичком оптиком (скраћено PDM), описана у претходном поглављу, може се схватити као надградња методе момената (МоМ- а). Она преузима читаву поставку из МоМ-а, додајући затим два своја концепта: макро функције базиса и PDM матрицу система. Методологија изложена у претходном поглављу је општа, али детаљи реализације зависе од усвојене МоМ реализације. Зато у овом поглављу најпре уводимо МоМ поставку, а затим се бавимо конструисањем макро функција базиса и израчунавањем PDM матрице система. Кад говоримо о МоМ поставци, у виду имамо коришћене функције базиса (које ћемо овде прецизности ради звати МоМ функцијама базиса) и геометријске елементе над којима су оне дефинисане. У овој тези усвојено је познато решење детаљно описано у [74], а биће изложен део неопходан за разумевање решења примењених у PDM-у. Подсетимо се сада смисла макро функција базиса. Показује се да итеративни поступак у коме се корекционе струје на читавој површи множе једним тежинским коефицијентом (пре него што се додају постојећим струјама) споро конвергира. Да бисмо убрзали конвергенцију, потребно је корекционе струје поделити на више целина – макро функција базиса - које се множе различитим тежинским коефицијентима. У граничном случају број макро функцији базиса једнак је броју МоМ функција базиса, а итеративни поступак своди се на једну итерацију у којој добијамо МоМ решење. У реалном случају број макро функција базиса је много већи од 1, али много мањи од броја МоМ функција базиса. До макро функција базиса долазимо поделом корекционих струја према изразу (3.24). Претходно је потребно да дискретне површинске струје добијене коришћењем израза (3.5) представимо линеарном комбинацијом МоМ функција базиса, као у изразу (3.34), како бисмо добили корекционе коефицијенте ( )nka у PO зони. 56 Систем једначина (3.23) код расејачких, односно (3.33) код антенских проблема, може се компактно записати у матричном облику. Ако систем садржи PDMN једначина, са десне стране матричне једначине фигурисање квадратна матрица познатих чланова, реда PDMN , коју називамо PDM матрицом (система). Од интереса је утврдити колико је операција потребно за њено рачунање и колико је меморијског простора потребно за њено смештање. 4.1 Геометријско моделовање и апроксимација струја МоМ функцијама базиса Под геометријским моделовањем подразумевамо математички опис структуре чији утицај на електромагнетско поље анализирамо. Пошто решавање проблема, изложених у другом поглављу, у крајњој линији сводимо на одређивање еквивалентних површинских струја, структуру ћемо моделовати граничним површима између различитих средина. Као површинске елементе користићемо тзв. билинеарне површи. Слика 4.1. Билинеарна површ облика четвороугала описана је векторима положаја своја четири темена, према изразу 4.1; у општем случају, реч је о неравној површи. 57 Билинеарна површ, приказана на слици 4.1, је у општем случају нераван четвороугао чија се површ може описати параметарском једначином 1,1,])1(1][)1(1[ 4 1),( 2 1 2 1 ≤≤−−+−+= ∑∑ = = spspsp i j ji ijrr . (4.1) Са слике 4.1 види се да су вектори положаја задати у глобалном координатном систему (важећем за читаву структуру, тј. и све остале билинеарне површи), док је површ посматраног елемента описана у локалном sp, координатном систему (који важи само за посматрани површински елемент). Преуређивањем израза (4.1) добијамо другачији приказ билинеарне површи 1,1,),( psspc ≤≤−+++= sppsspsp rrrrr , (4.2) где су cr , pr , sr и psr константни вектори, дати изразима ( )22211211c 4 1 rrrrr +++= (4.3а) ( )22211211p 4 1 rrrrr ++−−= (4.3б) ( )22211211s 4 1 rrrrr +−+−= (4.3в) ( )22211211ps 4 1 rrrrr +−−= . (4.3г) Ако је p константно, израз (4.2) се своди на 11,)( ss0 ≤≤−+= sss arr , (4.4) где су s0r и sa константни вектори, дати изразима pcs0 rrr p+= (4.5а) 58 psss rra p+= . (4.5б) Пошто израз (4.4) описује праву у функцији параметра s , закључујемо да су све s координатне линије праве. На сличан начин, за константно s , израз (4.2) своди се на 11,)( pp0 ≤≤−+= ppp arr , (4.6) где су p0r и pa константни вектори, дати изразима scp0 rrr s+= (4.7а) pspp rra s+= . (4.7б) Наравно, на основу израза (4.6), и закључак је сличан: све p координатне линије су праве. Вектори sa и pa представљају правце координатних линија. Елементарна површина у локалном sp, координатном систему, према слици 4.2, може се изразити у облику spS ddd sp aa ×= . (4.8) Слика 4.2. Елементарна површина у локалном sp, координатном систему једнака је интензитету векторског производа елементарних вектора положаја. 59 Електричне и магнетске струје које постоје на билинеарној површи разлажу се на p компоненту и s компоненту у локалном sp, координатном систему. Функције базиса су исте за електричне и магнетске струје, а p компонента може се добити из s компоненте формалном заменом p и s координата и векторских компоненти. Израз за s компоненту електричне струје на једној билинеарној површи узимамо у облику ∑ ∑ − = = = 1 0 0 s p s ),(),( n i n j ijij spasp PJ , (4.9) где је ( )          >− =+ =− × = − 11 12/)1( 02/)1( ),( 22sp s jss js js psp j i ij aa aP , (4.10) а pn и sn тзв. редови апроксимације по p и по s координати. Функције базиса за 1>j једнаке су нули за 1±=s . Њима описане струје стога се не „преливају“, преко заједничке ивице, у суседну билинеарну површ. У апроксимацији датој изразом (2.26), свака од ових функција базиса ijP представља једну функцију базиса kf , а сваки од одговарајућих коефицијената ija један коефицијент ka . Функције базиса за 0=j различите су од нуле за за 1−=s , а функције базиса за 1=j различите су од нуле за за 1=s . Њима описане струје „преливају“ се, преко заједничке ивице, у суседну билинеарну површ. Једна заједничка ивица две суседне билинеарне површи приказана је слици 4.3. 60 Слика 4.3. Дуж заједничке ивице две билинеарне површи постављена је p -оса, док свака површ има своју s -осу. Обе површи леже у равни цртежа. Приказане су вођице све три осе, као и нормалне компоненте вођица обе s -осе Према изразу (4.10), за сваку од ових површи функције 0iP су различите од нуле на заједничкој ивици. Изрази за ове функције на површи 1, односно 2, су: ( ) ( ) ( ) ( ) 2 1 1 1 sp 1 s1 0 spii − × = aa aP (4.11а) ( ) ( ) ( ) ( ) 2 1 2 2 sp 2 s2 0 spii − × = aa aP . (4.11б) Пошто је дуж заједничке ивице на слици 4.3 ( ) ( ) 121 −== ss , интензитети нормалних компоненти ове две функције дуж заједничке ивице су ( ) ( ) ( ) i i p1 sp 1 norms1 norm0 aa a P × = (4.12а) ( ) ( ) ( ) i i p2 sp 2 norms2 norm0 aa a P × = . (4.12б) 61 На основу слике 4.3 можемо написати изразе за интензитете компоненти вектора ( )1 sa и ( )2 sa нормалних на заједничку ивицу, ( ) ( ) p 1 sp1 norms a aa a × = (4.13а) ( ) ( ) p 2 sp2 norms a aa a × = . (4.13б) Уврштавањем израза за ( )1normsa и ( )2normsa из (4.13а-б), изрази (4.12а-б) постају ( ) i i p p 1 norm0 1 a P = (4.14а) ( ) i i p p 2 norm0 1 a P = . (4.14б) Из израза (4.14а-б) видимо да су нормалне компоненте функција ( )10iP и ( )20iP истих интензитета, а са слике 4.3 види се да су оне супротних смерова. Ако су коефицијенти 0ia уз ове функције у изразу (4.9) разликују само по знаку, тада ће на заједничкој ивици билинеарних површи са слике 4.3 бити аутоматски задовољена једначина континуитета за стационарно струјно поље. Зато се функције ( )10iP и ( )20iP комбинују у јединствену функцију базиса, чији је израз према референтном смеру ( )1sa : ( ) ( ) ( )2 0 1 0 2,1 0 iii PPD −= . (4.15) У апроксимацији датој изразом (2.26), свака од функција ( )kji ,0D , где су j и k индекси две произвољне билинеарне површи које имају заједничку ивицу, представља једну функцију базиса kf , и свакој је придружен један коефицијент ka . 62 Због начина конструкције, функције базиса ( )kji ,0D називамо дублетима, док функције базиса ijP ( 1>j ) називамо синглетима. Поред билинеарних површи, које ћемо у наставку равноправно називати и плочицама, у МоМ решењу на које се овде надовезујемо као градивни елементи користе се и жице. Жице се геометријски описују правим цилиндрима чији се попречни пресек континуално мења између крајњих тачака осе цилиндра. Жице ћемо пре свега користити у антенским проблемима, ради прикључивања побудних генератора. Жице обавезно припадају МоМ области структуре, па њихов тачан геометријски опис и апроксимација струја дуж њих нису од интереса при опису PDM-а (потребне величине узимају се из постојећег МоМ решења). 4.2 Одређивање корекционих струја Под одређивањем корекционих струја подразумевамо налажење одговарајућих (корекционих) коефицијената уз МоМ функције базиса. Ови коефицијенти се у МоМ зони антенских проблема добијају директно, решавањем система једначина (3.36). Међутим, код расејачких проблема (код којих МоМ зона не постоји) и у PO зони антенских проблема корекциони коефицијенти одређују се посредно, на основу корекционих струја добијених у дискретним тачкама на површи, према изразима (3.4-5). Стога одређивање корекционих коефицијената у PO зони захтева додатне операције које нису једнозначне, а које ћемо размотрити у наставку. Магнетско поље дато изразом (3.4) рачуна се непосредно испод површи посматраног тела. Одређивање тачака испод површи носило би са собом извесна питања – на којој тачно „дубини“ те тачке треба да буду, да ли се налазимо на веома танком делу тела и слично. Срећом, ово магнетско поље може се израчунати и користећи тачке на самој површи тела. Наиме, интеграл у изразу (3.4) представља главну Кошијеву вредност, односно искључује део површи у непосредној близини тачке у којој се поље рачуна. Зато се његова вредност практично не мења ако тачку рачунања „подигнемо“ на саму површ. Други део израза за поље, векторски производ са десне стране знака једнакости у изразу 63 (3.4), подразумева да је тачка рачунања непосредно испод површи, али не укључује и њен тачан положај. Зато се магнетско поље испод површи тела, дато изразом (3.4), може израчунати за тачке задате на самој површи тела. То су исте оне тачке на површи за које се рачунају корекционе струје дате изразом (3.5). Површинске струје чије вредности знамо у коначном броју тачака на површи тела S треба, према изразу (3.6), апроксимирати линеарном комбинацијом функција базиса (тежински коефицијенти ове линеарне комбинације су корекциони коефицијенти које тражимо). Резидијум ове апроксимације у свакој од задатих тачака на површи тела је: ( )( ) ( ) ( ) . 1 ss ∑ = −∆=∆ N k k n k nn aR fJJ (4.16) Непознате коефицијенте ( )nka одредићемо тако да минимизирамо квадратну грешку резидијума ( )( ) ( )( )∫ ∆=∆ S nn SRE d 2 ss JJ , (4.17) тј. из услова ( )( ) ( ) Nja E n j n ,...,1,0s == ∂ ∆∂ J , (4.18) чија примена резултује системом једначина ( ) ( ) NjSSa S n j N k S kj n k ,...,1,dd s 1 =∆⋅=⋅ ∫∑ ∫ = Jfff . (4.19) Због природе функција базиса kf описаних у претходном потпоглављу (које су дефинисане или само на једној билинеарној површи, у случају синглета, или на две суседне билинеарне површи, у случају дублета) систем једначина (4.19) биће растресит (тј. са пуно нултих чланова), али би и поред тога његово решавање за велики број функција базиса било рачунарски захтевно. Зато се систем једначина 64 облика (4.19) поставља за сваку билинеарну површ независно, тј. за билинеарну површ i решавамо систем: ( ) ( ) i S n j N k i S kj n k NjSSa i i i ,...,1,dd s 1 =∆⋅=⋅ ∫∑ ∫ = Jfff , (4.20) где iS означава билинеарну површ i , iN је број МоМ функција базиса дефинисаних на тој површи, а решавањем система добијамо тежинске коефицијенете ( )nka тих функција базиса. Коефицијенти ( )nka уз дублете тако ће бити израчунати два пута – по једном за сваку од две билинеарне површи на којима су дефинисани, па се узима њихова средња вредност. Интеграли у изразу (4.20) рачунају се нумерички, коришћењем Гаус- Лежандрове интегралне формуле [82]. Том приликом се вредност ових интеграла апроксимира двоструком сумом по p и s координатама билинеарне површи. У суми учествују вредности подинтегралних функција, у тачкама са тачно одређеним p и s координатама, помножене за ову врсту интеграције познатим тежинским коефицијентима. Вредности функција базиса kf које фигуришу у изразу (4.20) су познате у произвољној тачки, али вредности корекционе струје ( )nsJ∆ морају се израчунати у тачкама са координатама које одговарају Гаус-Лежандровој интеграцији. Те координате зависе од укупног броја тачака по p и s оси. Зато најпре треба утврдити колико је чланова у интегралној суми потребно за задовољавајућу тачност интеграције, па затим на основу тога одредити и тачне вредности p и s координата у којима ћемо рачунати корекционе струје ( )nsJ∆ . Сваки од iN коефицијената ( )n ka у систему (4.20) описује или p или s компоненту струје, док струја ( )nsJ∆ у једној тачки садржи и p и s компоненту. Зато, да би једначине система (4.20) биле независне, потребно је да укупан број тачака у којима се рачуна ( )nsJ∆ буде једнак или већи од половине броја МоМ функција базиса на посматраној билинеарној површи. На основу израза за 65 апроксимацију s компоненте струје (4.9) закључујемо да на једној билинеарној површи имамо ( )1sp −nn синглета и по pn дублета на свакој од ивица са координатом 1±=s . Међутим, како су дублети (ако постоје) заједнички за две билинеарне површи, њихов број по једној билинеарној површи је 2pn , па је број МоМ функција базиса једне билинеарне површи за s компоненту струје ( ) psp 1 nnn +− . Аналогно, број МоМ функција базиса једне билинеарне површи за p компоненту струје је ( ) sps 1 nnn +− . Сабирањем броја потребних функција базиса за p и s компоненту струје добијамо укупан број МоМ функција базиса по једној билинеарној површи, spMoM 2 nnN = . (4.21) Нумерички екесперименти показују да се стабилан резултат интеграције добија узимајући 1p +n тачака по p координати и 1s +n тачака по s координати, чиме је укупан број тачака у којима се рачунају корекционе струје ( ) ( )11 spJ ++= nnN . (4.22) Резимирајмо налажење корекционих коефицијената у PO зони. Најпре рачунамо расејано магнетско поље према изразу (3.4) у дискретним тачкама површи одређеним након анализе модела и усвајања редова апроксимације струја. Затим у истим тачкама рачунамо корекционе струје према изразу (3.5). Коначно добијамо корекционе коефицијенте ( )nka решавањем система једначина (4.20) за сваку билинеарну површ структуре. За коефицијенте који се појављују у две билинеарне површи узимамо средњу вредност из две израчунате. 4.3 Креирање макро функција базиса Корекционе струје у антенским проблемима делимо на МоМ и PO део, према изразу (3.24). Код расејачких проблема не постоји МоМ део, а корекционе струје делимо на исти начин као у PO делу антенског проблема. 66 Једина макро функција базиса у МоМ зони антенског проблема добија се као линеарна комбинација свих МоМ функција базиса из МоМ зоне, према изразу (3.37), где су тежински коефицијенти корекциони коефицијенти из МоМ зоне, добијени решавањем система једначина (3.36). Изузетак је прва итерација где се као макро функција базиса у МоМ зони користи почетно решење делимичног проблема, сходно изразу (3.38). Макро функције базиса у PO зони настају „поделом“ корекционих струја из PO зоне на M делова, као у изразу (3.26). Изузетак је нулта итерација у анализи расејача, где се макро функције базиса добијају поделом PO струја, према изразу (3.2). Овде се постављају питања избора вредности за M и начина на који поделу на макро функције базиса треба извршити. Према изразу (3.27), који важи за антене, али и расејаче под условом 0MoM =N , свака макро функција базиса у PO зони је линеарна комбинација свих МоМ функција базиса из PO зоне. Ово је математички погодан опис којим на јединствен начин представљамо било коју макро функцију базиса. Израз (3.27) може се приказати у облику ( ) Mlab N Nk k n k n kl n l ,...,1, 1 )()( MoM == ∑ += fF , (4.23) где је члан ( )nklb релативна тежина са којом корекциони коефицијент ( )nka МоМ функције базиса kf улази у l -ту макро функцију базиса (у n -тој итерацији). Поређењем израза (4.23) са изразом (3.27) очито је )()()( n k n kl n kl aba = , (4.24) па заменом овог израза за ( )nkla у (3.28) добијамо 1 1 )( =∑ = M l n klb . (4.25) Бавићемо се само случајевима када се тежине не мењају у различитим итерацијама, тј. ( ) klnkl bb = , одакле се изрази (4.23) и (4.25) своде на 67 ( ) Mlab N Nk k n kkl n l ,...,1, 1 )( MoM == ∑ += fF , (4.26) 1 1 =∑ = M l klb . (4.27) Свакој макро функцији базиса придружићемо групу МоМ функција базиса (над којма је дефинисана). За МоМ функцију базиса kf кажемо да припада l -тој групи ако је тежина са којом та МоМ функција базиса улази у l -ту макро функцију базиса различита од нуле, 0≠klb . Поступак одређивања група зваћемо груписањем. У свакој итерацији над истом групом креира се нова макро функција базиса (са другачијим вредностима корекционих коефицијената ( )nka ). Та макро функција базиса у апроксимацију коначног решења улази помножена једним тежинским коефицијентом. Тако је тежински коефицијент за све МоМ функције базиса из једне групе исти и не утиче на релативан однос између корекционих коефицијената тих МоМ функција базиса. Што је овај однос боље „погођен“, то ће решење бити боље. А то даље упућује на закључак да МоМ функције базиса из исте макро функције базиса треба да припадају физички блиским и глатким површима (с обзиром на „РО порекло“ корекционих струја). С обзиром на такву просторну условљеност група, МоМ функције базиса дефинисане на истој билинеарној површи – плочици – треба да припадају истој групи, као и МоМ функције базиса околних плочица. Стога је једноставније и визуелно подесније груписати плочице структуре, па на основу група плочица одредити групе МоМ функција базиса. Групу плочица по потреби ћемо звати макро површ. У складу са оваквом поставком, потребно је и плочице структуре поделити на МоМ зону и PO зону. У PO зони могу бити само металне плочице које су део затворене површи. Жице и све остале плочице морају бити у МоМ зони. Код антенских проблема све функције базиса које припадају жицама и плочицама из МоМ зоне смештају се у једну групу. За груписање плочица из PO зоне приказаћемо два аутоматска поступка. Један се одликује једноставношћу поступка и брзим радом, од другог очекујемо већу тачност при истом броју група. 68 Збога начина рада први ћемо звати запреминским, а други површинским груписањем. 4.3.1 Издвајање групе плочица за МоМ зону Корекционе површинске струје, дате изразом (3.5), имају смисла само на металним површима који су део затворене структуре. Плочице које чине такве површи могу се налазити у PO зони. Жице и све остале плочице морају припадати МоМ зони. Пошто одређивање корекционих струја у МоМ зони, путем решавања система једначина (3.36), може да се примени на произвољни део структуре, МоМ зони могу припадати и плочице које би могле да буду у PO зони. Могуће је, дакле, поред плочица за које је то обавезно, МоМ зони придружити и део околних плочица за које то није обавезно. Пошто ће од плочица из МоМ зоне (и жица) бити начињен посебан (делимични) модел ради издвојене анализе, потребно је да плочице које су повезане у читавом моделу остану повезане и у делимичном (да би се добио смислен резултат за делимични модел), тј. група плочица МоМ зоне у том случају треба да је повезана. 4.3.2 Запреминско груписање плочица из PO зоне Овде паралелепипед у који се може сместити („уписати“) структура делимо на одговарајући број мањих (елементарних), запремински равномерно распоређених паралелепипеда, да би за сваки елементарни паралелепипед формирали групу плочица чији се центри налазе у паралелепипеду. Затим се тако добијене групе деле на повезане делове без оштрих углова. У наставку ће нам бити потребни појмови суседства плочица и орта плочице. За две плочице кажемо да су суседне ако имају бар једну заједничку ивицу. Ортом плочице зваћемо спољашњи јединични вектор нормале у њеном центру, где је вектор положаја „центра“ плочице ( )cccc ,, zyx=r дат изразима (4.2) и (4.3а). Алгоритам запреминског груписања може се исказати наредним корацима. 1. Усвајамо Декартов координатни систем и координате свих чворова (темена плочица) структуре представљамо у том координатном систему. 69 2. За сваку плочицу одређујемо њен орт. 3. Међу координатама свих чворова налазимо минималне и максималне вредности по све три координатне осе: minx , maxx , miny , maxy , minz и maxz . 4. Усвајамо ( ) ( ) ( )2minmax2minmax2minmax zzyyxxD −+−+−= за максималну димензију структуре. 5. Задајемо проценат максималне димензије структуре, a ( 1000 << a ). 6. Одређујемо почетну дужину корака груписања, 100Dad = . 7. Одређујемо број корака по координатним осама, ( ) 1/minmaxx +−= dxxn , ( ) 1/minmaxy +−= dyyn , ( )  1/minmaxz +−= dzzn , где угласта заграда означава функцију „цео део“. 8. Одређујемо коначну дужину корака по координатним осама, ( ) xminmaxx / nxxd −= , ( ) yminmaxy / nyyd −= , ( ) zminmaxz / nzzd −= . 9. Формирамо zyx nnnM = група плочица, које су иницијално празне. 10. Сваку плочицу (која не припада МоМ зони) придружујемо групи zzyzyx ininnii ++= , где су xcx dxi = , ycy dyi = , zcz dzi = . 11. Бришемо све празне групе. 12. Задајемо критични угао, δ ( pi<δ<0 ). 13. На основу „посматране“ групе, у привремену групу иницијално издвајамо прву плочицу посматране групе, затим у итеративном поступку привременој групи прудружујемо све плочице посматране групе које су суседне са неким плочицама привремене групе при чему орт те плочице ни са једним ортом суседних плочица не заклапа угао већи од задатог критичног угла. Ако након овог итеративног поступка у посматраној групи постоје преостале плочице, креирамо нову групу и смештамо их у њу, а плочице привремене групе смештамо у 70 посматрану групу. Поступак издвајања понављамо док свака од група није била на „посматрању“. 4.3.3 Површинско груписање плочица из PO зоне Овде најпре издвајамо ивице структуре и формирамо ивичне групе са сваке стране ивице. Затим преостале плочице групишемо у повезане групе. Важе појмови уведени код запреминског груписања. Алгоритам површинског груписања може се исказати наредним корацима. 1. Усвајамо Декартов координатни систем и координате свих чворова (темена плочица) структуре представљамо у том координатном систему. 2. За сваку плочицу одређујемо њен орт. 3. Међу координатама свих чворова налазимо минималне и максималне вредности по све три координатне осе: minx , maxx , miny , maxy , minz и maxz . 4. Усвајамо ( ) ( ) ( )2minmax2minmax2minmax zzyyxxD −+−+−= за максималну димензију структуре. 5. Задајемо проценат максималне димензије структуре, a ( 1000 << a ). 6. Задајемо „тродимензионални“ критични угао, D3α ( pi<α< D30 ). 7. Заједничку ивицу сваког пара суседних плочица чији ортови заклапају угао већи од D3α проглашавамо ивичним сегментом. 8. Задајемо „дводимензионални“ критични угао, D2α ( pi<α< D20 ). 9. Формирамо низ међусобно спојених ивичних сегмената – ивицу, проглашавајући ивицом први неискоришћени ивични сегмент, и додајући тој ивици са обе стране нове суседне ивичне сегменте, све док има суседних сегмената, или док не наиђемо на више од једног суседног сегмента који настављају ивицу, или док не наиђемо на сегмент који са полазним ивичним сегментом заклапа угао већи од 71 D2α (за потребе налажења овог угла сегменте посматрамо као векторе усмерене према заједничком референтном смеру кроз ивицу). 10. За новоформирану ивицу рачунамо њено тежиште (средњу вредност збира вектора положаја центара њених сегмената), налазимо средишњи сегмент ивице - онај чији је центар најближи тежишту ивице, из ивице избацујемо све сегменте чији је центар на растојању већем од 100Da од центра средишњег сегмента, и за тако добијену ивицу формирамо две ивичне групе плочица; једну ивичну групу чине плочице које се налазе са једне стране ивице (а садрже њене ивичне сегменте), другу ивичну групу чине плочице које се налазе са друге стране ивице; сегменте ивице након формирања две ивичне групе сматрамо искоришћеним. 11. Ако има још неискоришћених ивичних сегмената, враћамо се на корак 9, у противном прелазимо на корак 12. 12. Задајемо глобални критични угао, D3β . 13. Међу негруписаним плочицама налазимо „најравнију“ - ону која има највећу вредност угла γpi -n , где је γ сума углова између орта посматране плочице и ортова са њом повезаних негруписаних плочица чији центри нису на растојању већем од 100Da од центра посматране плочице, а n је број таквих плочица – и око ње формирамо групу са тих n плочица. 14. Ако има још негруписаних плочица, враћамо се на корак 13, у противном је груписање завршено. 4.3.4 Креирање група МоМ функција базиса За МоМ зону креирамо само једну групу МоМ функција базиса. У њу улазе најпре све МоМ функције базиса којима се апроксимира струја дуж жица. Затим додајемо МоМ функције базиса којима је апроксимирана струја на плочицама које припадају МоМ зони. Приметимо да издвајањем плочица МоМ зоне из читавог модела настају слободне ивице дуж границе са PO зоном. Стога дублети дуж ових 72 ивица, дефинисани у читавом моделу, у делимичном моделу не постоје (струја коју представљају у читавом моделу у делимичном је једнака нули), па улазе у групе PO зоне. Треба рећи и да постоји могућност да се плочице из МоМ зоне поделе на више група, или да се чак свака МоМ функција базиса из МоМ зоне прогласи једном групом. Нумерички екесперименти, међутим, показују да овакав поступак у случају антене са диелектричним деловима може довести до нумеричких нестабилности и дивергирања итеративног поступка, те смо се у овој дисертацији задржали на само једној групи МоМ функција базиса за читаву МоМ зону. За сваку групу плочица у PO зони креирамо по једну групу МоМ функција базиса. Сви синглети дефинисани над плочицама једне групе улазе у одговарајућу групу МоМ функција базиса. Дублет дефинисан над двема плочицама из исте групе улази у одговарајућу групу МоМ функција базиса. Дублет дефинисан над двема плочицама из различитих група PO зоне улази у сваку од две групе МоМ функција базиса. Дублет дефинисан над двема плочицама од којих је једна у МоМ зони улази само у групу плочице из PO зоне. 4.3.5 Конструкција Макро функција базиса коришћењем група МоМ функција базиса По одређивању група МоМ функција базиса, нове макро функције базиса одређујемо у свакој итерацији коришћењем корекционих коефицијената ( )nka . У МоМ зони уводимо по једну нову макро функцију базиса у свакој итерацији. Сматрајући да су МоМ функције базиса које припадају нултој групи (тј. МоМ зони) нумерисане редним бројевима од 1 до MoMN , ту макро функцију базиса одређујемо према изразу (3.37). Макро функције базиса у PO зони добијамо према изразу (4.26), при чему је тежина ( klb ) МоМ функције базиса kf која припада само једној групи 1, а тежина МоМ функције базиса kf која припада двема различитим групама 0,5. На тај начин задовољена је релација (4.27). 73 4.4 Рачунање матрице система Чланови матрице система дати су изразом у загради са леве стране једнакости (3.23) за расејачке проблеме, односно изразом у загради са леве стране једнакости (3.33) за антенске проблеме. Видимо да се у оба случаја рачунају на исти начин, једино се број тих чланова разликује. Наиме, код расејача се користе и макро функције базиса из нулте итерације, док се код антена једина макро функција базиса из нулте итерације користи касније, у првој итерацији. Поред тога, у свакој итерацији после нулте у антенским проблемима уводи се једна макро функција базиса више – она из МоМ зоне, које нема код расејачких проблема. Како очекујемо да добијемо решење у свега неколико итерација, тај вишак макро функција базиса у антенским проблемима након нулте итерације је занемарљив. С друге стране, код расејачких проблема макро функције базиса у нултој итерацији добијамо на основу PO решења, које даје струје само у „осветљеном“ делу структуре, док у остатку структуре усвајамо да су струје једнаке нули. То значи да је отприлике половина тако добијених макро функција базиса једнака нули, те можемо сматрати да је број макро функција базиса у нултој итерацији упола мањи него у наредним итерацијама, односно износи 2M . Матрица система (3.23) има, дакле, за око 2M виши ред од матрице система (3.33). Пошто у расејачким проблемима нема МоМ зоне, она је једноставнија за разматрање. При томе, присуство МоМ зоне не мења суштински особине матрице. Стога ћемо посматрати само матрицу система (3.23). Коефицијенти матрице система (3.23) рачунају се према изразу ( ) ( )∑ = N j k jm i jl ZZ 1 * (4.28) и садрже четири различита индекса: индекси i и k односе се на редни број итерације, а индекси l и m на редни број групе МоМ функција базиса. У матричном приказу редни број колоне формира се на основу индекса i и l , а редни број врсте на основу индекса k и m . Помоћна променљива ( )ijlZ рачуна се према изразу (3.20). Видимо да у том изразу фигуришу корекциони коефицијенти 74 l -те групе у i -тој итерацији, ( )ikla , Nk ,..,1= , као и коефицијенти j -те врсте МоМ матрице, jkz , Nk ,..,1= . При решавању проблема коришћењем „стандардног“ МоМ-а читава МоМ матрица мора бити на располагању, што у случају када је матрица високог реда представља значајно меморијско заузеће. За ефикасно израчунавање коефицијената PDM матрице, према изразу (4.28), довољно је да истовремено имамо на располагању само j -ту врсте МоМ матрице, jkz , Nk ,..,1= , јер помоћу ње можемо да израчунамо ( )ijlZ за све вредности параметара i и l , а на основу тога и j -ти члан у суми (4.28), за све вредности i , l , k и m . Тако PDM матрицу можемо израчунати у N циклуса, Nj ,..,1= , додајући одговарајућем коефицијенту PDM матрице допринос j -те врсте МоМ матрице у j -том циклусу. PDM матрица је далеко нижег реда од МоМ матрице, па се најчешће може лако сместити у оперативну меморију рачунара. Долажење до PDM матрице, за расејачке проблеме, у једној итерацији укључује неколико временски захтевних рачунарских процедура: 1. Израчунавање магнетског поља, у тачкама површи PO зоне, према изразу (3.4), 2. Израчунавање чланова МоМ матрице (2.34) према изразу (2.32), 2. Израчунавање помоћних чланова ( )ijlZ према изразу (3.20), и 3. Израчунавање чланова PDM матрице (3.23) према изразу (4.28). Број тачака у којима се рачуна магнетско поље мења се од случаја до случаја. Ако претпоставимо да имамо PN плочица, и да је за сваку од њих ред апроксимације по p координати pn и по s координати sn , тада је укупан број МоМ функција базиса, N , ( )spP 2 nnNN = . (4.29) 75 Довољан број тачака у којима рачунамо магнетско поље, TN , је на основу изреченог у претходном поглављу ( )( )11 spPT ++= nnNN . (4.30) На основу претходна два израза је ( )( ) sp sp T 2 11 nn nn NN ++ = . (4.31) Број тачака TN је, дакле, пропорционалан са N . Пошто је при рачунању магнетског поља у свакој од ових тачака потребно одредити допринос сваке од N МоМ функција базиса, можемо писати да је укупно време потребно за рачунање магнетског поља 2 H NAT = , (4.32) где је A коефицијент пропорционалности, који зависи од начина рачунања и карактеристика рачунара. МоМ матрица је реда N , па је укупно време потребно за њено израчунавање 2 M NBT = , (4.33) где је B коефицијент пропорционалности, који, као и A , зависи од начина рачунања и карактеристика рачунара. При израчунавању помоћних чланова ( )ijlZ индекси узимају вредности ni ,..,0= , Nj ,..,1= , Ml ,..,1= , па је укупан број оваквих чланова ( ) NMn 1+ , где је M број група МоМ функција базиса, а n број PDM итерација. Чланови ( )ijlZ рачунају се формално према изразу (3.20), а практично број ненултих чланова ( )ikla је у просеку MN , па је укупно време потребно за рачунање чланова ( )ijlZ ( ) ( ) 2Z 1 NnCnT += , (4.34) 76 где је C збирно време потребно за једно множење и једно сабирање комплексних бројева. Имајући у виду да у нултој итерацији имамо отприлике 2M макро функција базиса, PDM матрица система (3.23) има приближно ( )[ ] 25.0 Mn + чланова. Пошто је за рачунање сваког од чланова потребно N множења и N сабирања комплексних бројева, укупно време потребно за рачунање чланова PDM матрице је ( )[ ] NMnCT 2P 5.0+= . (4.35) Број нових чланова које PDM матрица добија у свакој итерацији је ( )[ ] ( )[ ] 222N 25.05.0 nMMnMnN =−−+= . (4.36) Посматрањем израза (3.20) и (4.28) закључујемо да чланови из претходне итерације остају непромењени, па је могуће ускладишти их након што су први пут израчунати, тј. не рачунати их у каснијим итерацијама. У том случају време потребно за рачунање само нових чланова PDM матрице је ( ) 2N 2nMCnT = . (4.37) Време потребно за решавање система једначина (3.23), тј. инверзију PDM матрице, директном методом је ( ) ( )[ ]3R 5,03 1 MnCnT += . (4.38) Тако је укупно рачунарско време потребно за извршење n PDM итерација ( ) ( ) ( )∑∑∑ === ++++= n i n i n i iTiTiTnTnTT 1 R 1 N 1 ZMHPDM . (4.39) Сматрајући да је nn ≈+1 у изразу (4.34), односно да је nn ≈+ 5,0 у изразу (4.38), израз (4.39) постаје 77 ( ) nNBACMCNnCMCNnCMT       ++++        ++= 22 2 22 2 3 3 PDM 222 . (4.40) Како је у пракси 10,,, <≈>> nNMCBA , (4.41) PDMT коначно можемо да представимо као ( ) nNBAnCNnCMT 22233PDM 22 +++= . (4.42) С обзиром на услове (4.41), PDMT углавном линеарно расте са бројем итерација. Квадратна зависност постаје значајна тек при врло великом броју итерација, а кубна зависност постаје битна не само при врло великом броју итерација, већ и при довољно великом броју група МоМ функција базиса, M . С обзиром да се, иначе велика, МоМ матрица изнова рачуна у свакој итерацији, и да је у једном тренутку довољно да у опеартивној меморији буде само једна њена врста, једина меморијски захтевна структура је PDM матрица. Пошто је њен ред сигурно мањи од ( )Mn 1+ , меморија потребна за њено складиштење у n -тој PDM итерацији (под условом да се реални број представља са 4 бајта) је [ ] ( ) 2 2 1000 18MB       +< M nMEM . (4.43) Изрази за потребно рачунарско време (4.42) и меморијско заузеће (4.43) изведени су за случај PDM анализе расејача, али важе практично у неизмењеном облику и у случају када се анализирају антене на великим платформама. 78 5. Примери примене методе момената вођене физичком оптиком на расејаче У овом поглављу методу момената вођену физичком оптиком (скраћено PDM) употребићемо за анализу расејача, а у наредном за анализу антена на великим платформама. Резултате ове анализе упоредићемо са резултатима одговарајуће анализом применом методе момената (МоМ-а), као и са теоријским предвиђањима. Као меру за поређење са МоМ-ом користићемо нормализовани резидијум ∑ = = N j j n n v N RR 1 2 )( )( norm 1 , (5.1) добијен на основу резидијума )(nR датог изразом (3.21) код расејача, односно изразом (3.32) код антенских проблема. У изразу (5.1) је N број МоМ функција базиса, а jv је дато изразом (2.33) и представља колону слободних чланова у МоМ систему једначина (2.34). Ако у изразима за )(nR , (3.21) или (3.32), усвојимо 0)( =nilc , нормализовани резидијум дат изразом (5.1) биће једнак 1. За решење идентично МоМ решењу овај резидијум треба да је 0, а за било које смислено решење које одступа од МоМ решења очекујемо вредност између 0 и 1. У даљем тексту под резидијумом подразумевамо нормализовани резидијум дат изразом (5.1). Такође, као меру за поређење користићемо и карактеристичне резултате електромагнетске анализе: радарски одзив – RCS (Radar Cross Section) – код расејача, односно појачање - Gain – код антена на платформама. Математички опис структуре биће дат у Декатровом координатном систему, док ћемо RCS и Gain приказати у за њега везаном сферном координатном систему са слике 5.1 ( pi≤φ≤ 20 , 22 pi≤θ≤pi− ). 79 Слика 5.1. Координатни систем за математички опис структуре и приказ резултата анализе. Већина анализа извршена је на персоналном рачунару са четворојезгарним процесором такта GHz66,2 , RAMGB4 меморије такта MHz1066 , графичком картицом са 560GTX графичким процесором и диском од TB1 . 5.1 Савршено проводна сфера у побудном електромагнетском пољу Савршено проводна сфера, моделована са 24.936 плочица (као на слици 5.2), налази се у вакууму у пољу равног кружно поларизованог електромагнетског таласа који наилази из правца одређеног угловима 4pi=φ и 4pi−=θ у координатном систему са слике 5.1. Полупречник сфере је 13,3 таласних дужина. За МоМ решење коришћен је други ред апроксимације за сваку плочицу, што резултује са укупно 103.872 МоМ функција базиса. За PDM решење користили смо три различита запреминска груписања, са 352, 706 и 1442 групе плочица. Ово су и бројеви група МоМ функција базиса. Површинско груписање у случају сфере даје врло сличне групе плочица, те нема смисла радити поређење два начина груписања. 80 Слика 5.2. Савршено проводна сфера, моделована са 24.936 билинеарних површи, налази се у пољу равног кружно поларизованог електромагнетског таласа који наилази из правца одређеног угловима 4pi=φ и 4pi−=θ у координатном систему са слике 5.1. Полупречник сфере је 13,3 таласних дужина. Груписања су приказана на слици 5.3. Групи припадају све повезане плочице исте боје. а) б) в) Слика 5.3. Три различита PDM модела сфере са слике 5.2 са а) 352, б) 706 и в) 1442 групе плочица. Урађена је PDM анализа за сва три модела са слике 5.3. На слици 5.4 приказана је зависност броја макро функција базиса од броја итерација, за сва три модела. (Слово V које стоји уз број група плочица у 81 легендама слика означава да је реч о „запреминском“ груписању плочица.) У свакој итерацији додајемо онолико макро функција базиса колико група има модел. Изузетно, из нулте итерације узимају се само макро функције базиса које су претежно „осветљене“ побудним таласом. Тај број је отприлике једнак половини укупног броја група, што се може и видети са слике 5.4. Слика 5.4. Зависност броја макро функција базиса од броја итерација за моделе са слике 5.3. На слици 5.5 приказана је зависност резидијума од броја итерација за сва три модела. Итерација 0 представља PO решење, које је исто за све моделе, па је и резидијум исти. Очекивано, резидијум се смањује у свакој итерацији, али је и поправка у свакој следећој итерацији мања. Примећујемо да са већим бројем група можемо постићи мањи резидијум у истом броју итерација (тј. бржу конвергенцију). Ово је такође разумљиво, с обзиром да више група значи и већи број макро функција базиса. Истовремено, свака макро функције базиса покрива мању област на којој треба да апроксимира релативан однос струја, па је и та апроксимација боља. 82 Оно што нас при анализи расејача пре свега интересује јесте његов радарски одзив, RCS. Тродимензионални приказ RCS -а није погодан за визуелно поређење резултата, па ћемо посматрати дводимензионални пресек 4pi=φ , који пролази кроз максимум RCS -а. На слици 5.6 приказано је поређење RCS -а добијеног на основу МоМ анализе и RCS -а добијеног на основу PO решења. Видимо да је слагање у главном лобу до dB25 испод максимума готово бескпрекорно, али да испод тих нивоа долази до значајног мимоилажења. Прва итерација код сва три PDM модела приказана на слици 5.3 даје врло сличне RCS резултате, а на слици 5.7 приказан је онај са 706 група плочица, и упоређен са МоМ резултатом. Очито је веома сличан МоМ резултату и представља значајно унапређење у односу на PO решење. Почевши од друге итерације, резултати добијени PDM-ом практично се не могу разликовати од МоМ резулатата. Са слике 5.5 се види да се резидијум посматраног модела у другој итерацији практично спустио на -310 , па се оквирно може рећи да PDM резултат са резидијумом од -310 обезбеђује RCS који се одлично слаже са оним добијеним МоМ-ом. Слика 5.5. Зависност резидијума од броја итерација за моделе са слике 5.3. 83 Слика 5.6. RCS сфере са слике 5.2 добијен на основу МоМ решења и добијен на основу PO решења. Слика 5.7. RCS сфере са слике 5.2 добијен на основу МоМ решења и добијен на основу прве итерације PDM анализе модела са слике 5.3б. 84 Рачунарско време потребно за PDM анализу у функцији броја итерација приказано је на слици 5.8. За моделе са 352 и 706 група пораст времена је приближно линеаран, а тек са повећањем броја група на 1441 видније долази до изражаја нелинеарна зависност времена од броја итерација, према изразу (4.42). Време потребно за прву итерацију, која даје веома добар резултат, је око 45 минута, а до решења након друге итерације стижемо за око сат и по. Одговарајуће време потребно за МоМ решавање (коришћењем паралелизоване LU- декомпозиције) дуже је од пет и по сати. Слика 5.8. Рачунарско време потребно за извршавање PDM итерација за моделе са слике 5.3. Рачунарско време потребно за PDM анализу у функцији резидијума приказано је на слици 5.9. Време за сваки од модела измерено је током једне анализе, уз уобичајен начин рада рачунара и укључен антивирусни софтвер. Са слике се види да су модели са већим бројем група временски ефикаснији, односно обезбеђују исти резидијум за краће време. 85 У вези са мерењем времена треба истаћи једну напомену. Наиме, два временски најзахтевнија процеса PDM анализе, рачунање магнетског поља и решавање PDM система линеарних једначина, изведена су, сваки посебно, коришћењем паралелног извршавања на језгрима графичког процесора (такође и „пуњење“ и „инверзија“ матрице код МоМ-а). Таква реализација обезбеђује знатно убрзање у односу на стандардне процесоре, али спорадично може довести и до нелинеарног понашања. Такође, на диску се складишти постојећи део PDM матрице (који се не мења у наредним итерацијама), као и помоћна структура Z (ради рачунања резидијума), па постоје случајности при упису и читању података које могу да доведе до мањих варијација у времену извршавања. У том смислу, измерена времена не могу се сматрати лабораторијским узорком, тј. постоје варијације у времену потребном за PDM анализу истог модела у више узастопних покретања. Међутим, разлике су (у овом случају) реда неколико процената, а главни циљ мерења времена јесте оквирна провера апроксимативног израза (4.42), па су у том смислу приказани резултати задовољавајуће тачности. Слика 5.9. Рачунарско време потребно за добијање одређеног резидијума за моделе са слике 5.3. 86 Број макро функција базиса употребљених за PDM анализу у функцији резидијума приказан је на слици 5.10. Види се да модели са мањим бројем група углавном обезбеђују идентичан резидијум са мањим бројем макро функција базиса. Види се, међутим, и тенденција да при довољно ниским резидијумима модели са већим бројем група преузму преимућство у ефикасности (исти резидијум са мањим бројем макро функција базиса). Слика 5.10. Број макро функција базиса потребан за добијање одређеног резидијума за моделе са слике 5.3. 5.2 Савршено проводна коцка у побудном електромагнетском пољу Облик коцке уноси ефекат оштрих ивица, који код сфере нисмо имали. У [80] приказан је пример коцке у вакууму у пољу равног линијски поларизованог електромагнетског таласа који наилази из правца одређеног угловима 4pi=φ и 4pi=θ у координатном систему са слике 5.1. Дужина ивице коцке је 4 таласне дужине. Свака странице коцке моделована је са 3636× једнаких плочица. За МоМ 87 решење коришћен је први ред апроксимације за сваку плочицу, што резултује са укупно 15.552 МоМ функције базиса. На слици 5.11 приказана су три различита начина креирања 216 група плочица. Најпре је направљено 216 група са по 36 наусмично извучених плочица. Алгоритам за бојење покушава да за групу одабере једну од 15 боја које немају њој суседне групе. Ако су све боје заузете, групи се додељује сива боја. Тако је настао приказ „насумичног“ груписања на слици 5.11а. На слици 5.11б приказано је 216 „униформних“ група, свака са 66× плочица. „Оивичено“ груписање сачињено је од група са 22× плочице у угловима, група са 28× плочице дуж ивица, и група са по 88× плочица у унутрашњости страница, као што је приказано на слици 5.11в. а) б) в) Слика 5.11. Различити начини креирања 216 група плочица коцке чија је свака страница моделована са 3636× једнаких плочица: а) насумично, б) униформно, в) оивичено. Резидијуми након примене PDM анализе на сваки од три модела приказани су на слици 5.12. Број група у сва три модела је исти, па разлика у резидијумима потиче од начина груписања. Очекивано, насумичан избор плочица резултује најслабијом конвергенцијом. Наиме, просторно блиске и повезане плочице имају сличну расподелу струја, па им одговара и сличан тежински коефицијент, што оваквим начином груписања није узето у обзир. С друге стране, код просторно повезаних група, оивичавање показује преимућство у односу на униформни модел, посебно при вишим итерацијама. На основу овог и многих других нумеричких експеримената закључено је да оивичавање увек доноси одређену предност у односу на сличан модел без оивичавања, те је оно примењено у површинском аутоматском алгоритму за груписање плочица. 88 Слика 5.12. Резидијуми PDM анализе за моделе са слике 5.11. Може се поставити и питање смисла одређивања корекционих струја у PO стилу. Да бисмо то проверили, корекционе коефицијенте за модел са униформним груписањем, који би иначе били добијени на основу PO корекционих струја, усвојили смо на случајан начин: модул сваког коефицијента као случајан реалан број између 0 и 1, а његову фазу као случајан реалан број између 0 и pi2 . Резидијум приказан на слици 5.12 („Насумичне корекционе струје“) јасно указује на веома слабу конвегенцију и потврђује смисао одређивања корекционих струја у PO стилу. Посматрајмо сада коцку која је у смислу електричне величине слична претходно анализираној сфери. Коцка је савршено проводна, моделована са 25.872 плочица (као на слици 5.13), и налази се у вакууму у пољу равног кружно поларизованог електромагнетског таласа који наилази из правца одређеног угловима 4pi=φ и 4pi−=θ у координатном систему са слике 5.1. Дужина ивице коцке је 20 таласних дужина. За МоМ решење коришћен је други ред апроксимације за сваку плочицу, што резултује са укупно 106.368 МоМ функција 89 базиса. Може се видети да плочице местимично одступају од квадратног облика – ово је допуштено с намером да сегментација буде равноправна и за правилне и за неправилне облике. Слика 5.13. Савршено проводна сфера, моделована са 25.872 билинеарне површи, налази се у пољу равног кружно поларизованог електромагнетског таласа који наилази из правца одређеног угловима 4pi=φ и 4pi−=θ у координатном систему са слике 5.1. Дужина ивице коцке је 20 таласних дужина. Посматраћемо по три PDM модела за оба аутоматска начина груписања. На слици 5.14 приказани су модели са 384, 720 и 1427 група добијених запреминским поступком груписања. На слици 5.15 приказани су модели са 378, 720 и 1422 група добијених површинским поступком груписања. Видимо да се број група се са сваким наредним моделом отприлике удвостручује. 90 а) б) в) Слика 5.14. Три различита PDM модела коцке са слике 5.13 са а) 384, б) 720 и в) 1427 група плочица, добијених запреминским поступком груписања. а) б) в) Слика 5.15. Три различита PDM модела коцке са слике 5.13 са а) 378, б) 720 и в) 1422 групе плочица, добијених површинским поступком груписања. Промена резидијума по итерацијама за моделе са запреминским груписањем приказана је на слици 5.16, а за моделе са површинским груписањем на слици 5.17. Резидијуми за површинско груписање конвергирају за нијансу брже, али је та разлика у пракси занемарљива. Понашање и редови величина резидијума веома су слични онима за сферу, као и већина осталих тамо приказаних резултата, па нема потребе да их приказујемо и за коцку. RCS у равни 4pi=φ , добијен на основу МоМ решења, упоређен је са RCS -ом добијеним на основу PO решења на слици 5.18, а са RCS -ом добијеним на основу прве итерације PDM анализе модела са слике 5.15в на слици 5.19. Видимо да већ PO решење даје добру апроксимацију главних лобова, док прва PDM итерација побољшава резултат у области са нивоима знатно мањим од максималног. 91 Слика 5.16. Зависност резидијума од броја итерација за моделе са слике 5.14. Слика 5.17. Зависност резидијума од броја итерација за моделе са слике 5.15. 92 Слика 5.18. RCS коцке са слике 5.13 добијен на основу МоМ решења и на основу PO решења. Слика 5.19. RCS коцке са слике 5.13 добијен добијен на основу МоМ решења и на основу прве итерације PDM анализе модела са слике 5.16в. 93 Разлика између PO решења и решења добијеног у првој итерацији PDM-а, и каснијим, јаснија је када се погледа блиско поље. На слици 5.20 приказан је интензитет електричног поља у једном пресеку коцке, тако да појединим нивоима поље одговарају боје у легенди слика. У коцки нема поља, па би унутрашњи квадрат требало да буде потпуно плав, као на слици 5.20в, која приказује резултат пете PDM итерације. PO апроксимација струја не урачунава струје у зони сенке, па тај део контуре на слици 5.20а (горња и лева странице) није добро дефинисан, али се и поље ван коцке у околини те контуре разликује од оног на слици 5.20в. Већ прва PDM итерација, чији је резултат приказан на слици 5.20б, представља значајно побољшање. Резултат након пете PDM итерације, приказан на слици 5.20в, практично се не разликује од резултата добијеног на основу МоМ решења. а) б) в) Слика 5.20. Интензитет електричног поља у једном пресеку у околини коцке, добијен на основу а) PO решења, и решења добијеног PDM анализом модела са слике 5.15в у б) првој и в) петој итерацији. 5.3 Авион у побудном електромагнетском пољу Авион, моделован са 13.070 савршено проводних плочица (као на слици 5.21), налази се у вакууму у пољу равног кружно поларизованог електромагнетског таласа који наилази из правца одређеног угловима 4pi=φ и 4pi−=θ у координатном систему са слике 5.1. Дужина авиона је 60,8 таласних дужина. За МоМ решење коришћен је други ред апроксимације за сваку плочицу, што резултује са укупно 104.276 МоМ функција базиса. 94 За PDM решење користили смо по четири различита запреминска и површинска груписања, са 160, 318, 669 и 1333 групе плочица добијених запреминским поступком груписања, приказаних на слици 5.22, односно 157, 317, 660 и 1322 групе плочица добијених површинским поступком груписања, приказаних на слици 5.23. Слика 5.21. Савршено проводни авион, моделован са 13.070 билинеарних површи, налази се у пољу равног кружно поларизованог електромагнетског таласа који наилази из правца одређеног угловима 4pi=φ и 4pi−=θ у координатном систему са слике 5.1. Дужина авиона је 60,8 таласних дужина. Промена резидијума по итерацијама за моделе са запреминским груписањем приказана је на слици 5.24, а за моделе са површинским груписањем на слици 5.25. Резидијуми за површинско груписање конвергирају за нијансу брже, али је та разлика у пракси занемарљива. Као и код сфере и коцке, брзина конвергенције опада са порастом броја итерација. Очито, макро функције базиса које уводимо на почетку итеративног поступка „богатије“ су новим садржајем у односу на текуће решење (математичари би рекли „ортогоналније“) него оне у каснијим итерацијама. Да ли постоји могућност да итеративни поступак почне да дивергира? Сем у случају нумеричких нестабилности (нпр. коришћење полиномских функција базиса високог реда уз представљање реалног броја са 95 четири, уместо са осам бајтова), очекујемо сигурну конвергенцију. Начин формирања нових макро функција базиса статистички обезбеђује минимум ортогоналности у односу на постојеће макро функција базиса, довољан за сигурну (премда све спорију) конвергенцију. Шта се дешава за врло велики број итерација није од практичног интереса – из итеративног поступка изаћи ћемо много пре. а) б) в) г) Слика 5.22. Четири различита PDM модела авиона са слике 5.21 са а) 160, б) 318, в) 669 и г) 1333 група плочица, добијених запреминским поступком груписања. 96 Можемо приметити и да су резидијуми, добијени у истом броју итерација и за сличан број група, код авиона безмало за ред величине већи него при анализи сфере и коцке. Ово јасно наговештава да је метода боље прилагођена облицима без оштрих ивица и/или са великим равним површима, што је, с обзиром на „PO порекло“, очекивано. а) б) в) г) Слика 5.23. Четири различита PDM модела авиона са слике 5.21 са а) 157, б) 316, в) 660 и г) 1322 група плочица, добијених површинским поступком груписања. 97 Слика 5.24. Зависност резидијума од броја итерација за моделе са слике 5.22. Слика 5.25. Зависност резидијума од броја итерација за моделе са слике 5.23. 98 На слици 5.26 упоређен је RCS у равни 4pi=φ добијен на основу МоМ решења са: RCS -ом добијеним на основу PO решења (слика 5.26а), и RCS -ом добијеним на основу PDM анализе модела са слике 5.23г након: прве итерације (слика 5.26б), друге итерације (слика 5.26в) и пете итерације (слика 5.26г). За разлику од коцке, овде је прва PDM итерација сасвим јасно побољшање у односу на PO решење, видан је и помак између прве и друге итерације, док је резултат након пете PDM итерације практично идентичан резултату добијеном на основу МоМ решења. а) б) в) г) Слика 5.26. RCS авиона са слике 5.21 добијен на основу МоМ решења упоређен са RCS-ом добијеним на основу (а) PO решења, и на основу PDM анализе модела са слике 5.23г након (б) прве итерације, (в) друге итерације и (г) пете итерације. Током PDM итеративног поступка можемо пратити вредност резидијума, али, у реалној ситуацији, RCS који даје текуће PDM решење не можемо упоредити са МоМ решењем (које немамо). Да ли на основу резидијума можемо проценити колико ће нам RCS резултат одступати од оног добијеног МоМ-ом? С тим у вези, 99 погледајмо слику 5.27, на којој је приказана RCS грешка у функцији резидијума за модел са слике 5.23г. Ова је грешка израчуната као средња вредност апсолутног одступања RCS -а добијеног PDM анализом од RCS -а добијеног МоМ анализом, у 1441 тачки по θ оси (када ова средња вредност јасно конвергира). На преломним тачкама ове сегментне криве означене су итерације на које се грешка односи. (И остали модели са слике 5.23 дају врло сличне резултате.) Увидом у слику 5.26 можемо видети како се „визуелна“ разлика пресликава у само један број. Заправо, већ PO решење даје добре главне лобове, али се видно одступање у минимумима RCS-а пресликава у преко четири децибела грешке (PO решење је нулта итерација). Коме су само главни лобови важни, могао би да стане са итерацијама већ овде. RCS након пете PDM итерација се визуелно готово не разликује од оног добијеног из МоМ решења. Грешка – око пола децибела. Груб закључак – резидијум од -310 даје одличан резултат, већ -210 је довољно са практичног становишта. Слика 5.27. RCS грешка за модел са слике 5.23г добијена је као средња вредност апсолутног одступања од RCS-а добијеног МоМ-ом, у 1441 тачки. 100 Рачунарско време потребно за PDM анализу модела са слике 5.23, у функцији резидијума, приказано је на слици 5.28, а број макро функција базиса употребљених за PDM анализу ових модела (такође у функцији резидијума) приказан је на слици 5.29. Графици су слични одговарајућим графицима код сфере. Оно што овде вреди додати је да се избор оптималних параметара PDM анализе креће у „троуглу“ Резидијум-Време-Број макро функција базиса, где је могуће одабрати два параметра, али је тада онај трећи одређен двама одабранима. На пример, одговор на питање „како да за најкраће време добијем жељени резидијум?“ је, гледајући слику 5.28, „изабери што већи број група функција базиса“, а са слике 5.29 се тада види да то плаћамо (бар у приказаном опсегу вредности резидијума) већим бројем макро функција базиса (тј. већим меморијским заузећем). Слика 5.28. Рачунарско време потребно за добијање одређеног резидијума за моделе са слике 5.23. 101 Слика 5.29. Број макро функција базиса потребан за добијање одређеног резидијума за моделе са слике 5.23. 5.4 Хеликоптер у побудном електромагнетском пољу Хеликоптер, моделован са 12.821 савршено проводном плочицом (као на слици 5.30), налази се у вакууму у пољу равног кружно поларизованог електромагнетског таласа који наилази из правца одређеног угловима 43pi−=φ и 4pi−=θ у координатном систему са слике 5.1. Дужина хеликоптера је 52,8 таласних дужина. За МоМ решење коришћен је други ред апроксимације за сваку плочицу, што резултује са укупно 101.886 МоМ функција базиса. За PDM решење користили смо по три различита запреминска и површинска груписања, са 404, 674, и 1343 групе плочица добијених запреминским поступком груписања, приказаних на слици 5.31а-в, односно 384, 673 и 1343 групе плочица добијених површинским поступком груписања, приказаних на слици 5.31г-ђ. 102 Слика 5.30. Савршено проводни хеликоптер, моделован са 12.821 билинеарном површи, налази се у пољу равног кружно поларизованог електромагнетског таласа који наилази из правца одређеног угловима 4pi=φ и 4pi−=θ у координатном систему са слике 5.1. Дужина авиона коцке је 52,8 таласних дужина. Промена резидијума по итерацијама за моделе са запреминским груписањем приказана је на слици 5.32, а за моделе са површинским груписањем на слици 5.33. Као и у случају коцке и авиона, резидијуми за површинско груписање конвергирају за нијансу брже, али је та разлика у пракси занемарљива. Начелно се може закључити да оба поступка груписања обезбеђују сличне резултате, и да је сама метода стабилна по питању избора група, уколико су оне физички повезане и припадају глатким деловима површи. У наредном поглављу, при анализи антена на великим платформама, користићемо само запреминско груписање. То очито неће утицати на резултат, али знатно убрзава сам поступак груписања, будући да је запреминско груписање знатно једноставније, па се због тога и брже извршава на рачунару. 103 а) г) б) д) в) ђ) Слика 5.31. PDM модели хеликоптера са слике 5.30, са а) 404, б) 674 и в) 1343 група плочица добијених запреминским поступком груписања, г) 384, д) 673 и ђ) 1343 група плочица добијених површинским поступком груписања. На слици 5.34 је RCS у равни 43pi−=φ , добијен на основу МоМ решења, упоређен са: RCS -ом добијеним на основу PO решења (слика 5.34а), и са RCS -ом добијеним на основу PDM анализе модела са слике 5.31ђ након: прве итерације (слика 5.34б), друге итерације PDM (слика 5.34в) и пете итерације (слика 5.34г). Резултати су веома слични онима код авиона. 104 Слика 5.32. Зависност резидијума од броја итерација за моделе са слике 5.31а-в. Слика 5.33. Зависност резидијума од броја итерација за моделе са слике 5.31г-ђ. 105 а) б) в) г) Слика 5.34. RCS хеликоптера са слике 5.30 добијен на основу МоМ решења упоређен са RCS-ом добијеним на основу (а) PO решења, и на основу PDM анализе модела са слике 5.31ђ након (б) прве итерације, (в) друге итерације и (г) пете итерације. 106 6. Примери примене методе момената вођене физичком оптиком на антене на платформама У овом поглављу методу момената вођену физичком оптиком (скраћено PDM) употребићемо за анализу антена на великим платформама. Као платформа послужиће летелице анализиране у претходном поглављу. На авион ћемо монтирати полуталасни дипол, а на хеликоптер микрострип печ антену. Реч је о слабо усмереним антенама, а резултати су слични и ако ове антене размене платформе. „Квалитет“ решења пратићемо преко резидијума, датог изразом (5.1), и преко појачања антене на платформи. Добијено појачање упоредићемо са оним добијеним анализом применом методе момената. Изузетак је завршни пример, са безмало милион и по непознатих (тј. МоМ-ом добијених линеарних једначина), за који је решавање LU-декомпозицијом у тренутку писања овог текста готово неизводљиво на персоналним рачунарима. 6.1 Микрострип печ антена монтирана на хеликоптер Микрострип печ антена монтирана је на хеликоптер, као на слици 6.1. Нос хеликоптера је у Декартовом координатном систему усмерен дуж x− осе, а z оса нормална је на подножје и усмерена ка елиси на крову. Структура је сачињена од 12.022 плочице, а хеликоптер је на радној учестаности дуг 58,7 таласних дужина. Диелектрик печ антене моделован је као савршен, релативне пермитивности 4,3. Усвојен је други ред апроксимације за струју на свакој плочици што резултује са укупно 95.730 МоМ функција базиса. При изради PDM модела одлучујемо се за величину МоМ зоне и број група у PO зони. Најпре ћемо фиксирати величину МоМ зоне - у њу ћемо укључити плочице које макар делом припадају помоћној сфери полупречника четири таласне дужине, чији је центар у средишту печ антене, а која је приказана на слици 6.2а. На слици 6.2б су плочице из МоМ зоне обојене у сиво, а плочице из 107 PO зоне у жуто. Од плочица МоМ зоне и жица креирамо делимичан МоМ модел, приказан на слици 6.2в. Слика 6.1. Микрострип печ антена монтирана на подножје хеликоптера. Структура је сачињена од 12.022 плочице, а хеликоптер је на радној учестаности дуг 58,7 таласних дужина. Усвојен је други ред апроксимације за струју на свакој плочици што резултује са укупно 95.730 МоМ функција базиса. Потребно је још да групишемо плочице из PO зоне. Видели смо у претходном поглављу да између површинског и запреминског поступка груписања нема значајне разлике у резултатима PDM анализе, па ћемо овде користити запреминско груписање, које је једноставније. Да бисмо испитали утицај броја група, посматраћемо шест различитих модела, код којих је PO зона подељена на 342, 698, 1345, 2694, 5016 и 10064 групе. Модел са 698 група приказан је на слици 6.2г. Зависност резидијума од броја итерација приказана је на слици 6.3. Нулта итерација представља МоМ решење делимичног модела са слике 6.2в. За моделе са већим бројем група урађено је мање итерација, тако да укупан број Макро функција базиса не прелази значајно 10.000. Резидијуми за све моделе опадају у свакој итерацији, а брзина опадања је већа код модела са већим бројем група. Као и код расејача, брзина конвергенције опада са порастом броја итерација. Очито је и овде (као и код расејача) да „PO правац поправке“ најбоље резултате даје у почетним итерацијама. 108 а) б) в) г) Слика 6.2. (а) Помоћна сфера полупречника λ4 (четири таласне дужине) са центром у средишту печ антене. (б) Сиве плочице су макар делом у помоћној сфери и издвајамо их у МоМ зону, док жуте плочице додељујемо PO зони. (в) Од плочица из МоМ зоне и жица креирамо делимичан МоМ модел. (г) Плочице из PO зоне групишемо запреминским поступком; овде је приказан модел са 698 група. Пошто као критеријум за конвергенцију посматрамо резидијум, и овде се, као и код расејача, поставља питање да ли је резидијум верна слика других резултата од интереса – код антене на платформи тај „други“ резултат је, пре свега, појачање. Ово питање можемо разложити на два: прво, да ли смањивањем резидијума „поправљамо“ резултат за појачање (приближавамо га „МоМ решењу“), и друго, да ли се та поправка, у функцији резидијума, понаша слично за различите моделе. Ова су питања код антена на платформама „осетљивија“ него код расејача, с обзиром да су тежине функција базиса у зони напајања обично знатно веће у односу на остатак модела. 109 Слика 6.3. Зависност резидијума од броја итерација за PDM моделе са различитим бројем група у PO зони. На слици 6.4а приказано је појачање у симетралној равни хеликоптера ( 0=φ ) добијено на основу резултата нулте итерације PDM анализе (тј. МоМ анализе делимичног модела са слике 6.2в). Овај резултат исти је за сваки од разматраних PDM модела, а резидијум је око 0,15. На сликама 6.4б и 6.4г приказана су појачања добијена на основу резултата прве и треће итерације PDM анализе модела са 2964 групе. Чак се и голим оком види да је смањивање резидијума овде праћено приближавањем криве појачања оној добијеној након МоМ анализе. Исто важи и за све остале моделе, чији резултати овде нису приказани. Закључујемо да је смањивање резидијума праћено приближавањем појачања добијеног на основу PDM анализе ономе добијеном на основу МоМ анализе. Посматрајући слику 6.3 видимо да сличну вредност резидијума (око 0,015) дају модел са 5016 група у другој итерацији, модел са 2694 групе у трећој итерацији, модел са 1345 група у четвртој итерацији и модел са 698 група у осмој итерацији. Криве појачања за ове моделе у поменутим итерацијама приказане су 110 на сликама 6.4б-ђ. На основу овог, и бројних других примера који овде нису приказани, закључујемо да сличне вредности резидијума обезбеђују и сличне резултате за појачање. а) б) в) г) д) ђ) Слика 6.4. Појачање у симетралној равни хеликоптера добијено на основу резултата МоМ анализе упоређено са појачањима добијеним на основу резултата PDM анализе у различитим итерацијама и за различите PDM моделе. Зависност броја макро функција базиса од броја итерација приказана је, за три модела са мањим бројем група, на слици 6.5. Овде је извршена мала измена у 111 односу на PDM теорију изложену у трећем поглављу. Наиме, свака МоМ функција базиса на граници PO зоне и МоМ зоне у првој се итерацији (и само у првој) уводи као нова макро функција базиса. То значи да се ове МоМ функције базиса у првој и свакој следећој итерацији независно подешавају (имају сопствени тежински коефицијент, за разлику од груписаних МоМ функција базиса, где је тежински коефицијент заједнички за целу групу). На тај начин даје се већи степен слободе функцијама базиса које „премошћавају“ две зоне. Тако је број новоуведених макро функција базиса у првој итерацији GNM + , где је M број група у PO зони, а GN број „граничних“ МоМ функција базиса. Пошто је GN обично знатно већи од један, у првој се итерацији додаје знатно више макро функција базиса него у свакој наредној (у којима се додаје 1+M макро функција базиса, по једна за сваку од M група + једна у МоМ зони). Ово се види и на слици 6.5, где између нулте и прве итерације видимо већи нагиб него између осталих итерација. Слика 6.5. Зависност броја макро функција базиса од броја итерација за PDM моделе са различитим бројем група у PO зони. 112 Зависност потребног рачунарског времена од броја итерација приказана је на слици 6.6. За мерење времена важе напомене изречене у претходном поглављу, код анализе расејача: за сваки од модела приказано је време измерено током једне анализе, уз уобичајен начин рада рачунара и укључен антивирусни софтвер. Измерена времена у складу су са процењеним у поглављу четири. Код модела са 1345 група нелинеарност постаје уочљива при вишим итерацијама, а при даљем повећању броја група уочљива је и при мањем броју итерација. Слика 6.6. Укупно рачунарско време потребно да се изврши одређени број итерација за PDM моделе са различитим бројем група у PO зони. На слици 6.7 приказано је укупно рачунарско време потребно да поједини PDM модели достигну одређену вредност резидијума. Видимо да сваку вредност резидијума постоји „оптимална“ крива (тј. модел) која даје резидијум близак траженом за најкраће време. Тако, на пример, модел са 10.064 група даје резидијум од око 0,023 брже него модел са 698 група, али је за око 25% спорији од модела са 2964 групе, који при томе даје и бољи резидијум, од око 0,021. 113 Слика 6.7. Укупно рачунарско време потребно да се достигне одређени резидијум за PDM моделе са различитим бројем група у PO зони. Слика 6.8. Укупан број макро функција базиса потребних да се достигне одређени резидијум за PDM моделе са различитим бројем група у PO зони. 114 На слици 6.8 приказан је укупан број макро функција базиса потребних да поједини PDM модели достигну одређену вредност резидијума. За веће резидијуме ефикаснији су модели са мањим бројем група, а са смањивањем резидијума модели са већим бројем група постепено преузимају преимућство. Погледајмо сада утицај величине МоМ зоне у PDM моделу на брзину конвергенције. Полазећи од модела са слике 6.2г, направили смо два нова, који се од полазног разликују само по величини МоМ зоне – у једном моделу свели смо је на област сфере полупречника две таласне дужине, у другом на област елипсоида чији је полупречник по z оси остао на четири таласне дужине, а по преостале две осе повећан на осам таласних дужина. Зато ћемо два новонастала модела означити са λ2 и λ8 , а полазни са λ4 . Одговарајући делимични модели, начињени од жица и плочица из МоМ зоне, приказани су на сликама 6.9а-в. Параметре груписања нисмо мењали, па се број група симболично разликује - код λ2 модела незнатно се повећао (за групе плочица изузетих из МоМ зоне), а код λ8 модела незнатно смањио (за групе плочица придодатих МоМ зони). Резидијум за ова три модела, у функцији броја итерација, приказан је на слици 6.10. Видимо да су криве практично само померене, при чему се већом МоМ зоном постижу нешто ниже вредности резидијума. Сличне криве добијају се и за већи број група, што овде није приказано. Закључујемо зато да се повећањем МоМ зоне крива резидијума у функцији броја итерација практично само спушта наниже у извесној мери, зависно од процента увећања МоМ зоне. а) б) в) Слика 6.9. Делимични модели настали издвајањем плочица МоМ зоне из модела са слике 6.1, када је МоМ зона ограничена (а) сфером полупречника две таласне дужине – модел λ2 , (б) сфером полупречника четири таласне дужине – модел λ4 , и (в) елипсоидом чији је полупречник по висини хеликоптера четири таласне дужине, а по преостале две ортогоналне осе осам таласних дужина – модел λ8 . 115 Слика 6.10. Зависност резидијума од броја итерација за PDM моделе са различитом величином МоМ зоне и приближно истим бројем група у PO зони. На сликама 6.11а-б приказан је резидијум делимичних модела λ2 и λ8 , док је резидијум делимичног модела λ4 приказан на слици 6.4а. Видимо да се повећањем МоМ зоне појачање приближава оном које произилази из МоМ решења, али чак ни λ8 модел не обезбеђује добар резултат ван главног лоба криве, премда он, као што се да видети са слике 6.9в, обухвата значајан део плочица модела у околини печ антене. а) б) Слика 6.11. Појачање у симетралној равни хеликоптера добијено на основу почетног решења PDM анализе – тј. МоМ решења делимичног модела (а) са слике 6.9а, и (б) са слике 6.9в. 116 Стижемо и до „троугла“ резидијум-време-број макро функција базиса, којим смо се само летимично бавили код расејача. Јасан одговор на питање „како одабрати МоМ зону и број група да бисмо за најкраће време, или са најмањим бројем макро функција базиса, дошли до жељеног резидијума?“ тешко је пружити (ако уопште постоји). Претходна анализа говори у прилог постојању оптималног избора за задати проблем, али нама је одговор потребан без (МоМ) анализе, коју нећемо ни вршити. Имајући у виду да све анализе говоре да је поправка у првој итерацији (у односу на почетно решење) највећа и да за веома велике проблеме неће ни бити времена за другу итерацију (а ова метода је управо и предвиђена за веома велике проблеме), умесно делује покушај да оптимизујемо прву итерацију. Величину МоМ зоне треба одабрати тако да буде што већа, а да при томе МоМ анализа делимичног модела (тј. МоМ зоне) буде временски и меморијски ефикасна. За наш проблем, λ4 модел задовољава у том погледу. Повећање броја група повећава нагиб криве у првој итерацији. Да бисмо видели и колико, посматрамо резидијум након прве итерације за λ4 модел са различитим бројем група, приказан на слици 6.12. Овај нагиб ће, и очито и здраворазумски (јер, када се број група приближава броју МоМ функција базиса, а то је овде 95.730, резидијум тежи нули) постајати све стрмији, али то нама није битно, јер смо са 10.064 групе на граници прихватљивог – PDM анализа која нема макар за ред величине мање макро функција базиса од МоМ функција базиса нема претераног смисла. С друге стране, нумерички експерименти указују да минималан број група за ефикасну PDM анализу не треба да буде мањи од квадратног корена из броја МоМ функција базиса - овде модел са 342 групе представља ту доњу међу. Између модела са 342 групе и оног са 10.064 групе број макро функција базиса повећао се 20 пута, а резидијум се смањио мање од 3 пута. Дакле, велико повећање броја макро функција базиса зарад малог смањења резидијума. Време анализе се при томе повећало око 4 пута, што није критично ако је краће од та два времена око пола сата, као овде, али се ситуација мења ако је оно пар дана. 117 Као практично правило, за моделе са преко 000.100>N МоМ функција базиса, можемо да минимални број група M усвојимо као квадратни корен из ,N NM ≥ , а минимални полупречник сфере МоМ зоне, r , као једну таласну дужину, λ≥r . Стварну вредност усвојићемо у зависности од задатог модела, расположивих рачунарских средстава и расположивог времена. Слика 6.12. Резидијум решења добијеног у првој итерацији PDM анализе модела са слике 6.1а за PDM моделе са различитим бројем група (тј. за различит број макро функција базиса). 6.2 Полуталасна дипол антена на авиону Нос авиона је у Декартовом координатном систему усмерен дуж x осе, а z оса нормална је на подножје и усмерена ка крову авиона. Структура је сачињена од 13.072 плочице, а авион је на радној учестаности дуг 60,8 таласних дужина. Полуталасна дипол антена удаљена је од подножја авиона за четвртину таласне дужине, као на слици 6.13. Усвојен је други ред апроксимације за струју на свакој плочици што резултује са укупно 104.279 МоМ функција базиса. 118 Слика 6.13. Полуталасна дипол антена удаљена је четвртину таласне дужине од подножја авиона. Авион је сачињен од 13.072 плочице и на радној учестаности је дуг 60,8 таласних дужина. Усвојен је други ред апроксимације за струју на свакој плочици што резултује са укупно 104.279 МоМ функција базиса. а) б) в) г) Слика 6.14. (а) λ2 модел - 319 група плочица у PO зони, МоМ зона у сфери полупречника λ2 , (б) λ4 модел - 316 група плочица у PO зони, МоМ зона у сфери полупречника λ4 , (в) λ8 модел - 306 група плочица у PO зони, МоМ зона у елипсоиду полупречника λ4 по z оси и λ8 по преостале две осе, (г) елипсоид МоМ зоне за λ8 модел. 119 Направљена су три PDM модела са различитим величинама МоМ зоне – на слици 6.14а је модел са 319 група плочица у PO зони, чија је МоМ зона ограничена сфером полупречника две таласне дужине, на слици 6.14б је модел са 316 група, чија је МоМ зона ограничена сфером полупречника четири таласне дужине, а на слици 6.14в је модел са 306 група, чија је МоМ зона ограничена елипсоидом полупречницима четири таласне дужине дуж z осе и осам таласних дужина дуж преостале две осе, као што је приказано на слици 6.14г. На сликама 6.14а-в сиво су обојене плочице МоМ зоне. Груписање плочица извршено је запреминским поступком. Слика 6.15. Зависност резидијума од броја итерација за PDM моделе са слика 6.14а-в. Зависност резидијума од броја итерација за моделе приказане на сликама 6.14а-в приказана је на слици 6.15. Иако је број група на доњој граници према претходно усвојеном услову ( NM ≥ ), конвергенција је веома добра. Добијени резидијуми су за три реда величине нижи у односу на претходни пример печ антене. Такође, резидијум је знатно осетљивији на промену величине МоМ зоне. Другим речима, најбитнији део структуре је у околини дипола, и већ МоМ 120 анализа делимичног модела успева да обезбеди добро почетно решење. Увећањем МоМ зоне почетно решење се додатно побољшава, те је у том случају и конвергенција нешто спорија (тешко је битно поправити нешто што је већ добро), односно криве за резидијум нису само „спуштене“ као код модела са печ антеном. а) б) в) г) д) ђ) Слика 6.16. Појачање у симетралној равни авиона ( 0=φ ) добијено на основу решења (а) МоМ анализе делимичног λ2 модела, (б) прве итерације PDM анализе λ2 модела, (в) МоМ анализе делимичног λ4 модела, (г) прве итерације PDM анализе λ4 модела, (д) МоМ анализе делимичног λ8 модела, и (ђ) прве итерације PDM анализе λ8 модела. 121 Криве појачања приказане су на слици 6.16. Видимо оно што се дало закључити и посматрањем кривих за резидијум – повећање МоМ зоне значајно унапређује почетно решење („нулту“ итерацију), па тако λ8 модел већ у нултој итерацији даје сасвим добро решење, очито боље него λ2 модел у првој итерацији. 6.3 Планарни низ микрострип печ антена монтираних на хеликоптер Планарни низ од 44× микрострип печ антене монтиран је на хеликоптер, као на слици 6.17. Нос хеликоптера је у Декартовом координатном систему усмерен дуж x осе, а z оса нормална је на подножје и усмерена ка елиси на крову. Структура је сачињена од 84.773 плочице, а хеликоптер је на радној учестаности дуг 293,3 таласне дужине. Диелектрик печ антене моделован је као савршен, релативне пермитивности 3,38. Струја на плочицама представљена је другим и трећим редом апроксимације што резултује са укупно 1.489.732 МоМ функција базиса. Дакле, безмало милион и по. Слика 6.17. Планарни низ од 44× микрострип печ антене, десно, монтиран је на подножје хеликоптер. PDM модел, лево, има МоМ зону облика сфере полупречника 2,4 таласне дужине са центром у средишту планарног низа, и 217 група плочица у PO зони. У МоМ зону уврштене су плочице које су макар једним делом у сфери, полупречника 2,4 таласне дужине, чији је центар у средишту планарног низа. 122 Плочице из PO зоне подељене су, запреминским поступком груписања, на 217 група. PDM модел приказан је на слици 6.17, лево. Овај модел добили смо љубазношћу фирме WIPL-D d.o.o, где је PDM анализа извршена на рачунару са процесором Intel Core i7 930 на 2.8GHz, 24 GB RAM меморије, осам дискова капацитета по 2 ТB и три графичке картице GeForce GTX 480, свака са по 1536 MB VRAM меморије. Оперативни систем је Win 7 Pro 64-bit. За прву итерација PDM анализе било је потребно око 39 сати. Резидијум решења у првој итерацији је 0,0123. Тродимензионални дијаграм зрачења добијен на основу решења прве PDM итерације приказан је на слици 6.18. Слика 6.18. Тродимензионални дијаграм зрачења добијен на основу решења прве PDM итерације. 123 7. Закључак Кључни допринос ове дисертације је нова метода за итеративно решавање површинских интегралних једначина електромагнетског поља. Метода је названа метода момената вођена физичком оптиком (скраћено PDM, по енглеском називу). PDM је надградња методе момената (МоМ-а) примењене на површинске интегралне једначине поља, од које наслеђује читаву поставку: геометријско моделовање, апроксимацију струја (функције базиса) и елементе МоМ матрице система. Оперативни циљ обе методе је исти: израчунавање тежинских коефицијената функција базиса (тј. расподеле струје). Разлика је у начину њиховог израчунавања – PDM заобилази „инверзију“ МоМ матрице (операцију која ограничава електричну величину проблема који се могу ефикасно решити МоМ-ом) и тежинске коефицијенте израчунава итеративним поступком у коме се правац корекције одређује на оригиналан начин (различит нпр. од познатих итеративних поступака заснованих на Криловљевим просторима). У оквиру дисертације развијен је математички алат за решавање две врсте електромагнетских проблема: 1. расејача начињених од савршено проводног материјала, и 2. антена на платформама, начињених комбинацијом савршено проводног материјала и линеарних изотропних диелектрика. У оба случаја анализа се врши у комплексном домену, на једној учестаности и за једну побуду. Модел расејача мора бити физички затворена тродимензионална савршено проводна структура, док за антенске платформе то није обавезан услов за функционисање (али јесте за ефикасну анализу). Тачност PDM решења у односу на стандардно МоМ решење описујемо резидијумом. За решење идентично ономе добијеном МоМ-ом резидијум је 0, ако за решење усвојимо све нуле, резидијум је 1; за смислено, али не сасвим тачно, решење очекујемо резидијум много мањи од 1. 124 На брзину конвергенције PDM итеративног поступка можемо да утичемо преко броја група МоМ функција базиса ( M ) и начина груписања, а код антенских проблема и преко величине МоМ зоне. PDM се може се применити на било које постојеће решење базирано на МоМ- у. Детаљи реализације зависе од МоМ решења, а за резултате приказане у овој дисертацији користили смо нумеричко језгро софтверског пакета за електромагнетско моделовање WIPL-D Pro [75] (које користи површинску формулацију интегралних једначина електричног поља, геометријско моделовање билинеарним површима и конусним жицама, апроксимацију струја коришћењем полиномских функција, Галеркиново тестирање, те решавање система линеарних једначина паралелизованом LU-декомпозицијом). WIPL-D језгро је у неизмењеном облику искоришћено за тополошку анализу модела, рачунање чланова оригиналне МоМ матрице и рачунање магнетског поља – тим аспектима се нисмо бавили у дисертацији. Написан је програмски код за остале процедуре, па је све уклопљено у програмску компоненту овог софтверског пакета. Све приказане анализе урађене су коришћењем графичког интерфејса и нумеричког језгра те програмске компоненте. Примена PDM-а на расејаче показала је да итеративни поступак релативно брзо конвергира, али и да брзина конвергенције опада са порастом броја итерација. Ово тумачимо смањеном количином „информација“ које носе корекционе струје у каснијим итерацијама. Математички гледано, простор решења, који „обогаћујемо“ новим корекционим струјама, добија све мање ортогоналних компоненти. Са повећањем броја група МоМ функција базиса конвергенција се убрзава, но тада расте и утрошак меморије, па постоји горња граница за прихватљиви број група (при којем је PDM меморијски значајно ефикаснији од МоМ-а). Док год се држимо препоруке да број група буде сразмеран квадратном корену из броја МоМ функција базиса, меморијски утрошак је сразмеран броју МоМ функција базиса, односно није критичан. Такође смо видели да је груписање физички блиских функција базиса ваљано решење, а да се варијацијама у начину овог груписања само незнатно може добити на тачности. Резидијум код каноничких примера, сфере и коцке, је за ред величине мањи него резидијум за реалистичне примере авиона и хеликоптера. Ово је 125 очекивано, с обзиром да је итеративна корекција решења коришћењем модификованог PO-а погоднија за глатке и равне површи. Такође очекивано, мањи резидијум значи и мање одступање у резултатима за радарски одзив и блиско поље. Потребно рачунарско време са порастом броја итерација расте линеарно код модела са мањим бројем група, док се при већем броју група уочава нелинеарно понашање, у складу са теоријским предвиђањима. Конвергенција PDM методе код антена умногоме зависи од „спреге“ између МоМ зоне и PO зоне. Ако је спрега слаба, конвергенција је у почетним итерацијама веома брза (али исто тако брзо опада са порастом броја итерација), а почетно решење веома осетљиво на величину МоМ зоне. Ако су МоМ зона и PO зона „јако“ спрегнуте, почетно решење не мења се драстично са променом величине МоМ зоне, а конвергенција је значајно спорија него код слабе спреге. Наравно, више група МоМ функција базиса значи и бржу конвергенцију. Треба имати у виду и да је резидијум, као мера блискости PDM решења са МоМ решењем, мање индикативан код антена него код расејача (јер је примеренији ситуацијама у којима нема драматичних промена у расподели површинских струја, што је случај код расејача, али не и код антена, где зона напајања може имати значајно веће струје од остатка модела). Избор параметара PDM анализе је питање компромиса. Практично, број група не би требало да је мањи од корена из броја МоМ функција базиса, а полупречник сфере која описује МоМ зону (код антена на платформама) не би требало да је мањи од једне таласне дужине. Горње границе ових величина ограничене су расположивим рачунарским средствима. Резултат сваке PDM итерације је апроксимација одговарајућег МоМ решења. Грубо говорећи, PDM решење са резидијумом од 0,01 обезбеђује добру апроксимацију МоМ решења, а при резидијуму од 0,001 разлика између PDM и МоМ решења тешко се уочава. Снага PDM анализе лежи у ефикасном решавању електрички великих проблема, у којима МоМ дискретизација резултује системом са великим бројем једначина (у овој дисертацији електрички великим проблемима сматрали смо оне са преко 100.000 једначина). За PDM анализу довољно је у једном тренутку на 126 располагању имати само један ред МоМ матрице, а сама PDM матрица је много мањих димензија МоМ матрице. Зато при PDM анализи електрички великих проблема добијамо прихватљиву апроксимацију МоМ решења са више редова мањим меморијским заузећем и знатно краћим временом анализе. У дисертацији смо показали решење добијено након прве итерације PDM анализе модела од милион и по непознатих на „кућном“ рачунару, која је трајала мање од два дана, уз минималан меморијски утрошак. (Меморијско заузеће је практично у свим примерима занемарљиво у односу на одговарајуће МоМ решавање.) Процењено време трајања класичне МоМ анализе на истом рачунару, уз максималну GPU паралелизацију је двадесет дана, док је потребно меморијско заузеће преко 33 TB. За проблеме ове величине не очекујемо да се у пракси ради више од једне итерације. То је сасвим довољно да решење добије на детаљима које почетно решење не може да „ухвати“. Поправка у другој итерацији најчешће није довољна да би оправдала (предуго) рачунарско време потребно за њено извршавање. PDM итеративни поступак начелно је нумерички стабилан. У свакој итерацији се минимизира резидијум, а постоји вишe макро функција базиса него у претходној (итерацији), па очекујемо сигурно смањивање резидијума. До појаве линеарно зависних једначина у PDM систему (и краха при његовом решавању) може доћи само у случају да се у две различите итерације над истом групом креирају потпуно идентичне макро функција базиса, што је у нормалним околностима практично немогуће (чак и мала разлика у њиховим тежинским коефицијентима омогућила би нормалан прорачун). Ипак, у критичним ситуацијама (које су излазиле на видело током тестирања) показује се предност приказа бројева у „двострукој тачности“ (8 бајтова по реалном броју) уместо у „једнострукој“ (4 бајта по реалном броју). Полиномске функције базиса, које смо користили, дају слабије условљене матрице (са већим condition number) при вишим степенима полинома. Ово условљава опадање тачности (повећање резидијума) када се користе виши степени полинома (виши редови апроксимације). Томе у овом раду није посвећена пажња, јер је нагласак на општости методе. Сви приказани нумерички експерименти урађени су са претежено другим редом апроксимације, изузев 127 завршног примера од милион и по непознатих, где је апроксимација готово искључиво трећег реда. У даљим истраживањима, пажња ће бити посвећена побољшању конвергенције, убрзању креирања макро функција базиса и побољшању нумеричке стабилности PDM анализе при коришћењу МоМ функција базиса вишег реда. 128 Литература [1] Maxwell, J., “A Dynamical Theory of the Electromagnetic field,” Scientific Papers, pp. 459-512, 1864. [2] Whittaker, E.T., A History of the Theories of Aether and Electricity, Longmans, Green, and Co., London, 1910. [3] Darrigol, O., Electrodynamics from Ampere to Einstein, Oxford University Press, New York, 2000. [4] Maxwell, J., A Treatise on Electricity and Magnetism, Macmillan and Co. Ltd., London, 1873. [5] Heaviside, O., Electromagnetic Theory, Vol. I, The Electrician - Printing and Publishing Companu Ltd., London, 1893. [6] Hertz, H., Miscellaneous Papers, (Authorised English Translation by Jones, D.E.) Macmillan and Co. Ltd., London, 1896. [7] Stratton, J.A., Electromagnetic Theory, McGraw Hill, New York, 1941. [8] Jackson, J.D., Classical Electrodynamics, John Wiley&Sons, Inc., New York, 1999. [9] Surutka, J., Elektromagnetika, Građevinska knjiga, Beograd, 1965. [10] Popović, B.D., Elektromagnetika, Građevinska knjiga, Beograd, 1990. [11] Đorđević, A., Elektromagnetika, Akademska misao, Beograd, 2008. [12] Nersessian, N.J., “Constructing Meaning in Scientific Theories”, Kluwer, 1990. [13] Harrington, R.F, Time-Harmonic Electromagnetic Fields, New York, McGraw- Hill Book Co., 1961. [14] Harrington, R.F., Field Computation by Moment Methods, New Yorк, The MacMillian Co., 1968. [15] Harrington, R.F., “Origin and development of the method of moments for field computation,” Antennas and Propagation Magazine, IEEE , vol.32, no.3, pp.31- 35, Jun 1990. [16] Peterson, A.F., Ray, S;L., and Mittra, R., Computational Methods For Electromagnetics, IEEE Press, New York, 1998. 129 [17] Sarkar, T.K, Djordjevic, A.R., and Arvas, E., “On the Choice of Expansion and Weighting Functions in the Numerical Solution of Operator Equations,” Antennas and Propagation, IEEE Transactions on , vol.33, no.9, pp.988-996, Sept. 1985. [18] Djordjevic, A.R. and Sarkar, T.K, “A Theorem on the Moment Methods,” Antennas and Propagation, IEEE Transactions on , vol.35, no.3, pp.988-996, March 1987. [19] Abramowitz, M and Stegun. I, Handbook of Mathematical Functions, National Bureau of Standards, 1964. [20] Arfken, G.B. and Weber, H.J., Mathematical Methods for Physicists, Elsevier Inc., USA, 2005. [21] Golub, G.H. and Van Laon, C.F., Matrix Computations, The Johns Hopkins University Press, Baltimore, 1996. [22] Parlett, B.N., The Symmetric Eigenvalue Problem, SIAM, USA, 1998. [23] Press, W.H., Teukolsky, S.A., Vetterling, W.T., and Flannery, B.P., Numerical Recipies, Cambridge University Press, USA, 2007. [24] Saad, Y., Iterative Methods for Sparse Linear Systems, SIAM, 2003. [25] Strang, G., Computational Science and Engineering, Wellesley-Cambridge Press, USA, 2007. [26] Strang, G., Introduction to Linear Algebra, Wellesley-Cambridge Press, USA, 2009. [27] Trefethen, L. and Bau. D., Numerical Linear Algebra, SIAM, USA, 1997. [28] Knepp, D.L. and Goldhirsh, J., “Numerical analysis of electromagnetic radiation properties of smooth conducting bodies of arbitrary shape,” IEEE Trans. Antennas Propag., vol. AP-20, no. 3, pp. 383–388, May 1972. [29] Newman, E.H. and Pozar, D.M., “Electromagnetic modeling of composite wire and surface geometries,” IEEE Trans. Antennas Propagat., vol. AP-26, pp. 784– 789, Nov. 1978. [30] Wang, J. and Papanicolopulos, C., “Surface-patch modeling of scatterers of arbitrary shapes,” Antennas and Propagation Society International Symposium, 1979 , vol.17, no., pp. 159- 162, Jun 1979. [31] Singh, J. and Adams, T, “A nonrectangular patch model for scattering from surfaces,” IEEE Trans. Antennas Propagat., vol. AP-27, pp. 531–535, July 1979. [32] Glisson, A.W. and Wilton, D.R., “Simple and efficient numerical methods for problems of electromagnetic radiation and scattering from surfaces,” IEEE Trans. Antennas Propagat., vol. AP-28, pp. 593–603, Sept. 1980. 130 [33] Rao, S., Wilton, D., and Glisson, A., “Electromagnetic scattering by surfaces of arbitrary shape, ” Antennas and Propagation, IEEE Transactions on , vol.30, no.3, pp. 409- 418, May 1982. [34] Newman, E.H. and Tulyathan, P., “A surface patch model for polygonal plates,” IEEE Trans. Antennas Propagat., vol. AP-30, pp. 588–593, July 1982. [35] Newman, E.H., Alexandropoulos, P., and Walton, E.K., “Polygonal plate modeling of realistic structures,” IEEE Trans. Antennas Propagat., vol. AP-32, pp. 742–747, July 1984. [36] Wilkes, D.L. and Cha, C.-C. “Method of moments solution with parametric curved triangular patches,” in IEEE APS Int. Symp. Dig., London, Canada, July 1991, pp. 1512–1515. [37] Kolundzija, B.M. and Popovic, B.D., “Entire-domain Galerkin method for analysis of metallic antennas and scatterers,” Proc. Inst. Elect. Eng., vol. 140, pt. H, pp. 1–10, Feb. 1993. [38] Aberegg, K.R. Taguchi, A., and Peterson, A.F., “Application of higher order vector basis functions to surface integral equation formulations,” Radio Sci., vol. 31, no. 5, pp. 1207– 1213, Sep.-Oct. 1996. [39] Kolundzija, B.M., “Comparison of a class of subdomain and entire domain basis functions automatically satisfying KCL,” Antennas and Propagation, IEEE Transactions on , vol.44, no.10, pp.1362-1366, Oct 1996. [40] Graglia, R.D., Wilton, D.R., and Peterson, A.F., “Higher order interpolatory vector bases for computational electromagnetics,” IEEE Trans. Antennas Propagat., vol. 45, pp. 329-342, Mar. 1997. [41] Notaros, B.M., Popovic, B.D., Weem, J.P., Brown, R.A., and Popovic, Z.B., “Efficient large-domain MoM solution to electrically large EM problems,” IEEE Trans. Microwave Theory and Techniques, vol. 49, pp. 151–159, Jan. 2001. [42] Gang Kang; Jiming Song; Weng Cho Chew; Donepudi, K.C.; Jian-Ming Jin; , "A novel grid-robust higher order vector basis function for the method of moments," Antennas and Propagation, IEEE Transactions on , vol.49, no.6, pp.908-915, Jun 2001. [43] Cai, W. Yu, T., Wang, H., and Yu, Y. “High-order mixed RWG basis functions for electromagnetic applications,” IEEE Trans. Microwave Theory and Techniques, vol. 49, pp. 1295–1303, July 2001. [44] Notaros, B.M., “Higher Order Frequency - Domain Computational Electromagnetics,” Antennas and Propagation, IEEE Transactions on , vol.56, no.8, pp. 2251–2276, August 2008. [45] Zoric, D.P.; Olcan, D.I.; Kolundzija, B.M.; , “Solving electrically large EM problems by using out-of-core solver accelerated with multiple graphical processing units,” Antennas and Propagation (APSURSI), 2011 IEEE International Symposium on , vol., no., pp.1-4, 3-8 July 2011. 131 [46] Peterson, A.F. and Mittra, R., “Method of conjugate gradients for the numerical solution of large-body electromagnetic scattering problems,” J. Opt. Soc. Amer., vol. 2, pp. 971–977, 1985. [47] Sarkar, T., “The conjugate gradient method as applied to electromagnetic field problems,” Antennas and Propagation Society Newsletter, IEEE , vol.28, no.4, pp.4-14, August 1986. [48] Peterson, A.F., Smith, C.F., and Mittra, R., “Eigenvalues of the moment-method matrix and their effect on the convergence of the conjugate gradient algorithm [EM scattering],” Antennas and Propagation, IEEE Transactions on , vol.36, no.8, pp.1177-1179, Aug 1988. [49] Ray, S.L. and Peterson, A.F., “Error and convergence in numerical implementations of the conjugate gradient method [EM problems],” Antennas and Propagation, IEEE Transactions on , vol.36, no.12, pp.1824-1827, Dec 1988. [50] Smith, C.F., Peterson, A.F., and Mittra, R., “The biconjugate gradient method for electromagnetic scattering,” Antennas and Propagation, IEEE Transactions on , vol.38, no.6, pp.938-940, Jun 1990 [51] Pocock, M.D. and Walker, S.P., “The complex bi-conjugate gradient solver applied to large electromagnetic scattering problems, computational costs, and cost scalings,” Antennas and Propagation, IEEE Transactions on , vol.45, no.1, pp.140-146, Jan 1997. [52] Canning, F.X., “The impedance matrix localization method (IML) for moment- method calculations,” IEEE Antennas and Propagation Society Magazine, vol. 32, pp. 18-30, Oct. 1990. [53] Canning, F.X., “Improved impedance matrix localization method,” Antennas and Propagation, IEEE Transactions on , vol.41, no.5, pp.659-667, May 1993. [54] Jorgensen, E., Volakis, J., Meincke, P., and Breinbjerg, O., “Higher order hierarchical Legendre basis functions for electromagnetic modeling,” Antennas and Propagation, IEEE Transactions on , vol.52, no.11, pp. 2985- 2995, Nov. 2004. [55] Coifman, R., Rokhlin, V., and Wandzura, S., “The fast multipole method for the wave equation: A pedestrian prescription,” IEEE Antennas Propagat. Mag., vol. 35, no. 3, pp. 7–12, June 1993. [56] Song, J., Cai-Cheng Lu, and Weng Cho Chew, “Multilevel fast multipole algorithm for electromagnetic scattering by large complex objects,” Antennas and Propagation, IEEE Transactions on , vol.45, no.10, pp.1488-1493, Oct 1997. [57] Donepudi, K.C., Song, J., Jin, J.-M., Kang, G., and Chew, W.C.; , “A novel implementation of multilevel fast multipole algorithm for higher order Galerkin's method,” Antennas and Propagation, IEEE Transactions on , vol.48, no.8, pp. 1192- 1197, Aug 2000. 132 [58] Ergul, O. and Gurel, L., “Efficient Parallelization of the Multilevel Fast Multipole Algorithm for the Solution of Large-Scale Scattering Problems,” Antennas and Propagation, IEEE Transactions on , vol.56, no.8, pp.2335-2345, Aug. 2008. [59] Hesford, A.J., Chew, W.C., “On Preconditioning and the Eigensystems of Electromagnetic Radiation Problems,” Antennas and Propagation, IEEE Transactions on , vol.56, no.8, pp.2413-2420, Aug. 2008. [60] Canning, F.X. and Rogovin, K., “Fast direct solution of standard moment-method matrices,” Antennas and Propagation Magazine, IEEE , vol.40, no.3, pp.15-26, Jun 1998. [61] Shaeffer, J., “Direct Solve of Electrically Large Integral Equations for Problem Sizes to 1 M Unknowns,” Antennas and Propagation, IEEE Transactions on , vol.56, no.8, pp.2306-2313, Aug. 2008. [62] Umashankar, K.R., Nimmagadda, S.; and Taflove, A., “Numerical analysis of electromagnetic scattering by electrically large objects using spatial decomposition technique,” Antennas and Propagation, IEEE Transactions on , vol.40, no.8, pp.867-877, Aug 1992. [63] Olćan, D.I., Stevanović, I.M., Kolundžija, B.M., Mosig, J.R., Djordjević, A.R., “Diakoptic Surface Integral-Equation Formulation Applied to 3-D Scattering Problems,” 24th Annual Review of Progress in Applied Computational Electromagnetics (ACES), Niagara Falls, Canada, pp. 676-681, March 30 - April 4 2008. [64] Matekovits, L., Vecchi, G., Dassano, G., and Orefice, M., “Synthetic function analysis of large printed structures: The solution space sampling approach,” Proc. IEEE AP-S 2001, Boston, USA, pp. 568–571, Jul. 2001. [65] Matekovits, L., Laza, V.A., and Vecchi, G., “Analysis of Large Complex Structures With the Synthetic-Functions Approach,” Antennas and Propagation, IEEE Transactions on , vol.55, no.9, pp.2509-2521, Sept. 2007. [66] Prakash, V.V.S. and Mittra, R., “Characteristic basis function method: A new technique for efficient solution of Method of Moments matrix equations,” Microwave and Optical Technology Letters, Vol. 36, No. 2, pp. 95-100, Jan. 2003. [67] Ufimtsev, P., “New Insight into the Classical Macdonald Physical Optics Approximation,” Antennas and Propagation Magazine, IEEE , vol.50, no.3, pp.11-20, June 2008. [68] Thiele, G. and Newhouse, T., “A hybrid technique for combining moment methods with the geometrical theory of diffraction,” Antennas and Propagation, IEEE Transactions on, vol.23, no.1, pp. 62–69, Jan 1975. [69] Tae Kim and Thiele, G., “A hybrid diffraction technique--General theory and applications,” Antennas and Propagation, IEEE Transactions on , vol.30, no.5, pp. 888–897, Sep 1982. 133 [70] Al-Hekail, Z.O. and Burnside, W.D., “A hybrid method for computing the scattered fields from complex structures,” Antennas and Propagation, IEEE Transactions on , vol.42, no.6, pp.847-852, Jun 1994. [71] Jakobus, U. Landstorfer, F.M., “Improved PO-MM hybrid formulation for scattering from three-dimensional perfectly conducting bodies of arbitrary shape,” IEEE Trans. Antennas Propag., vol. 43, no. 2, pp. 162–169, Feb. 1995. [72] Jorgensen, E., Meincke, P., and Breinbjerg, O., “A hybrid PO-higher-order hierarchical MoM formulation using curvilinear geometry modeling,” Antennas and Propagation Society International Symposium, 2003. IEEE , vol.4, no., pp. 98–101 vol.4, 22-27 June 2003. [73] Djordjevic, M., and Notaros, B.M., “Higher order hybrid method of moments- physical optics modeling technique for radiation and scattering from large perfectly conducting surfaces,” Antennas and Propagation, IEEE Transactions on, vol.53, no.2, pp. 800–813, Feb. 2005. [74] Kolundžija, B.M. and Djordjević, A.R., Electromagnetic modeling of composite metallic and dielectric structures, Boston, Artech House, 2002. [75] WIPL-D Pro v9.0, 3D EM solver, http://www.wipl-d.com, 2011. [76] Tasic, M.S. and Kolundzija, B.M., “A PO driven iterative solution of MFIE for large scatterers,” Telecommunications in Modern Satellite, Cable and Broadcasting Services, 2005. 7th International Conference on , vol.1, no., pp. 24- 27 vol. 1, 28-30 Sept. 2005. [77] Tasic, M.S. and Kolundzija, B.M., “PO Driven Iterative Least Square Solution of MFIE,” Proc. of IEEE/ACES Conf. on Wireless Comm. and Appl. Commput. Electromag., pp. 470-475, (CD ROM Edition: s15p06.pdf), Miami, Florida, 12-16 March 2006. [78] Tasic, M.S. and Kolundzija, B.M., “PO Driven Iterative Galerkin Solution of Field Integral Equations,” Antennas and Propagation Society International Symposium 2006, IEEE , vol., no., pp.4073-4076, 9-14 July 2006. [79] Tasić, M., Kolundžija, B., “Comparison of PO Driven and Conjugate Gradient Iterative Solution of Field Integral Equations,” ETRAN 2006, Beograd, Srbija, 6-8 jun 2006. [80] Tasic, M.S. and Kolundzija, B.M., “Efficient Analysis of Large Scatterers by Physical Optics Driven Method of Moments,” Antennas and Propagation, IEEE Transactions on , vol.59, no.8, pp.2905-2915, Aug. 2011. [81] Chen, K.-M., “A mathematical formulation of the equivalence principle,” Microwave Theory and Techniques, IEEE Transactions on , vol.37, no.10, pp.1576-1581, Oct 1989. [82] Tošić, D.Đ., Uvod u numeričku analizu, Naučna knjiga, Beograd, 1991. БИОГРАФИЈА Миодраг Тасић рођен је 1972. у Београду. Дипломирао је на Електротехничком факултету у Београду (просечна оцена 8,56) 21. септембра 1998. са темом "Решавање дводимензионалних електростатичких проблема методом коначних елемената", одбрањеним на Катедри за општу електротехнику. Магистарску тезу "Ефикасно електромагнетско моделовање засновано на аутоматској сегментацији полигона на четвороуглове" одбранио је 22. октобра 2004. године на Смеру за примењену електромагнетику и оптоелектронику. Септембра 2000. године изабран је за асистента-приправника при Катедри за општу електротехнику Електротехничког факултета у Београду, а септембра 2005. године изабран је за асистента при истој Катедри. Држи вежбе на табли из предмета Основи електротехнике, Електромагнетика, Антене и простирање и Микроталасна мерења, као и лабораторијске вежбе из антенских и микроталасних мерења. Главне области истраживања су му електромагнетско моделовање и нумеричка електромагнетика. Коаутор је три рада у часописима и више конференцијских радова у области примењене електромагнетике. Коаутор је два софтверска пакета за електромагнетско моделовање и анализу. Учесник је четири научно- истраживачка пројекта под покровитељством Министарства за науку и технолошки развој.