ELEKTROMANYETĐK DALGALARIN BĐYOLOJĐK DOKULAR
ĐÇĐNDEKĐ YAYILIMINININ ZAMANDA SONLU FARKLAR
METODU ĐLE ANALĐZĐ VE SĐMÜLASYONU
Cafer BUDAK
YÜKSEK LĐSANS TEZĐ
(ELEKTRĐK-ELEKTRONĐK MÜHENDĐSLĐĞĐ ANABĐLĐM DALI)
DĐYARBAKIR AĞUSTOS - 2008
T.C.
DĐCLE ÜNĐVERSĐTESĐ
FEN BĐLĐMLERĐ ENSTĐTÜSÜ MÜDÜRLÜĞÜ DĐYARBAKIR
Cafer BUDAK tarafından yapılan " Elektromanyetik Dalgaların Biyolojik Dokular Đçindeki Yayılımınının Zamanda Sonlu Farklar Metodu Đle Analizi Ve
Simülasyonu " konulu bu çalışma, jürimiz tarafından Elektrik-Elektronik Mühendisliği Anabilim Dalında YÜKSEK LĐSANS tezi olarak kabul
edilmiştir.
Jüri Üyesinin Ünvanı Adı Soyadı Başkan: Prof.Dr. Süleyman DAŞDAĞ Üye : Doç. Dr. Mehmet AKIN
Üye : Yrd. Doç. Dr. M. Bahattin KURT
Tez Savunma Sınavı Tarihi : .../.../...
Yukarıdaki bilgilerin doğruluğunu onaylarım. .../.../...
Prof. Dr. HAMDĐ TEMEL ENSTĐTÜ MÜDÜRÜ
( MÜHÜR )
TEŞEKKÜR
Đlk olarak şunu ifade etmek isterim ki, bu çalışmanın başlatılması, geliştirilmesi ve tamamlanmasında gerekli sabrı, rehberliği, tavsiyeleri ve desteğini esirgemeyen değerli danışmanım Sayın Yrd. Doç Dr. M. Bahattin Kurt’a teşekkürlerimi sunarım.
Çalışmamın sürdürülmesinde her zaman yanımda olan değerli Arş. Gör. Hüseyin Acar, Arş. Gör. Ferhat Çıra, Ögr. Gör. Musa Yılmaz ve Ögr. Gör. M. Emin Asker’e teşekkür ederim.
Yüksek lisans eğitimim sırasında gösterdiği sürekli sabır, sevgi ve özverisiyle her zaman yanımda olan beni her konuda destekleyen eşim ve çocuklarıma teşekkürlerimi sunarım.
AMAÇ
Günlük yaşantımıza son on yılda giren cihazların sayısı şu ana kadar bulunan cihazların sayısından fazladır. Bu durum da her zamankinden daha fazla kontrolsüz elektromanyetik alan demektir. Bu artışla birlikte, kullanılan bu kaynakları insan vücudunu ne şekilde etkilediği soruları akla gelmiştir. Özellikle cep telefonları ve diğer mobil kaynaklar başta akla gelenlerdir.
Bu tez çalışmasında en kullanışlı sayısal yöntemlerden biri olan Finite Difference Time domain (FDTD) ayrıntılı olarak incelenmiştir. Elektromanyetik dalga yayınım problemleri için farklı ortamlarda dalga hareketlerinin bir iki ve üç boyutlu matlab programı kullanılarak simülasyonları yapılmıştır.
Üç boyutlu matlab kodları geliştirildikten sonra hesap uzayı içine kafa eşdeğeri dielektrik madde konularak düzlemsel bir dalganın bu yapı içinde EM alan enerjisinin dokularda yarattığı etkiler için Özgül Soğurma Oranı (SAR:Specific Absorption Rate) parametresi hesaplanmıştır.
Geliştirilen matlab kodlar sayesinde, profesyonel olarak piyasada satılan programlara alternatif geliştirilmeye çalışılmıştır. Mevcut kodlar daha da geliştirilebilir ve bunlara ara yüz ve grafik özellikleri eklendiği takdirde diğer profesyonel yazılımlardan hiçbir farkı olmayacaktır. Programların esnekliği sayesinde yapı içerisine daha kompleks modeller eklenip değişik sonuçlar gözlemlenebilir.
ÖZET
Bu çalışmada, güçlü sayısal yöntemlerden biri olan FDTD (Zamanda Sonlu Farklar Yöntemi) ayrıntılı olarak incelenmiştir. FDTD yönteminin hızlı bilgisayarlarda kullanılması, hafızanın artması hesaplama performansının gelişmesi FDTD kullanarak dalga simülasyonunu mümkün kıldı.
Başlangıçta FDTD ile ilgili matematiksel ve fiziksel özellikler verilerek konu ile ilgili teorik bilgi eksikliği giderilmiş. Ardından sayısal analiz için gerekli formülasyonlar oluşturulmuş. FDTD yöntemiyle ilgili sınır koşulları, parametre seçimi, kararlılık kriterleri hakkında bilgi verilmiştir.
Uygulama olarak da, Elektromanyetik dalga yayılım problemleri için geliştirilen, matlab programları ile farklı ortamlarda, farklı dalga tipleri ile dalga hareketlerinin analizleri yapılmış ve bunlarla ilgili bir iki ve üç boyutlu simülasyonlar gerçekleştirilmiştir. Son olarak da küp şeklinde kafa eşdeğeri dielektrik madde üzerindeki emilen EM enerji yani SAR grafikleri çizdirilmiştir.
SUMMARY
In this study, the FDTD (Finite-Difference Time-Domain), one of the strong numerical methods was examined. By using the FDTD method with fast computers, memory increment, computing performance, has offered FDTD usage for wave simulations.
At the beginning, mathematical and physical properties related to FDTD were given to remove theoretical knowledge deficiencies. Then the required formulations were presented for numerical analysis. Knowledge about boundary conditions, parameter selection, and stability criterion for FDTD was presented.
For electromagnetic waves diffusion problems, one, two, and three dimensional MATLAB applications were prepared for analysis of different waves and their movement at various positions. As a result, SAR illustrations about EM energy absorbed from dielectric substances which was equivalent to cubic head.
ĐÇĐNDEKĐLER Sayfa TEŞEKKÜR……….. ...i AMAÇ...ii ÖZET...iii SUMMARY ...iv ĐÇĐNDEKĐLER...v 1. GĐRĐŞ...1 2. ÖNCEKĐ ÇALIŞMALAR...3 3. MATERYAL VE METOT ...5
3.1. Sayısal Teknikler Ve FDTD Yöntemine Genel Bir Bakış...5
3.1.1. Sayısal Teknikler...5
3.1.2. Đletim Hattı Matrisi (TLM,Transmission Line Matrix) Yöntemi ...5
3.1.3. Parabolik Denklem (PE,Parabolic Equation) yöntemi ...5
3.1.4. Moment,(MoM,Moment metod) Yöntemi ...6
3.2. FDTD Ve Sonlu Farklar Yöntemi ...6
3.2.1.Sonlu farklar Tanımı...6
4. FDTD YÖNTEMĐ...9
4.1. FDTD Yönteminin Kullanım Alanları ...10
4.2. Üç Boyutlu FDTD Yöntemi ...11
4.3. Kararlılık (Courant) Koşulu ...14
4.4. FDTD yönteminde Sayısal Dağılma (Dispersiyon)...15
4.5. Parametre Seçimi...16
5. FDTD YÖNTEMĐNĐN HESAP UZAYINA UYGULANMASI ...19
5.1.FDTD Algoritmasında Uyarma Fonksiyonu ...21
5.2. Gauss Darbesinin Ortama Enjekte Edilmesi ...24
5.3. FDTD Algoritması ...27
5.4. Sınır Koşulları ...29
6. ELEKTROMANYETĐK PROBLEMLERĐN FDTD ĐLE ANALĐZĐ ...31
6.1. Bir Boyutlu Simülasyon ...31
6.2. Bir Boyutta Soğurulma Sınır Şartları ...37
6.3. Dalganın Dielektrik Bir Ortama Çarpması...38
6.4. Gauss Darbesi Yerine Sinüzoidal Kaynak Kullanmak...41
7. ĐKĐ BOYUTLU SĐMÜLASYON ...46
7.1. Farklı Kaynak Kullanarak Đki Boyutlu (2-D) Simülasyon ...49
7.2. Kayıplı Ve Sınır Şartları Kullanarak Đki Boyutlu (2-D) Simülasyonlar ...50
8. ÜÇ BOYUTLU SĐMÜLASYON ...52
9. BĐYOLOJĐK DOKULARDA SAR DEĞERĐNĐN HESAPLANMASI ...56
9.1. Dokulardaki Enerji Emilimini Belirleyen Faktörler...56
9.1.1. Dielektrik Özellikler...56
9.1.2. Doku Geometrisi Ve boyutu...57
9.1.3. Doku Yönelimi Ve Alan Pozisyonu ...57
9.1.4. Alan Frekansı ...57
9.1.5. Kaynak Konumu...58
9.1.6. Maruziyet Ortamı ...58
9.1.7. Zaman ve Şiddet Faktörleri ...58
9.2. Matlab kullanarak Sar incelenmesi...58
9.2.1. Eşdeğer Dokuya 30 Mhz Uygulanmış Sonuçlar...61
9.2.2. Eşdeğer Dokuya 900 Mhz Uygulanmış Sonuçlar...62
9.2.3. Eşdeğer Dokuya 1800 Mhz Uygulanmış Sonuçlar...62
9.2.4. Eşdeğer Dokuya 2200 Mhz Uygulanmış Sonuçlar...63
9.2.5. Eşdeğer Dokuya 3 Ghz Uygulanmış Sonuçlar ...63
9.2.6. Eşdeğer Dokuya 5 Ghz Uygulanmış Sonuçlar ...63
9.3. Kaynağın Farklı Uzaklıklardan Uyarılması...64
9.3.1. 1800 MHz Kaynak Konumu X=9 ...64
9.3.2. 1800 MHz Kaynak Konumu X=14 ...64
9.4. Zaman Adımının Etkisi ...65
9.4.1. 1800 MHZ Time Step=200 ...65 9.4.2. 1800 MHZ Time Step=400 ...65 10. SONUÇLAR VE TARTIŞMA ...66 11. EKLER ...68 12. KAYNAKLAR...126 ŞEKĐLLER LĐSTESĐ ...131 ÇĐZELGELER(TABLOLAR) LĐSTESĐ ...132
SĐMGELER VE SEMBOLLER LĐSTESĐ ...133
1.GĐRĐŞ
Günümüzde teknolojinin ulaşılmaz hızı baş döndürmektedir. Özellikle 2000’li yıllara doğru artan bu hız ile birlikte Elektromanyetik Alan (EM) kaynaklarının çeşitliliği ve sayısı da çığ gibi arttı. Elektronik cihazlar artık yaşamımızın vazgeçilmez baş aktörleri haline gelmeye başladı. Hemen hemen herkesin bir cep telefonu veya kişisel bilgisayara sahip olması v.b cihazları kamuya açık alanlarda kullanması artık günümüzde çok olağan sayılmaktadır. Bu durum da her zamankinden daha fazla Radyo frekans (RF) alanlara insanların maruz kalmasına sebep olmaktadır. Her teknolojik gelişimin avantajlarının yanında bazı dezavantajları da vardır. Bu cihazların yaşamımızı kolaylaştırmakla birlikte bazı sorunları da beraberinde getirdiği kesindir. Peki, maruz kalınan EM alanının miktarı neye bağlıdır? Bunları azaltan veya artıran etkenler varmıdır? uygulanan alanın frekansı, şiddeti, doğrultusu, dalga şekli, canlının pozisyonu, maruz kalma süresi canlının fiziksel büyüklüğü ne derece bunu etkiler? gibi sorular karşımıza çıkmaktadır. Kuşkusuz uğraştığı problemin fiziğini iyi bilen ve güçlü sayısal algoritmalar kullanabilen mühendislerin problemdeki başarısı oldukça yüksek olacaktır. Diğer önemli bir nokta kullanılan sayısal algoritmaların da geçerli oldukları parametreleri ve doğruluk derecelerini iyi bilmek gerektiğidir.
EM alan maruziyetinden kaynaklı alınan doz değerinin, doğrudan insan üzerinde ölçmek hem etik hem de mümkün olmaması nedeniyle ya laboratuar ortamında doku eşdeğeri sıvı ile dolu fantomlarda veya kobay farelerde ya da bilgisayar ortamındaki simülasyonlar kullanılarak tespiti yapılmaktadır.
Simülasyon (modelleme), teorik ya da gerçek fiziksel bir sisteme ait neden-sonuç ilişkilerinin bir bilgisayar modeline yansıtılması tekniğidir. Simülasyon çalışmaları değişik koşullar altında gerçek sisteme ait davranışların bilgisayar modelinde izlenmesini sağlar. Kısacası simülasyon, gerçek hayattaki olayların bilgisayar ortamına aktarılması işlemidir. Bu tanıma uygun olarak, simülasyon çalışmalarında gerçek hayatta karşılaşılması en olası durumlar ele alınmaktadır. Bu amaca uygun olarak yaptığımız çalışmada biyolojik dokuların emdiği EM alan miktarını bilgisayar simülasyonları kullanılarak incelenmiştir.
Bu simülasyonlar literatürde kısaca (FDTD) zamanda sonlu farklar yöntemi ingilizce (Finite Difference Time Domain) baş harfleri kullanılarak yapılmıştır. Simülasyonlar öncelikle tek boyutlu kayıpsız ortamlarda en basitten başlanarak üç boyutlu kayıplı birden fazla dielektrik bulunan ortamlara uygulanmıştır. Hem analitik hem de FDTD çözümlerinden elde edilen denklemler Matlab kullanılarak yazılıma dönüştürülmüş ve simülasyonlar gerçekleştirilmiştir.
2. ÖNCEKĐ ÇALIŞMALAR
1966 Đlk defa YEE [1] tarafından Maxwell’in rotasyonel (curl) denklemleri doğrudan zaman domeninde FDTD sayısal tekniği ile tanımlandı.
1975 Taflove ve Brodvin YEE’nin algoritmaları için doğru sayısal kararlılık kriterini elde ettiler[2,3].
1981 MUR Yee hücreleri için ilk kararlı emici sınır koşullarını(ABC) yayınladı[4].
1982 Taflove ve Umashankar ilk FDTD elektromanyetik dalga saçılımını geliştirdi. Modeller yakın alan uzak alan iki ve üç boyutlarda hesaplandı[5,6].
1986 Choi ve Hoefer ilk FDTD dalga kılavuzu simülasyonunu yayınladı[7].
1987 Kriegsmann ve Moore ABC teorisi üzerine ilk makaleyi yayınladılar. IEEE Trans. Antennas and propagations [8,9].
1988 Sullivian Đnsan vücudu tarafından emilen sabit sinüzoidal dalganın üç boyutlu modelini yayınladı[10].
1988 Zhang tarafından FDTD mikroşerit modeli sunuldu [11].
1990 Frekansa bağlı elektriksel geçirgenlik Kashiwa ve Fukai tarafından sunuldu[12].
1990 Antenlerin FDTD modeli Maloney tarafından sunuldu[13].
1990 Optoelektronik anahtarların FDTD modeli Sano ve Shibata tarafından sunuldu[14].
1991 Doğrusal olmayan dağıtıcı ortamların optik darbelerin yayınımının FDTD modeli sunuldu[15].
1992 Elektronik devre elemanlarının iki boyutlu FDTD modeli Sue tarafından sunuldu[16].
1994 Berenger tarafından Perfect Matched Layer (PML) iki boyutlu olarak sunuldu. Katz tarafından üç boyuta genişletildi[17,18].
3. MATERYAL VE METOT
3.1. Sayısal Teknikler Ve FDTD Yöntemine Genel Bir Bakış
Her ne kadar sayısal yöntemler burada anlattıklarımızla sınırlı değilse de, en güçlü ve en yaygın olarak kullanılanlardır.
3.1.1. Sayısal Teknikler
Elektromanyetik problemlerin çözümünde birçok yöntem bulunmaktadır. Birden çok yönteme gereksinim duyulmasının sebebi her problem için farklı çözüm gereksinimidir. Yani, her yöntem belli koşullarda, belli problemlerde iyi sonuç verirken, farklı problemlerde aynı yöntem istenilen hassas sonuçları verememektedir.
3.1.2. Đletim Hattı Matrisi (TLM, Transmission Line Matrix) Yöntemi
Đletim hattı eşdeğeri yöntemi, yapısal olarak FDTD tekniğinden çok farklı görünse de, benzer mantık ve sayısallaştırma ilkeleriyle zamanda iteratif çözümler üreten diğer bir yöntemdir. FDTD tekniğinde alan denklemleriyle yola çıkılırken TLM tekniğinde devre mantığı kullanılır. Her iki teknikte de ele alınan üç boyutlu yapı birim hücrelerin toplamı şeklinde düşünülür. Yapı yüz binlerce eş hücreye bölündükten sonra FDTD'de alan denklemleri, TLM'de ise iletim hatları denklemleri ile zamanda iteratif çözümler üretilir. FDTD tekniğinde elektrik ve manyetik alan bileşenleri bağımsız değişkenlerdir. Oysa TLM tekniğinde gerilim ve akımlar bağımsız değişken olarak kullanılır[19].
3.1.3. Parabolik Denklem (PE, Parabolic Equation) Yöntemi
PE yöntemi, parabolik formdaki dalga denkleminin Fourier dönüşümü kullanılarak sayısal çözümüne dayanır. Özellikle iki boyutlu (2-D) dalga iletimi problemlerinde çok kullanışlı bir yöntemdir. Bu yöntem hem akustik hem de elektromanyetik dalga
iletimi problemlerine başarı ile uygulanmaktadır. Parabolik formdaki dalga denklemi, tek yönde ilerleyen dalga hareketlerini modellediğinden, PE yöntemi ile geriye yansımaların olmadığı ya da önemsiz sayılabildiği problemler ele alınabilir[19].
3.1.4. Moment (MoM, Method of Moment) Yöntemi
MoM tekniği, özellikle anten tasarımı ve saçılma-ters saçılma problemlerinde yaklaşık otuz yıldır yaygın olarak uygulanan bir frekans domeni çözümüdür. Yöntem ana hatlarıyla Huygens eşdeğer kaynaklar prensibine dayanır. Ele alınan diğer üç yöntemden farklı olarak MoM tekniğinde problemin Green fonksiyonu bulunur. Yani problem aslında analitik olarak çözülür. Ancak Green fonksiyonu bilinen yada sayısal olarak hesaplanabilen bir problemde tersten gidilerek özellikle iletken yüzeyler sonsuz küçük parçalara bölünür ve bu parçalardaki yüzey akımları sayısal olarak çözülür. Bu yüzey akımları hesaplandığında, Green fonksiyonu yardımıyla istenen noktadaki elektromanyetik etki hesaplanır[19].
3.2. FDTD Ve Sonlu Farklar Yöntemi
FDTD yöntemi, diferansiyel formdaki Maxwell denklemlerinin doğrudan zaman domeninde ayrıklaştırılıp çözülmesidir. Maxwell denklemindeki diferansiyel operatörler sayısallaştırılıp sonlu farklar eşdeğerleri ile değiştirilir.
3.2.1. Sonlu Farklar Tanımı
FDTD yöntemi analitik türev operatörünün sayısallaştırılmasına dayanır ve FD (sonlu farklar) yöntemi olarak isimlendirilir. Kendisi ve türevleri her noktada sürekli olan bir f(x) fonksiyonunun bir xo noktasındaki türevi, fonksiyonun o noktadaki eğimidir ve
Şekil 3.1. Sonlu farklar gösterimi. x xo f x xo f x dx df x f ∆ − ∆ + → ∆ = = ( ) ( ) 0 lim ) ( ' (2.1)
olarak tanımlanır. (2.1) bağıntısında Taylor açılımı kullanılır. Gerekli düzenlemeler yapılır ve limit operatörü de kaldırılırsa
) ( ) ( ) ( ) ( ' x Q x x f x x f x f + ∆ ∆ − ∆ + = (2.2)
şeklinde sayısal türevin tanımı kolaylıkla elde edilir. Đleri Sonlu Farklar (ĐSF) adı verilen (2.2) bağıntısındaki Q(∆x) ifadesi, katkılarının küçük olduğu varsayılan üst mertebeli terimleri göstermektedir. Sayısal türev tanımı, yine üst mertebeli terimlerin katkılarının küçük olduğu varsayımı ile
) ( ) ( ) ( ) ( ' x Q x x x f x f x f + ∆ ∆ ∆ − − = (2.3)
(2.3) bağıntısına Geriye Sonlu Farklar (GSF) yöntemi adı verilmektedir. Sayısal türevdeki hata miktarını azaltmak için Merkezi Farklar Denklemi (MFD) adı verilen üçüncü bir yöntem ise;
) ( 2 ) ( ) ( ) ( 2 ' x Q x x x f x x f x f + ∆ ∆ ∆ − − ∆ + = (2.4)
Şeklindedir MFD açılımı ise
) ( ) 2 1 ( ) 2 1 ( ) ( 2 ' x Q x x x f x x f x f + ∆ ∆ ∆ − − ∆ + = ( 2.5)
şeklinde de ifade edilebilmektedir.
MFD yöntemindeki hata (2.4) bağıntısında da görüldüğü gibi diğer iki yöntemden daha küçüktür. Türev açılımında daha fazla nokta kullanılarak hata istenilen mertebeye indirilebilir, fakat bu durumda bilgisayardaki hesaplama süresi ve bellek gereksinimi artacaktır. Maxwell denklemlerinin FD ile yazılması ve zamana göre türevlerin de sayısallaştırılarak genelleştirilmesi yöntemi FDTD adıyla özel olarak adlandırılmıştır. Böylece elde edilen ayrık denklemler ele alınan yapıya ait sınır koşulları da sağlanarak iteratif olarak çözülür[20,21].
4. FDTD YÖNTEMĐ
Đlk defa 1966 yılında Kano Yee[22] tarafından ortaya atılan bu yöntem, 1990’lı yıllarda bilgisayarların hız ve kapasitelerindeki büyük artışı sayesinde hemen her türlü probleme uygulanabilir hale gelmiştir. Genel halde üç elektrik alan ve üç manyetik alan bileşenlerinin uzayın seçilen ayrık noktalarında ayrık zaman aralıklarında hesaplanması ilkesine dayanır. Seçilen hesap uzayının her noktasında sürekli elektrik ve manyetik alanlar hesaplandığı için hızlı ve yüksek kapasiteli bilgisayarlara gerek duyulur.
Kaynak olarak sinüzoidal (dar bandlı) ya da darbesel (geniş bandlı) kaynaklar seçilebilir. Đteratif olarak her küçük zaman artımında söz konusu alan bileşenlerini, Maxwell denklemlerinden çözmeye dayanan FDTD yönteminde yakınsama, dispersiyon ve doğruluk problemleri çok önemlidir. FDTD yönteminin uygulanmasında en önemli ölçüt, çözülmeye çalışılan problemin fiziğinin iyi anlaşılması ve her adımda izlenmesidir.
Tek tonlu sinüzoidal işaretlerin simülasyonunu yapabildiği gibi, Gauss darbesi gibi darbesel işaretlerin de ele alınabilmesine olanak sağlar. Bu nedenle, bir simülasyonda çok geniş band da frekans yanıtları elde edilebilir[19].
FDTD yöntemi ile ilgilenilen yapılar yüksek doğrulukla ve kolaylıkla modellenebilmektedir. FDTD yöntemi, diğer bazı nümerik yöntemler gibi çok büyük matris yapıları içermemekte ve matris tersi alma gereksinimi olmadığı için ulaşılan çözümler oldukça güçlü olmaktadır. Zaman boyutunda yapılan çözümler sonucu elektromanyetik dalgaların görsel olarak istenilen zamanlarda hareketleri izlenebilmekte ve incelenen yapının dalga yayılımı açısından anlaşılması kolay olmaktadır[23].
FDTD metodunda problem uygun ızgara koordinat sistemine yerleştirilir. Ele alınan problemde üç boyutlu hesap uzayının, eş özellikli Nx x Ny x Nz adet dikdörtgen prizmadan oluştuğu varsayılır. Birbirinin özdeşi olan bu küçük prizmaların boyutları,
∆ x x ∆ y x ∆ z olup, FDTD hücresi olarak isimlendirilir. Hücre numaraları (i,j,k) sırasıyla, x'te, y'de ve z'deki hücre numarası olmak üzere tam sayılarla belirtilir. Maxwell denklemlerindeki diferansiyel operatörleri sonlu farklarla hesaplanır. Izgara düğümlerindeki alanlar ∆t zaman aralığı, n tamsayı olmak üzere ayrık zaman adımlarında (n∆t) bulunur. Bu işlem zamanda ilerleme olarak adlandırılır. Herhangi bir düğümde, t anındaki alan, bir önceki adımda düğümdeki ve komşu düğümlerdeki alanlardan hesaplanır[23].
Hesap uzayının ızgaralanması sırasında, ızgaranın boyutları modeldeki en küçük dalga boyundan daha küçük olmalıdır bu da çok büyük bir hesap uzayı ve çözüme ulaşmak için uzun zaman gerektirir. FDTD metodu ile elektrik ve manyetik alan bileşenleri hesap uzayının her noktasında doğrudan bulunur.
4.1. FDTD Yönteminin Kullanım Alanları
FDTD Yönteminin birçok uygulama alanı bulunmaktadır. Örnek olarak;
• Birçok ortamın modellenmesi için gerekli olan frekansa bağlı yapısal parametrelerin hesaplanması (kayıplı dielektrik, manyetik ortamlar),
• Radar simülasyonları,
• Sistemlerin geniş band cevaplarının incelenmesi, • Üç boyutlu geometrik yapıların incelenmesi,
• Herhangi bir iletkenliğe sahip nesnenin, düşük veya sıfır iletkenliğe sahip nesneyle etkileşiminin hesaplanması,
• Yakın alanlardan türetilen uzak alanları da hesaplayan (saçılmış alanlar, anten düzenleri, radar etkin kesiti, yüzey cevabı, akımlar, güç yoğunluğu gibi) konuların incelenmesi,
• Mikroşerit hatlı devre analizi,
• Elektromanyetik kaynak etkisinde kalmış insan kaslarındaki özgül emilim oranının hesaplanması(SAR),
• Açık ya da kapalı dalga kılavuzlarında dalga iletimi ve süreksizlikleri, • Mikrodalga fırın simülasyonu verilebilir.
4.2. Üç Boyutlu FDTD Yöntemi
Üç boyutlu (3-D) problemlerde uzaydaki ayrıklaştırma, şekil 4.1'de gösterildiği gibi Yee tarafından önerilen birim hücre kullanılarak gerçekleştirilmektedir.
Şekil 4.1 3-D birim YEE hücresi
3-D FDTD için Yee birim hücresi kullanıldığında elektrik ve manyetik alanların bileşenleri iteratif denklemlerle ele alınır. Kayıplı bir ortamda;
E E X t Hρ ρ ρ
σ
µ
=−∇ − ∂ ∂ 0 (3.1.a) H X t Eρ ρ ∇ = ∂ ∂ε
(3.1.b)şeklinde verilen ilk iki Maxwell denklemi (rotasyonel denklemleri) ayrıklaştırıldığında, n; zaman adımı ve (i,j,k) da, sırasıyla (x,y,z)'deki konum adımları olmak üzere alan bileşenleri
)] , 1 , ( ) , , ( [ ) 2 ( 2 )] 1 , , ( ) , , ( [ ) 2 ( 2 ) , , ( 2 2 ) , , ( 1 k j i H k j i H y t t k j i H k j i H z t t k j i E t t k j i E n z n z n y n y n x n x − − ∆ ∆ + ∆ + − − ∆ ∆ + ∆ − ∆ + ∆ − = −
σ
ε
σ
ε
σ
ε
σ
ε
(3.2)şeklinde bulunur. Burada n = n + l/2'dir. Hücre yapısı nedeniyle 3-D FDTD için, zamanın tam katlarında elektrik alanlar, kesirli katlarında ise manyetik alanlar hesaplanmaktadır. Burada verilen (3.1) denklemleri (µ=µo) boşluğun manyetik geçirgenlik katsayısı sahip, kayıplı dielektrik ortamlarda,
(
ε
=ε
0.ε
r)
geçerlidir[20].)] , 1 , ( ) , , ( [ ( )] 1 , , ( ) , , ( [ ) , , ( ) , , ( 1 k j i E k j i E y o t k j i E k j i E z o t k j i H k j i H n z n z n y n y n x n x − − ∆ ∆ + − − ∆ ∆ − = ≈− ≈
µ
µ
(3.3)Şekil 4.1 de verilen Yee Hücresi incelendiğinde aşağıdaki noktalar önemle göze çarpar. • Her birim Yee hücresinde üç elektrik ve üç manyetik alan bileşeni bulunur. Her
hücre (i, j, k) etiketi ile anılır. Bunlar sırasıyla x, y ve z'de ki hücre indisleridir. • Bir hücre içerisindeki altı bileşen de aynı (i,j,k) etiketi ile gösterilir şekil 4.1’den görüldüğü gibi, bu bileşenlerin hücre içi yerleşimi farklıdır. Örneğin, Ex(i,j,k) hücrenin x-kenar ortasında iken, Hz(i,j,k) hücrenin xy-yüzey ortasında bulunur. Yani elektrik alanlar hücre kenarlarında, manyetik alanlar hücre yüzeylerindedir.
• Aynı hücrede, elektrik ve manyetik alanların yerleşimleri gibi, hesaplandıkları zaman adımları da farklıdır. Elektrik ve manyetik alanlar birbirinden ∆t/2 kadar farklı zamanlarda hesaplanırlar. Yani elektrik alan bileşenleri t=0, t, 2∆t, 3∆t, vb. adımlarında hesaplanırken, manyetik alanlar t=∆t/2, 3∆t/2, 5∆t/2, vb. adımlarında hesaplanmaktadır. Böylece hesaplama bir elektrik alanlar bir manyetik alanlar diye iteratif olarak sürdürülür.
• Aynı hücrede belli bir noktada elektrik ve manyetik alanlarından söz edebilmek için konumda ve zamanda ortalama almak yeterlidir. Örneğin alan bileşenlerini hücre merkezine ötelemek için iki manyetik alan bileşeni
[
(, , ) ( 1, , )]
* 2 1 ) , , (i j k H i j k H i j k Hx = x + x + (3.4)yeterlidir. Ancak elektrik alan bileşenleri için komşu dört bileşene gerek vardır.
)] , 1 , 1 ( ) , 1 , ( ) , , 1 ( ) , , ( [ * 4 1 ) , , (i j k E i j k E i j k E i j k E i j k Ez = z + z + + z + + z + + (3.5)
Şekil 4.2 Tek birim FDTD hücresindeki ortalama elektrik alan
• Đteratif FDTD denklemlerinde herhangi bir ortam, üç ortam parametresi ile temsil edilir. Bunlar; dielektrik sabiti ε, manyetik geçirgenlik µ ısıl kayıpları temsil eden iletkenlik σ‘dır. Bunlardan ε ve σ elektrik alan bileşenlerinin hesaplandığı denklemlerde, µ ise manyetik alan bileşenlerinin hesaplandığı denklemlerde görünmektedir.
• Hücrelerde elektrik ve manyetik alan bileşenleri için farklı ortam parametreleri (ε, µ, σ) belirlenerek ince teller ya da değişik cisimler modellenebilmektedir.
• FDTD simülasyonu boyunca elektrik ve manyetik alanlar
güncellenmektedir. Buna karşın, bu alan değerlerinden herhangi bir noktada gerilim yada akım hesabı kolayca yapılabilir. Herhangi bir (i, j, k) hücresindeki gerilim ve akım Gauss ve Amper yasasından elde edilir. Örneğin Vz ve Iz
z k j i E k j i Vz(, , )=− z( , , )*∆ y k j i H k j i H x k j i H k j i H k j i Iz(, , )=[ x(, −1, )− x(, , )]*∆ +[ y( , , )− y( −1, , )]*∆ (3.6) kullanılarak hesaplanır.
• FDTD simülasyonunda hem sinüzoidal hem de darbesel kaynaklar modellenebilir. Örneğin, (i, j, k) hücresinde Vz (t) kaynağı için
yazmak yeterlidir. Burada genliği Ao frekansı f0 olan sinüzoidal kaynak için ) 2 sin( * ) * (n t A f n t Vz ∆ = o
π
o ∆ (3.8)M ve L zamanda öteleme ve darbe süresini ayarlayan tamsayılar olmak üzere Gauss tipi darbesel kaynak için ise
) / ) ( exp( * ) * ( 2 2 L M n A t n Vz ∆ = o − − (3.9) kullanılır.
4.3. Kararlılık (Courant) Koşulu
Çözümler iteratif olduğundan zaman ve konum aralıkları arasında belli bir kararlılık kriteri söz konusudur. Yani, konumda ayrıklaştırma yapıldıktan sonra zaman aralığı keyfi seçilemez.
FDTD tekniğinde zaman adımı ∆t 'nin konum adımları ∆x, ∆y ve ∆z 'den bağımsız değildir. Đteratif denklemlerin kararlı (sonsuza gitmeyen) sayısal çözümleri garantilemesi için hücre boyutları arasında bir ilişki sağlanmalıdır. Đşte Courant kriteri [1,2 ] ile konum ve zaman adımlan ışık hızı (c) ile birbirine
2 2 2 ( 1 ) ( 1 ) ) 1 ( 1 z y x c t ∆ + ∆ + ∆ ≤ ∆ (3.10)
şeklinde bağlıdır. Courant kriterini açıklayacak olursak; bir an için sadece zamana ve x 'e bağlı dalga problemi için ele alalım. Tek boyutlu FDTD de bu durumda (3.10) bağıntısı
n c x t . ∆ ≤ ∆ ; n=simülasyon boyutu (3.11)
0 0 1
ε
µ
= c ; (c=3.108m/sn) (3.12) z t c.∆ ≤∆ (3.13)şekline dönüşür. Benzer Şekilde 3-D FDTD içinde
2 1 2 2 2 ( 1 ) ( 1 ) ) 1 ( − ∆ + ∆ + ∆ ≤ ∆ z y x t c (3.14 )
olarak bilinmektedir [3]. Eşitliğin sol tarafı, c x ∆t (hızxzaman) yani dalganın aldığı yolu verir. Işık hızı dalganın maksimum hızı olduğuna göre dalganın alacağı maksimum yoldur. Sağ tarafında ise hücre boyutu vardır. O halde FDTD iteratif denklemlerinin kararlı olabilmesi için seçilen zaman adımında dalganın maksimum ilerlemesi hücre boyutunu aşmamalıdır. Başka bir deyişle, dalga hareketinin bir zaman adımında hücre içerisinde kalabilmesi için zaman adımı yeterince küçük seçilmelidir. Aksi takdirde dalgayı simule etmek imkansız olacaktır [19]. Burada Courent Kriterinin sadece kararlılığı sağladığı ancak doğruluğu garanti etmediğini belirtmek gerekir.
4.4. FDTD Yönteminde Sayısal Dağılma (Dispersiyon)
FDTD bağıntılarının sayısal simülasyonunda hücreler arası elektromanyetik dalga iletiminde dispersiyon olmaktadır. Yani, hücre içindeki elektromanyetik dalgaların faz hızları boşluktaki c hızından farklı değerler almaktadır. Bu farklılık frekansa, hücre boyutlarına ve dalga ilerleme yönüne göre değişmektedir. Sayısal dispersiyon, FDTD analizinin en önemli konularından birisidir. Özellikle geniş bandlı analizler yapılırken sayısal dispersiyonun etkileri doğruluğa, simülasyonun güvenilirliğine etkilerinin iyi anlaşılması gerekmektedir[20].
Sayısal dispersiyon Simüle edilen işaret içersinde en küçük dalga boyuna (en yüksek frekansa) sahip bileşenin konumda kaç hücre ile örnekleneceğine bağlı bir tanımdır.
Zaman-frekans ilişkisi içersinde biline Nyquist örnekleme teoremine çok benzeyen bir tanımdır. Nyquist teoremine göre bir işaretin bilgi kaybı olmadan tekrar elde edilebilmesi için zamanda örnekleme hızının içerdiği en yüksek frekansın iki katı olması gerekir. Benzer şekilde FDTD yönteminde sayısal dispersiyon minimum dalga boylu işaretin kaç konum örneği ile örnekleneceğine bağlıdır[20].
Şekil 4.3'te (a) şıkkında λmin = 2∆x (b) şıkkında ise λmin = 8∆x ile örneklenmiştir. Bu durumda (b) şıkkında FDTD'nin λmin bileşeninin dalga yayılımını daha iyi simüle edeceği açıktır.
(a) (b)
Şekil 4.3 Sayısal dispersiyonun fiziksel yorumu, (a) periyot başına iki örnek, (b) periyot başına sekiz
örnek
eğer zamanda ve konumda FDTD örneklemesi uygun boyutlarda yapılırsa, sayısal dispersiyonun etkisi istenilen dereceye indirilebilir[20].
4.5. Parametre Seçimi
FDTD yönteminde önemli unsurlardan birisi de uygun parametrelerin seçilmesidir. Yöntemde birbirini doğrudan etkileyen birçok parametre söz konusudur. Ele alınan yapıya, gerçeklenmek istenen simulasyona ve beklenen sonuçlara göre bu parametrelerin adım adım optimize edilmesi gerekir.
Birçok elektromanyetik problemde frekans domeni çözümleriyle ilgilenilir. FDTD ile elde edilen zaman domeni davranışından Fourier dönüşümü ile frekans davranışı elde edilir. Bu nedenle, özellikle ayrık Fourier dönüşümünde karşılaşılan örtüşme ve spektral sızıntı[24] gibi sorunlara karşı önlem almak zorunludur. FDTD
simulasyonunda darbesel kaynak kullanılarak geniş frekans bandlarında analiz yapılabilir. Đlgilenilen band içerisinde en yüksek frekans fmax ve frekans çözünürlüğü
∆f, zaman parametrelerinden T: gözlem süresi ve ∆t örnekleme aralığına bağlıdır. En yüksek frekans ile zamanda örnekleme aralığı birbirine
max 2 1 f tNyq = ∆ (3.15)
şeklinde bağlıdır. Görüldüğü gibi frekans analizi açısından parametreler keyfi seçilememektedir. Đlgilenilen frekans yükseldikçe (Nyquist örnekleme hızı) zaman domeninde işareti daha sık örneklemek gerekmektedir. Benzer şekilde
T f = 1
∆ (3.16)
olduğundan, daha iyi frekans çözünürlüğü için işaretin zaman domeninde daha uzun süre gözlenmesi (simulasyonu) gerekmektedir. FDTD için gerekli ∆t zaman adımı pratikte olarak seçilmektedir[25].
K t t≤ ∆ Nyq
∆ ; K>1 (3.17)
Zaman adımı bu şekilde belirlendikten sonra konum adımı için gereken alt ve üst sınırlar iki koşula göre hesaplanmaktadır. Courant kararlılık kriteriyle 3-D FDTD için konum adımı alt sınırı
z t
Vmax∆ ≤∆
3 (3.18)
şeklinde belirlenmektedir[6]. x ve y boyutları için de (3.18) denkleminde ∆z yerine sırasıyla ∆xve ∆y konulmaktadır.
Konum adımının üst sınırı, sayısal dispersiyon koşulu ile hesaplanmaktadır. Konum-işaret dalga boyu arasında da zaman-frekans ilişkisindeki Nyquist örnekleme
hızına benzer bir koşul bulunmaktadır. Konumda seçilen adım uzunluğuna bağlı olarak, işaret içindeki en küçük dalga boyu, yani en yüksek frekanslı bileşen, konumda en az iki düğüm ile gösterilmelidir. Üç boyutlu ve ∆z adımlarına sahip bir yapıda konumda en uzak iki nokta √3 ∆z olduğundan, üst sınır
2
3∆z≤λmin (3.19)
olarak hesaplanmaktadır [25]. Uygulamada bu değer √3 ∆z ≤ λmin/200 değerine kadar inebilmektedir.
5. FDTD YÖNTEMĐN HESAP UZAYINA UYGULANMASI
1- Önce uygulanacak kaynak, örneğin Gauss darbesi, seçilir. Kaynak darbe süresi analiz yapılacak en yüksek frekansa göre ayarlanır.
2-Gauss darbesi ve darbe süresi seçildikten sonra, ikinci adım zamanda ve konumda ayrıklaştırmadır.
3-Đteratif denklemlerinden de görüleceği gibi her hücrenin yapısı ε, µ ve σ parametreleri ile belirlenmektedir. Uygulama da her hücrede her elektrik alan bileşeni ayrı bir kimlik numarası ile anılır (i,j,k). Böylece bu alan bileşenlerinin olduğu noktalar ε ve σ ile modellenebilir. Kayıpsız ortamlarda (σ=0) , (ε= εo)’dur
3- Gauss darbesinin frekans bandındaki en büyük frekans bileşeni, Nyquist örnekleme teoremi gereği zamanda ayrıklaştırmayı belirler. Bu ayrık Fourier dönüşümü için gereklidir.
4-Yani, istenen en yüksek frekans bileşeni fmax ise işarette bilgi kaybı olmaması için
zaman adımı,
(
)
max 2 1 f tFFT = ∆ olmalıdır. Bu en üst sınırdır.5- FDTD tekniğinin uygulanmasında ise konumda hücre boyutları, ∆x, ∆y ve ∆z ile zaman adımı, ∆tFDTD yine ilgilenilen en yüksek frekanstan (fmax) başlayarak belirlenir.
6- FDTD açısından iki sınırlama söz konusudur. Birincisi sayısal dispersiyon, ikincisi kararlılık kriteridir. Yukarıda kısaca sözü edilen bu iki nokta hem FDTD simülasyonu hem de simülasyon sonrası ayrık Fourier analizi açısından önemlidir.
7- Ele alınan zaman işareti içersinde en yüksek frekanslı (fmax) yani en küçük dalga boylu (
λ
min) bileşen en az birkaç hücre ile örneklenebilmelidir. Yani FDTD hücre boyutları işaretin en küçük dalga boyuna sahip bileşeninden bile bir kaç misli küçük olmalıdır. Uygulamadaλ
min /100 ileλ
min/4 arasında hücre boyutlarınarastlanmaktadır. Genelde
λ
min /10-λ
min /20 arasındaki değerler uygun seçimlerdir.8- Bu şekilde hücre boyutları seçildikten sonra, kararlılık kriteri gereği zaman adımı,
FDTD
t
∆ belirlenir.
9-Genelde, simülasyon sonrası frekans analizi için gerekli ∆tFFT değeri simülasyon
adımı ∆tFDTD değerinden çok daha küçüktür ve FDTD zaman adımını belirlemek
diğerini otomatik olaraksağlar.
10- Ancak, simülasyon süresinin belirlenmesi frekans analizi açısından önemlidir. Zaman domeninde FDTD simülasyonu ne kadar sürecektir? Bu sorunun da iki açıdan yanıtlanması gerekir. Birincisi simülasyon sonrası frekans analizidir. Frekans analizinde hangi frekans çözünürlüğü (∆f) ile sonuçlar istenmektedir? Đstenen ∆f çözünürlüğü işaretin zaman domeninde hangi sürede (Tmax) gözleneceğini belirler. Bu ikisi birbiri ile ters orantılıdır. Yani, ∆f = 1 / Tmax olduğundan ne kadar hassas frekans ayrımı isteniyorsa o kadar uzun süre işaret gözlenmelidir (simülasyon sürdürülmelidir).
11- Tmax ve ∆tFDTD belirlendikten sonra simülasyon adımı NFDTD= Tmax/∆tFDTD olarak
seçilir.
12- Simülasyonda, ele alınan yapıda işaretin gözlenen zaman davranışı daha kısa sürede bitiyorsa FDTD simülasyonunu kesip işaretin öncesine veya sonrasına sıfır eklemeyle gerekli Tmax sağlanır. Tersine, Tmax süresi işaretin zaman davranışının sona ermesine yetmiyorsa simülasyon zaman davranışı sona erinceye dek sürdürülür. Örneğin rezonatör gibi çınlamaya neden olan yapıların incelenmesinde FDTD simülasyonunda, işaret zamanda çok yavaş söndüğünden, uzun süreler gerekir. Uygulamada bazen simülasyonu belli bir süre sonra kesip, işaret işleme teknikleri ile sonraki anlarda yapının davranışını kestirmek yoluna gidilir.
Özellikle FDTD hücre yapısının daha kolay anlaşılabilmesi ve zaman domeninde, örneğin darbe iletimi gibi, simülasyon sonuçlarının daha kolay gözlenebilmesi için tek yada iki boyutlu konum uzayındaki problemlerle başlamak daha uygun olabilir.
Kaynak noktasında zamanda darbe oluşmasının izlenebilmesi için tek ya da iki boyutlu yapılar daha basit olacaktır[19].
5.1. FDTD Algoritmasında Uyarma Fonksiyonu
FDTD yöntemi için zaman domeninde sonlu ve sınırlı frekans spektrumuna sahip dalga biçimi uygundur. Ancak, fiziksel olarak bu mümkün değildir, çünkü zaman domeninde sonlu bir işaretin frekans spektrumu sonsuzdur. Aynı şekilde frekans domeninde sınırlı bir işaret de zaman domeninde sonsuz olmaktadır.
Genel olarak işaretin zaman domenindeki şekli ile frekans band genişliği arasında ters ilişki vardır; işaretin süresi zaman domeninde kısaldıkça frekans spektrumu genişler. Buna göre, herhangi bir f(t) işaretinin etkin süresi ( Teff ) ile etkin band genişliği
(Beff) çarpımı bir sabite eşittir[25].
Teff x Beff = sabit. (4.1)
Zaman fonksiyonu f(t) ile bu fonksiyonun Fourier dönüşümü F(w) için etkin süre ve etkin band genişliği tanımları.
∫
∫
∞ ∞ − ∞ ∞ − − = dt t f dt t f t Teff 2 2 2 ) ( ) ( (4.2)∫
∫
∞ ∞ − ∞ ∞ − = dw w F dw w F w Beff 2 2 2 ) ( ) ( (4.3)şeklindedir. FDTD yöntemi için geniş bandlı kaynak modellemesinde en uygun işaret, zaman-band genişliği (Teff x Beff) çarpımı minimum olan işarettir. Bu tip bir
işaretin hem zaman süresi hem de frekans band genişliği oldukça küçüktür; bu sayede Fourier dönüşümündeki örtüşme hataları azalmaktadır. Teorik olarak zaman-band genişliği çarpımı
(TefJ x Beff ) ≥ ½ (4.4)
değerini almaktadır [25]. Bu koşulları en iyi sağlayan fonksiyon Gauss darbesidir. Gauss darbesi iki farklı şekilde tanımlanmaktadır
[
/(4 )]
exp ) ( 2 2 eff T t t f = − (4.5.a) ya da[
2 2]
exp ) (t B t f = − eff (4.5.b) (4.5) denklemlerinde tanımlanan Gauss darbesi sonsuz sürelidir. FDTD'de sonlu değerlerle çalışıldığı için bu Gauss darbesi belli bir To süreli ve birim genlikli bir işaretle (p(t/To)) çarpılarak sınırlandırılmaktadır. Bu şekilde, sonsuz süreli Gauss darbesi zamanda dikdörtgen bir darbe ile çarpılarak sonlu hale getirilmektedir. Bu fonksiyong(t) = f(t )x p(t/T0) (4.6)
şeklinde tanımlanmaktadır[25]. Gauss darbesini zamanda sınırlı yapmak Teff değerini azaltıp Beff 'i arttırmaktadır. Ancak uygun T0 değerleri yardımıyla bu
değişiklik çok küçük tutulabildiğinden, sınırlı Gauss darbesinin Teff ve Beff değerleri
orijinal Gauss darbesininkilerle eşit kabul edilebilir.
Sınırlandırılmış darbenin T0 değeri seçilirken istenen veya ulaşılabilir kesinlik
değerleri göz önüne alınmaktadır. emin, darbenin sıfırdan farklı sayılacağı değer olarak kabul edilirse, bu değerin altında gauss darbesinin sıfır olduğu söylenebilir ve Gauss darbesi T0'da
g(T0/2) = emin (4.7)
şeklinde sınırlandırılabilir[25]. Normalize Gauss darbesinin tepe değeriyle emin
arasındaki oranın logaritmik değerine dinamik sınır (dynamic range) adı verilmekte ve; ) ( log 20 ) 1 ( log 20 10 min min 10 e e Rdyn = =− [dB] (4.8)
şeklinde ifade edilmektedir. Sınırlandırılmış işaretin darbe süresi T0 istenilen
dinamik sınıra bağlı olarak hesaplanabilmektedir. Öncelikle, emîn dinamik sınır cinsinden 20 / 10 min dyn R e = − (4.9) şeklinde tanımlanmakta, buradan
[
2]
/20 0/2) 10 ( exp Rdyn effT B = − − (4.10)denklemi To için çözülerek
2 1 2 96 . 0 ) ( log 20 2 10 0 dyn eff dyn eff dyn eff R B R B e R B T = ≈ ≈ ( 4.11) eşitliği elde edilmektedir [25]. Bu denklem, yaklaşık olarak
dyn eff R T To 2 ≈ (4.12)
şeklinde yazılabilir. (Rdyn [dB] cinsindedir.) (4.6) denklemi yardımıyla herhangi bir
dinamik sınır için emin = 10‾6 elde edilmektedir; yani Gauss darbesinin genliği 10‾6 değerinden küçük olduğunda sıfır olarak kabul edilirse dinamik olarak 120 [dB] bulunmaktadır. Bu değer (4.6) denkleminde kullanıldığında
eff T T0 ≥16 Veya eff B T0 ≥ 8 (4.13) elde edilmektedir [25].
Sınırlandırılmış Gauss darbesi normalde zamanda T0/2 kadar ötelenmiş olarak
yazılmaktadır. Bu sayede darbenin t < 0 anlarındaki tüm değerleri sıfıra eşit olacağı için nedensellik koşulu sağlanmış olmaktadır. Bu durumda sınırlandırılıp ötelenen Gauss darbesi
[
]
[
0 0]
2 0 2 ( /2) ( /2)/ exp ) (t B t T p t T T g = − eff − − (4.14) olarak yazılmaktadır. 120 [dB]'den büyük Rdyn için To = 8/Beff seçilirse[
Befft t T] [
p t Beff Beff]
t
g( )=exp−( −4)2( − 0/2)2 ( −4/ )/(8/ (4.15)
elde edilmektedir.
5.2. Gauss Darbesinin Ortama Enjekte Edilmesi
Elektromanyetik dalganın düzlem içerisinde uyarılması birkaç değişik şekilde gerçekleştirilmedir. Normali ve en basiti düzlem dalganın konumda belli bir anda yerleştirilmesi ve simülasyonun başlatılmasıdır. Ancak düzlem dalga konumda düzlem karakteri taşırken zamanda sinüzoidal ya da darbesel olabilmektedir. Zamanda bir Gauss darbesi[23];
− − = 2 0 exp ) ( T t t t g (4.16)
Şeklinde ifade edilir. Burada t darbenin simülasyon başlangıcından ne kadar sonra 0
enjekte edileceğini T ise darbenin genişliğini belirleyen parametrelerdir. Zamanda darbe ne kadar dar (yani T ne kadar küçük) ise işaretin frekans bandı ( zaman band genişliği sabit olduğundan ) o kadar geniş olacaktır. Benzer şekilde konumda bir Gauss darbesi ise[23]; − − = 2 exp ) ( g s z z z z g (4.17)
Şeklindedir. Burada Zs darbenin uygulanacağı konum noktasını ve Zg ise konumda darbe genişliğini belirleyen parametrelerdir. Konumda darbe ne kadar dar ( yani Zg ne kadar küçük) ise işaretin dalga sayısı zaman domeninde kapsadığı alan o kadar geniş olur. Zaman simülasyonu sırasında düzlem dalga adım enjekte edilmiş olup yani[23],
) ( ) ( ) (k E k g n Ekn s = kn s + (4.18) eşitliği kullanılmıştır. Burada Ks konumda Gauss darbesinin uygulanacağı noktayı göstermekte ve K=1 ile K=Ke arasında seçilebilmektedir ( Ke maksimum nokta olmak üzere) zamanda ayrık gauss fonksiyonu g(n) ise ( to=n0 x ∆t ve T =nTx∆t) olma
üzere[23] − − = 2 0 exp ) ( T n n n t g (4.19)
şeklindedir. Gauss darbesini kaynak olarak kullanmanın bir diğer yolu da Eşitlik 4.17’de ifade edildiği gibi doğrudan konumda uyarmak ve simülasyonu başlatmaktır. Bu durumda zaman simülasyonu başlamadan önce konumda düzlem dalganın oluşturulması gerekir. Bu amaçla ( Zs=ks x∆z ve Zg=kg x ∆z olmak üzere)
− − = 2 0( ) exp g s k k k k k E (4.20) 0 0( ) ( ) η k E k H o x y = (4.21) Ω = = 0 0 377 0 ε µ η kullanılır.
Darbesel işaretlerin FDTD uzayında simülasyonu için önemeli sorunlardan bir tanesi gauss darbesi için darbe süresi ve gecikme parametrelerinin ayrık ∆z ve ∆t cinsinden belirtilmesidir. Genelde istenen darbenin konumda toplam uzayın ancak belli bir oranını kapsaması, bu sayede dalganın zamanda yayılmasının net olarak izlenmesidir. Bu amaçla değişik yollar izlenilir:
1) Operasyonel parametreler cinsinden ∆z verildiğinde ve kararlılık koşulundan ∆t belirtildiğinde darbe süresi doğudan ∆z cinsinden belirlenebilir. Bunun için en yüksek frekanslı işaret bileşeninin dalga boyu ∆z cinsinden belirlenebilir. Örneğin B=ƒmax=c/(10x∆z) seçilebilir. Sonrasında iki parametre α=3.3β² ve
α ϕ= 4
Şeklinde seçilerek gauss darbesi g(t)=exp
{
−α(
1−ϕ)
2}
olarak uygulanır. Bu sayedekonumda λmin/10 örnekleme ile sayısal dispersiyon koşulu sağlanmış olur.
2) Đkinci yol Gauss darbe parametrelerinin en yüksek frekansa (ya da bant genişliğine) göre belirlenmesi ve ayrıklaştırma parametrelerinin yani (∆z ve ∆t’nin) verilen banda göre belirlenmesidir. Bunun için B veya ƒmax verildiğinde ∆z=c/(10xB) seçilerek
ayrıklaştırma gerçeklenir ve diğer iki parametre yine α=3.3β² ve ϕ = 4 α şeklinde seçilerek Gauss darbesi birinci şıktaki gibi g(t)=exp
{
−α(
1−ϕ)
2}
uygulanır.3) Bir diğer yol ise α =16
(
ϕx∆t)
2 ve ϕ =32 (ϕ =16 veya ϕ =64) seçerek Gaussdarbesini ∆t cinsinden g(t)=exp
{
−α(
1−ϕx∆t)
2}
olarak belirlenmektedir[23]. Buuygulamada ilk yöntemdeki yani; B=ƒmax=c/(10x∆z), α=3.3β² ve α ϕ= 4 ve
(
)
{
2}
1 exp ) (t = −α −ϕg Şeklindeki Gauss fonksiyonu kullanılmıştır.
5.3. FDTD Algoritması
FDTD algoritmasında döngü zaman döngüsüdür. Seçilen maksimum zaman adımı tamamlanıncaya kadar ana zaman döngüsü çalıştırılır. Ayrı ayrı işleri yapan alt programlar şeklinde düzenlenmesi, FDTD algoritmasının gerek işleyişini, gerekse kullanımı ve anlaşılmasını kolaylaştırır.
Şekil 5.1'den de görüleceği gibi, ana FDTD algoritması önce ilgili yapıyı yerleştirir, gerekli parametrelerin hesabını yapar yani program boyunca kullanılacak sabit parametreler girilir. FDTD parametre seçim kısmında anlatılan koşullar sağlandıktan sonra, yani kararlılık, sayısal dispersiyon v.b koşullarda tanımlandıktan sonra zaman döngüsünü başlatır. Ana FDTD döngüsü t=0 anında başlar. Önce elektrik alan bileşenleri bütün hücrelerde hesaplanır. Bitiminde sınır yüzeylerde ışınım koşullarının (ABS;açık bölge simülasyonu) sağlanması gerçeklendikten sonra manyetik alanların hesabına geçilir. Ancak, elektrik ve manyetik alanların hesaplandıkları anlar birbirinden ∆t/2 kadar farklı olduğundan önce zaman adımı ∆t/2 kadar artırılır. Bunların bitiminden tekrar zaman ∆t/2 kadar artırılır. Böylece simülasyon iteratif olarak sürdürülür. Yani, bir elektrik alanlar, bir manyetik alanlar hesaplanarak işleme devam edilir. Arada sınır yüzeylerde ABS gerçeklenir. Programda simüle edilen kaynak istenen çıkışların yazdırılması, her türlü gerilim, akım ya da güç hesaplarının yapılması ve (isteniyorsa) uzak alan dönüşümlerinin hesaplanması bu döngü içerisinde yerini alır[26].
5.4. Sınır Koşulları
Hiç bir bilgisayarın, sonsuz bir uzay için gerekli hesaplamayı ve veri depolamasını yapabilmesi mümkün değildir. Bu nedenle FDTD’deki hesap uzayının da sınırlandırılması gereklidir. FDTD metodunun en büyük sorunlarından biri EM dalga etkileşim probleminin sınırlı olmayan bölgelerde etkin ve doğru çözümüdür[27]. Bir metodun serbest uzayı simule edebilmesi için hesap uzayının sınırsız olduğunu da simule edebilmesi gerekmektedir[28]. Bu tip problemlerde sonsuza uzanan döngüyü simule edebilmek için dış döngü sınırlarında Emici Sınır Kosulları (Absorbing Boundary Condition, ABC) tanımlanmalıdır[27]. Emici ya da açık sınır koşulları mevcut olmadığında ızgaranın dış sınırındaki dalgalar hesap uzayının içine geri yansırlar[29]. Dolayısı ile emici sınır koşulu FDTD hesaplamasının doğruluğunu doğrudan etkilemektedir[30].
FDTD denklemleri incelendiğinde, ilk hücre için hesaplamada kullanılacak olan aynı noktadaki bir önceki alan değerini barındıran hücrenin mevcut olmadığı, son hücre için ise yine hesaplamada kullanılması gereken bir sonraki alan değerlerini içeren hücrenin de mevcut olmadığı görülecektir. Bu hücreler mevcut olmadığı için bu değerler sıfır gibi algılanacak ve tam yansıma söz konusu olacaktır. Dolayısı ile FDTD hesap uzayının sınır yüzeylerinde iteratif eşitlikleri hiçbir önlem almaksızın aynen uygulamak tam yansımalara neden olacaktır. Bu ise ele alınan yapının değişmesi demektir[19]. Bu yansıma sorununu engellemek amacı ile kullanılan emici sınır koşulları arasında en verimli kullanılabileni Mükemmel Uyumlu Katmanlar (Perfectly Matched Layer, PML) yöntemidir. Kendisinden önce kullanılan yöntemlere kıyasla büyüklük olarak yaklaşık 3-4 kat daha az hata vermektedir[31]. Jean Pierre Berenger tarafından 1994’te önerilen PML yöntemi, gelen dalgaların hesap uzayının dış sınırına konumlandırılan metal plakalar tarafından yansıma olmaksızın tamamen emilmesi temeline dayanmaktadır. Burada dikkat edilmesi gereken PML katmanının sadece teorik olarak mükemmel yansıtmasız olduğudur. Gerçek FDTD uygulamalarında PML materyalinin kalınlığına ve iletkenlik profiline bağlı olarak bu katmanlardan bazı yansımalar söz konusudur[31]. Bu olayın fiziksel uygulama olarak bir benzeri yalıtılmış
(anechoic)odaların duvarlarıdır[27]. PML’de performansı geliştirmek amacı ile Tek Eksenli Mükemmel Uyumlu Katmanlar (Uniaxial Perfectly Matched Layer, UPML)[32] ve Kıvrımlı PML (Convolutional Perfectly Matched Layer, CPML)[33] gibi farklı modifikasyonlar önerilmektedir.
6. ELEKTROMANYETĐK PROBLEMLERĐN FDTD YÖNTEMĐYLE ANALĐZĐ
Bu kısımda şu ana kadar teorisi verilmiş olan FDTD ile ilgili uygulamalara geçilecektir. Özellikle zaman domeninde düzlem dalga ya da küresel dalga ayrımının gözlenebilmesi için tek boyutlu simülasyonlardan başlanıp 2-D ve 3-D simülasyonlara yer verilecektir.
6.1. Bir Boyutlu simülasyon
a) Serbest uzayda maxvell in curl denklemleri
∂E ∂t = 1
ε
0 ∇ × H (6.1 a) ∂H ∂t = − 1µ
0 ∇ × E (6.1 b)olarak yazılır. Burada E ve H vektörleri Üç boyutludur εo=dielektrik sabiti (uzayın) =8.8542x10−12F /m
µo=manyetik geçirgenlik (uzayın)=4πx10−7H/m
Eşitlikleri vektör formunda yazarsak, Del burada Türevsel operatördür.
z k y j x i ∂ ∂ + ∂ ∂ + ∂ ∂ = ∇ (6.2) z y x E E E z y x k j i xE ∂ ∂ ∂ ∂ ∂ ∂ = ∇ (6.3)
Benzer şekilde elektrik alan için uyguladığımız işlemi manyetik(H) alan için yazarsak
∂ ∂ − ∂ ∂ ∂ ∂ − ∂ ∂ ∂ ∂ − ∂ ∂ = ∂ ∂ ∂ ∂ ∂ ∂ = ∂ ∂ ∂ ∂ ∂ ∂ y H x H x H z H z H y H H H H z y x k j i t E t E t E x y z z y z z y x z y x 0 0 1 1
ε
ε
(6.4.a) ∂ ∂ − ∂ ∂ ∂ ∂ − ∂ ∂ ∂ ∂ − ∂ ∂ − = ∂ ∂ ∂ ∂ ∂ ∂ − = ∂ ∂ ∂ ∂ ∂ ∂ y E x E x E z E z E y E E E E z y x k j i t H t H t H x y z z y z z y x z y x 0 0 1 1µ
µ
(6.4.b)Bu denklemleri daha da basitleştirerek 1 boyuta indirgersek elektrik alan enine salınım yapar ve z yönünde ilerler o zaman Ey=0’dır. Ve biz biliyoruz ki elektik alan ile manyetik alan birbirine dik olarak hareket eder onların nokta çarpımı sıfır olmalıdır. 0 ) 0 , , ).( 0 , 0 , ( .H = Ex Hx Hy =EzHz = E (6.5) ∂ ∂ = ∂ ∂ 0 0 1 0 0 0 z H t Ez y ε (6.6.a) ∂ ∂ − = ∂ ∂ 0 0 1 0 0 0 t E t H z y
µ
(6.6.a) Böylece H=(0,Hy,0) (6.4) eşitlikleri şu hale dönüşür.Yada eşdeğer olarak, bir boyutlu simülasyon için aşağıdaki ifade yazılabilir. ∂Ex ∂t = − 1 ε0 ∂Hy ∂z (6.7.a) ∂Hy ∂t = − 1 µ0 ∂Ex ∂z (6.7.b) Bu durum aşağıdaki gibi görselleştirilebilir.
Şekil 6.1 Ex ve Hy’nin uzayda birbirine göre konumları
Geçici ve uzaysal türevlerin her ikisi ile beraber merkezi farklar yaklaşımı kullanılarak aşağıdaki denklemler elde edilir.
Exn+1/ 2 (k)− Exn−1/ 2 (k ) ∆t = − 1 ε0 Hyn( k+1 / 2) − Hyn( k −1 / 2) ∆x (6.8.a) Hyn+1( k+1 / 2) − Hyn( k+ 1/ 2) ∆t = − 1 µ0 Ex n+1/ 2( k+ 1) − Ex n+1/ 2( k) ∆x (6.8.b)
Burada n zaman belirtir ve t = ∆t.n şeklinde ifade edilir. Simülasyonda gerekli kod yazılımı için her şeyin ayrıştırılması zorunludur n+1 terimi bir sonraki zaman adımını ifade eder, k mesafeyi belirtir ve kullanımının daha fazla duyarlılık sağladığı görülebilir. Bununla birlikte uzaysal artışlarda genel olarak ∆x kullanılır ve burada da ∆x kullanılmıştır. (6.8.a) ve (6.8.b) denklemleri E ve H alanlarının mesafeyle ve zamana bağlı olarak birbirinden ayrıldığını farz eder. E alan değerleri arasına yerleştirilmiş
olarak kabul edilen H alan değerlerini göstermek için k+1/2 ve k–1/2 ifadeleri kullanılır. Bu şekil 6.2’de gösterilmektedir. Benzer olarak n+1/2 ve n–1/2 ifadeleri de sırasıyla n 'den önceki ve n 'den sonraki değerleri göstermektedir. (6.8.a) ve (6.8.b) iterasyonla oluşturulan bir algoritmayla yeniden düzenlenebilir[34].
Ex n+1/ 2 (k)= Ex n−1/2 (k )− ∆t ε0⋅ ∆x Hy n (k+1 / 2) − Hy n (k −1 / 2)
[
]
(6.9.a) Hy n+1 (k + 1 / 2) = Hy n (k+ 1/ 2) − ∆t µ0⋅ ∆x Ex n+1/ 2 (k+1) − Ex n+1/ 2 (k)[
]
(6.9.a)Şekil 6.3 Hesaplama yapılan hücrelerde elektrik ve manyetik alan gösterimi
Dikkat edilmelidir ki hesaplamalar zaman ve mesafe olarak birbirinden ayrılmıştır. Örneğin 6.9.a denkleminde Ex'in yeni değeri, önceki Ex değerinden ve Hy 'nin son değerlerinden hesaplanır. Bu, zamanda sonlu farklar metodunun temel paradigmasıdır. 6.9.a ve 6.9.b denklemleri birbirine çok benzerdir fakat buradaki fark εo ve µo’ın genlik değerinde oluşturduğu farktır, dolayısıyla Ex ve Hy genlik yönünden farklı olacaktır. Bundan aşağıdaki değişiklik yapılarak kaçınılır.
˜
E =
ε
0µ
0E (6.10)
Burada E% normalize elektrik alan değeridir.
˜ E x n+1/ 2 (k)= ˜ E x n−1/2 (k )− 1
ε
0µ
0 ∆t ∆x Hy n (k+ 1/ 2) − Hy n (k−1 / 2)[
]
(6.11.a)Hyn+1(k + 1 / 2) = Hyn(k+ 1/ 2) − 1
ε
0µ
0 ∆t ∆x E ˜ x n+1/ 2 (k+ 1) − ˜ E xn+1/ 2(k)[
]
(6.11.b)Böylece Maxwell denklemlerinden ikisi bir boyuta indirgenip sonlu farklar uygulanarak son şeklini alır. Bundan sonraki kısımda Önceden belirtilen kurallara uyularak simülasyonlar başlatılabilir. Hücre boyutu ∆x seçildiğinde, zaman adımı ∆t hesaplanabilir. ∆t = ∆x 2⋅ c0; Burada c0 ışığın hızıdır. (6.12) Mesafe x x c x c t c ∆ = ∆ = ∆ = ∆ = η 2 1 2 0 0 0 (6.13) 1
ε
0µ
0 ∆t ∆x = c0⋅ ∆x / 2 ⋅ c0 ∆x = 1 2 (6.14)Bu durum kararlılık kısmında açıklanmıştı. Son olarak denklemler bilgisayar koduna aşağıdaki gibi çevrilir.
ex[k] = ex[k] + 0.5*( hy[k-1] - hy[k] ) (6.15.a)
hy[k] = hy[k] + 0.5*( ex[k] - ex[k+1] ) (6.15.b)
denklemlerde; n, (n+1/2) ve (n–1/2) ifadeleri bulunmadığına dikkat edilmelidir. FDTD metotta zaman, ima edilir. (6.15.a) denkleminde eşitliğin sağ tarafındaki Ex değeri, n–1/2 de ki önceki Ex değerini, sol taraftaki Ex değeri ise n+1/2 anındaki son Ex değerini göstermektedir. Bununla beraber pozisyonlar açıkça belirtilmiştir. Tek fark program içinde gösterilimde (k+1/2) ve (k–1/2) yerine k ve (k–1) kullanılmasıdır.
0 50 100 150 200 250 -0.5 0 0.5 1 Hücre Sayısı E x
Elektrik pulsenin simülasyonu
0 50 100 150 200 250 -1 -0.5 0 0.5 1 Hücre Sayısı H y
Manyetik pulsenin simülasyonu
Şekil 6.4Bir boyutta elektrik ve manyetik alanın konumdaki hareketleri
Program; problem uzayının merkezinde bir Gauss darbesi üretir ve darbe şekil 6.4’de gösterildiği gibi her iki yönde yayılır. Her iki yönde de Ex alan pozitiftir, fakat Hy alanı negatif yönde negatif değerler alır. Hesap uzayının bittiği sınır noktalarında herhangi bir değer olmadığı için eğer önlem alınmazsa program tarafından sıfır olarak algılanacaktır, bu da sınır noktasında tam yansıma meydana gelmesine sebep olacaktır.
6.2. Bir Boyutta Soğurulma Sınır Şartları
Soğurulma sınır şartları problem uzayına yansıyan E ve H alanlarının hesabı için gereklidir. Normalde E alanlarının hesabında etrafındaki tüm H alanların değerlerinin bilinmesine ihtiyaç duyarız. Bu FDTD metodunda temel kabuldür. Bununla birlikte problem uzayı dışında hiçbir kaynak olmadığını bilmek gibi bir avantajımız bulunmaktadır. Bu yüzden sınır noktasına ulaşan alanlar dışarıya doğru yayılır sınırdaki değerleri bir sonraki değerlere bağlı olarak hesaplamak için bu iki durum kullanılacaktır.
K=0 durumunu incelediğimizi farz edelim, bir dalga serbest uzayda sınıra doğru ilerliyorsa bu dalganın hızı
0
c ’dır. Yani FDTD algoritmasında bir zaman adımında
Mesafe= 0 c x ∆t = 0 c .∆x/2 0
c = ∆x/2 kadar ilerler. Bu denklem bir dalganın bir hücreyi
geçmesinin iki zaman adımı alacağını gösterir yani sınır şartlarının ( )= n−2(1)
x n
x o E
E
olabileceği ortak bir yaklaşım olarak kabul edilebilir. Bunu geliştirmek çok kolaydır. Temel olarak 2 zaman adımı için Ex(1) de bir değer depolanır ve daha sonra bu
) 0 (
x
E da yerine konulur. Bu haliyle merkeze yerleştirilen darbe dışarıya doğru yayılır
ve problem uzayında geriye yansıma olmadan emilir[34].
0 20 40 60 80 100 120 140 160 180 200 -1 -0.5 0 0.5 1
Elektrik pulsenin simülasyonu
E x x hücreleri 0 20 40 60 80 100 120 140 160 180 200 -1 -0.5 0 0.5 1
Manyetik pulsenin simülasyonu
x hücreleri
H
y
Şekil 6.5 Sınır şartları uygulandıktan sonra elektrik ve manyetik alanın durumu
6.3. Dalganın Dielektrik Bir Ortama Çarpması
Dielektrik sabiti birden farklı olan, yani boşluk olmayan bir ortamdaki yayılımın
xH
t
E
r∇
=
∂
∂
ε
ε
01
(6.16.a)xE
t
H
∇
−
=
∂
∂
01
µ
(6.16.b)simülasyonunu gerçekleştirebilmek için öncelikle bağıl dielektrik sabiti
ε
rMaxwell denklemlerine ilave edilir. Bir boyutlu simülasyon için yeniden normalize edilerek yazılırsa; z t H t t E Y r x ∂ ∂ − = ∂ ∂~ ( ) 1 ( ) 0 0µ ε ε (6.17.a)z
t
E
t
t
H
Y x∂
∂
−
=
∂
∂
(
)
1
~
(
)
0 0µ
ε
(6.17.b)Çıkan son denklem sonlu farklar yaklaşımı kullanılırsa
x k H k H t k E k E yn n y r n x n x ∆ − − + = ∆ − − + ( ) ~ ( ) 1 ( 1/2) ( 1/2) ~ 0 0 2 / 1 2 / 1 µ ε ε (6.18.a) x k E k E t k H k H n x n x n y n y ∆ − + − = ∆ + − + + + + ) ( ~ ) 1 ( ~ 1 ) 2 / 1 ( ) 2 / 1 ( 1/2 1/2 0 1 µ (6.18.b)
sonlu farklara uygulanmış hali elde edilir denklem düzenlenirse
2 1 1 0 0 = ∆ ∆ x t µ ε
olduğuna göre bunu 6.18 denklemine uygularsak
)] 2 / 1 ( ) 2 / 1 ( [ 2 / 1 ) ( ~ ) ( ~ +1/2 = −1/2 + + − − k H k H k E k E ny yn r n x n x ε (6.19.a) )] ( ~ ) 1 ( ~ [ 1 ) 2 / 1 ( ) 2 / 1 ( 1/2 1/2 0 1 k E k E k H k Hny+ + = ny + − xn+ + − xn+ µ (6.19.b)
son olarak (6.19.a) ve (6.19.b) elde edilir dielektrik ortamın boşluktan farkının sadece
