Розробка методики для отримання числових рішень диференціальних рівнянь

Розробка методики для отримання числових рішень диференціальних рівнянь

1.     
Огляд існуючих програмних комплексів

OptiFDTD представляє собою потужне, високо інтегроване і зручне середовище CAD, яке дозволяє проектувати і моделювати передові пасивні і нелінійні фотонні компоненти.дозволяє розробляти, аналізувати і тестувати компоненти для поширення хвиль, розсіювання, відбиття, дифракції і поляризації. В основі програми OptiFDTD - метод кінцевих різниць у часовій області (FDTD) алгоритм другого порядку точності і самі передові граничні умови - одноосьовий ідеально узгоджений шар (UPML).

Алгоритм вирішує електричні і магнітні поля в тимчасових і просторових доменах. Це дозволяє для будь-яку геометрію моделі і не накладає ніяких обмежень на властивості матеріалу.

Застосування:

·        поверхневого плазмового резонансу (SPR);

·        нано-частинки і клітини тканин;

·        дифракційні мікро-оптичні елементи і лінзи;

·        комплекс інтегрованих структур оптики;

·        нелінійні матеріали, дисперсійні матеріали;

·        оптичні мікро-кільцеві фільтри і резонатори;

·        решітки на основі хвилеводних структур;

·        електромагнітні явища;

FDTDpro дозволяє спостерігати в процесі рахунку картини поширення електромагнітних хвиль, обчислювати струми, напруги, хвильовий опір, резонансні частоти, поглинену енергію, питому поглинену потужність (SAR), потужність випромінювання антен.

Програма FDTDpro являє собою один exe-файл, в який вбудовані:

·        обчислювальний модуль;

·        тривимірний редактор об'єкта;

·        імпорт об'єктів з формату dxf;

·        модуль візуалізації із записом на диск послідовної серії картинок поширення електромагнітних хвиль;

·        ряд інших сервісних функцій.

Програму FDTDpro можна використовувати:

·        у навчальному процесі при вивченні теорії електромагнітних хвиль;

·        в інженерній практиці для розрахунків різних пристроїв (хвилеводів, СВЧ-резонаторів, деяких типів антен);

·        для проведення розрахунків з електромагнітної сумісності;

·        для вивчення впливу синусоїдальних та імпульсних електромагнітних полів на живі організми, в т.ч. на людину (наприклад, СВЧ-випромінювання, випромінювання мобільних телефонів);

·        для розрахунків ефективної площі відбиття, полів в далекій зоні, діаграм спрямованості.

Програма має масу налаштувань та змогу ознайомитись з можливостями в дещо обмеженому режимі.

FDTD Solutions це 3D-застосунок, здатний аналізувати взаємодії ультрафіолетового, видимого та ІЧ-випромінювання зі складними структурами функції довжин хвиль. FDTD Solutions дає можливість точно враховувати матеріальну дисперсію в широкому діапазоні довжин хвиль за допомогою власного Multi-коефіцієнта моделювання.

Основні переваги FDTD Solution

·        зниження витрат на розробку продукту за допомогою високоточних алгоритмів з вбудованою оптимізацією, яка дозволяє швидке віртуальне прототипування, щоб уникнути дорогого фізичного прототипування.

·        скорочення часу виходу на ринок завдяки високо оптимізованого ядра моделювання, розробленого для високої пропускної оцінки передових обчислювальних систем.

·        підвищення продуктивності за допомогою інструментів проектування, використання на з метою сприяння швидкому навчанню і швидкого початку роботи.

Особливості FDTD Solution

·        моделювання довільної геометрії в 2D, 3D

·        паралельні обчислення на декількох комп'ютерах

·        дисперсійні, нелінійні, посилюючі і анізотропні матеріальні можливості моделювання

·        паралельні обчислення на багатоядерних і багатовузлових системах

·        оптимізація обчислювального ядра

·        розширена сітка алгоритмів

·        потужна скриптова мова

·        налаштування дизайну та ієрархічний макет

·        фільми Моделювання динаміки.

2.     
Теоретичні відомості

2.1 Метод FDTD

Метод кінцевих різниць у часовій області (англ. Finite Difference Time Domain, FDTD) - один з найбільш популярних методів чисельної електродинаміки, заснований на дискретизації рівнянь Максвелла, записаних в диференціальній формі.

У 1966 р. Йі (Yee) розробив техніку, що реалізує явну звичайно - різницеву схему другого порядку для вирішення вихрових рівнянь Максвелла у просторі та часі.

Вихідними є рівняння Максвелла в диференціальній формі:

(H) = ∂D/∂t + J    (E) = - ∂B/∂t                                              (2.1)

а також:

D = ε ε0 E

J = σ E                                                       (2.2)

B = μ μo H

Тут E - вектор напруженості електричного поля (В / м), Н - вектор напруженості магнітного поля, (А/м), ε, μ - відносні діелектрична і магнітна проникності (без розмірності), ε₀ - діелектрична постійна (Ф / м), μ₀ - магнітна постійна (Гн / м), B - вектор магнітної індукції (Тл), D - вектор електричного зміщення (Кл/м²), J - вектор густини струму (А/м²), σ - електрична провідність (См / м), і t - час у секундах.

ε₀ = 10⁷/(4πc²),

де с - швидкість світла в вакуумі (2,997925010∙10⁸ м/с).

μ₀ = 4π/10⁷.

У рівняннях Максвелла зміна електричного поля E (приватна похідна) залежить від розподілу в просторі магнітного поля H (ротор). Аналогічно, зміна поля H залежить від розподілу в просторі поля Е.

Е = Ex (t, x, y, z) X+ Ey (t, x, y, z) Y+ Ez (t, x, y, z) Z

H = Hx (t, x, y, z) X+ Hy (t, x, y, z) Y+ Hz (t, x, y, z) Z                (2.3)

На цьому спостереженні заснований алгоритм Йі. Сітки для полів E і H зміщені по відношенню один до одного на половину кроку дискретизації часу і по кожній із просторових змінних. Кінцево-різницеві рівняння дозволяють визначити поля E і H на даному часовому кроці на підставі відомих значень полів на попередньому.

Всі компоненти (Ex, Ey, Ez, Hx, Hy, Hz) знаходяться в різних місцях, тобто рознесені в просторі. Е - компоненти знаходяться посередині ребер, Н - компоненти - по центру граней. Всі компоненти незалежні один від одного, тобто кожній з них можна присвоїти свої унікальні електричні (для Е) і магнітні (для Н) параметри.

Просторові координати кожного вектора x, y і z виражаються в номерах осередків i, j і k відповідно, час t виражається в кроках n за часом:

= i∆x= j∆y= k∆z                                                (2.4)= n∆t  

Поля E і H обчислюються із зрушенням на півкроку за часом. Позначення, введені Yee, наступні: En - значення поля E на тільки що обчислення кроці; Eⁿ⁺¹ - значення поля E на обчислюваному зараз кроці за часом. H^n-1/2 - значення поля H на тільки що обчислення кроці; H^n+1/2 - значення поля на обчислюваному зараз півкроку за часом. З цих позначень випливає, що процедура обчислень починається з поля H^n+1/2, тому що в момент t = 0 (n = 0) встановлені початкові умови по всьому счетному обсягом: всі значення полів E і H дорівнюють нулю. Хоча в принципі це лише найбільш поширена умовність. Можна вважати, що просторова сітка проходить через вектор H, що процедура рахунку починається з поля E.

Тепер, коли введені основні позначення, покажемо висновок виразів, придатних для розрахунків за допомогою комп'ютера.

Поставимо (2.3) і (2.2) в (2.1). отримаємо:

(H) X = εε₀∂Ex /∂t + σEx; rot(E) Y = - μμ₀∂Hy /∂t       (2.5)

Застосовуючи кінцево-різницеву апроксимацію, перетворимо (5) у вирази для кроків n і n +1, враховуючи (4). отримаємо:

σEx^n+1/2 ≈ σ_(i+1/2,j,k)(Exⁿ_(i+1/2,j,k)+ Exⁿ⁺¹_(i+1/2,j,k))/2,

εεo∂Ex^n+1/2 /∂t ≈ ε_(i+1/2,j,k)ε_o(Exⁿ⁺¹_(i+1/2,j,k) - Exⁿ_(i+1/2,j,k))/∆t,

μμo∂Hyⁿ /∂t ≈ μ_{(i+1/2,j,k+1/2)}μo(Hy^n+1/2_{(i+1/2,j,k+1/2)} - Hy ^n-1/2_{(i+1/2,j,k+1/2)})/∆t, (2.6)(H^n+1/2) X ≈ (Hz^n+1/2_{(i+1/2,j+1/2,k)} - Hz^n+1/2_{(i+1/2,j-1/2,k)})/ ∆y -

(Hy^n+1/2_{(i+1/2,j,k+1/2)} - Hy^n+1/2_{(i+1/2,j,k-1/2)})/∆z,(En) Y ≈ (Exⁿ_{(i+1/2,j,k+1)} - Exⁿ_(i+1/2,j,k))/ ∆z - (Ezⁿ _{(i+1,j,k+1/2)} - Ezⁿ _(i,j,k+1/2))/∆x        

^n+1/2_{(i+1/2,j,k+1/2)} = Hy^n-1/2_{(i+1/2,j,k+1/2)} + C_{Hy(i+1/2,j,k+1/2)} *((Ezⁿ _{(i+1,j,k+1/2)} -ⁿ _(i,j,k+1/2))/ ∆x - (Exⁿ_{(i+1/2,j,k+1)} - Exⁿ_(i+1/2,j,k))/ ∆z),                    (2.7)_{Hy(i+1/2,j,k+1/2)} = ∆t/ (μ_{(i+1/2,j,k+1/2)} μo)ⁿ⁺¹_(i+1/2,j,k) = C1_{Ex(i+1/2,j,k)} Exⁿ_(i+1/2,j,k) +C2_{Ex(i+1/2,j,k)} * (Hz^n+1/2_{(i+1/2,j+1/2,k)} -^n+1/2_{(i+1/2,j-1/2,k)})/ ∆y - (Hy^n+1/2_{(i+1/2,j,k+1/2)} - Hy^n+1/2_{(i+1/2,j,k-1/2)})/ ∆z) (2.8)_{Ex(i+1/2,j,k)} = (ε_(i+1/2,j,k)εo - 0,5σ_(i+1/2,j,k)∆t)/(ε_(i+1/2,j,k)εo + 0,5σ_(i+1/2,j,k)∆t)_{Ex(i+1/2,j,k)} = ∆t/(ε_(i+1/2,j,k)εo + 0,5σ_(i+1/2,j,k)∆t)

Аналогічні вирази можна отримати для решти чотирьох компонент осередку Yee.

З виразів (2.7) і (2.8) видно, що значення μ, ε та σ задаються для кожного з векторів осередку і можуть бути різними в різних напрямках. Тобто при необхідності можна задати анізотропію матеріалів для Е і/або Н полів.

Вирази (2.7) і (2.8) є достатніми для багатьох вирішуваних завдань, але для розрахунків зосереджених елементів (джерел напруги, індуктивностей, транзисторів і т.п.), а також для розрахунків матеріалів з нелінійними властивостями потрібно їх модифікація. На закінчення слід згадати, що явні кінцево-різницеві схеми вимагають спеціальних умов для сталої роботи. Для методу FDTD це умова має вигляд:

Δt ≤ 1 / (v √ ((1/Δx2) + (1/Δy2) + (1/Δz2))),

де v - максимальна швидкість електромагнітних хвиль у рахунковому обсязі, а вираз (1/Δx2) + (1/Δy2) + (1/Δz2) знаходиться під знаком квадратного кореня.

Зазвичай v = c (швидкості світла у вакуумі).

При заданих початкових умовах алгоритм Йі дає еволюційне рішення в часі від початку відліку із заданим тимчасовим кроком.

Поля в комірці сітки FDTD. З таких осередків складається просторова тривимірна сітка Йі.

Аналогічна (розділена) сітка використовується при вирішенні задач гідродинаміки (для тиску і поля швидкості).

Як і в будь-якому іншому різницевому методі, в FDTD існує проблема неточного відображення границі тіла на обчислювальну сітку. Будь крива поверхня, що розділяє сусідні середовища і геометрично не узгоджена з сіткою, буде спотворюватися ефектом «сходового наближення». Для вирішення даної проблеми можна використовувати додаткову сітку з великою роздільною здатністю в тих областях простору, де розташовані тіла зі складною геометричною структурою. Також можна видозмінювати різницеві рівняння у вузлах сітки, що знаходяться поблизу кордону між сусідніми тілами. Менш витратним методом є введення ефективної діелектричної проникності поблизу кордону між тілами (subpixel smoothing).

Чисельна схема FDTD не припускає можливості табличного завдання залежності діелектричної проникності від частоти. Однак, її можна представити у вигляді апроксимації (фітинга) членами Дебая, Друде, Лоренца або Лоренца з поглинанням. Така апроксимація не обов'язково має фізичний сенс, і може бути отримана чисельно, наприклад за допомогою програми.

Переваги та недоліки

Як і будь-який інший чисельний метод, FDTD свої переваги і недоліки.

Переваги:

·        FDTD - це простий і інтуїтивно зрозумілий метод.

·        Оскільки FDTD працює в тимчасовій області, він дозволяє отримати результат для широкого спектру довжин хвиль за один розрахунок. Це може корисно при вирішенні завдань, в яких не відомі резонансні частоти або в разі моделювання широкосмугових сигналів.

·        FDTD дозволяє створювати анімовані зображення поширення хвилі в рахунковому обсязі.

·        FDTD зручний при завданні анізотропних, дисперсних і нелінійних середовищ.

·        Метод дозволяє безпосередньо моделювати ефекти на отворах, так само як ефекти екранування, причому поля всередині і поза екраном можуть бути розраховані як безпосередньо, так і немає.

Недоліки:

·        Величина кроку дискретизації по простору повинна бути мала в порівнянні зі спектром досліджуваних довжин хвиль і характерним розміром досліджуваної структури. У деяких випадках (інверсні опали з маленькими перегородками між кульками) це може зажадати сіток з великою роздільною здатністю, що означає великі витрати пам'яті і великий час розрахунку.

·        FDTD розраховує поля всередині лічильної області. Якщо потрібно знайти поле на великій відстані від джерела, це вимагає збільшення рахункової області і часу розрахунку. Існують розширення методу для знаходження дальніх полів, але вони вимагають постобробки.

2.2    Граничні умови PEC симетрії і АВС

У рахунковому обсязі кожен вектор Е або Н обчислюється через 4 сусідніх вектора. Так відбувається по всьому об'єму. Але на кордонах самі останні вектори Е мають на гранях паралелепіпеда рахункового обсягу тільки три сусідніх вектора Н з чотирьох необхідних на ребрах - два. Тому точно обчислити поле Е на кордонах неможливо.

Проблема обчислення граничних полів вирішується різними способами.

2.2.1 Умови PEC - ідеальний провідник

Умови PEC такі, що граничні вектора Е ніколи не обчислюються, а, значить, завжди дорівнюють нулю. Як відомо, поле Е завжди дорівнює нулю в ідеальному провіднику, тому такі кордони ведуть себе як ідеальний провідник: електромагнітні хвилі 100% відбиваються назад в обрахункову область.

2.2.2 Умови симетрії

У деяких випадках поле Е або поле Н може бути симетричним щодо деякої площини. Тоді в цій площині можна задати умову симетрії і тим самим вдвічі зменшити рахунковий обсяг. При цьому поряд з даної площиною симетрії буде проходити кордон рахункового обсягу з умовами симетрії.

Симетрія може бути парній або непарній.

При непарній симетрії площину симетрії проходить всередині рахункового обсягу паралельно грані на відстані пів-осередку від грані. Умови непарній симетрії для симетрії по Е виходять простим перенесенням значень найближчих до кордону векторів Е на саму кордон, а для симетрії по Н виходять так само, але при цьому вектор Е змінює свій знак. Нижче наведені приклади непарній симетрії для Е і Н.

Умови парної симетрії дещо складніше. Площина симетрії проходить на відстані цілої осередку від кордону, тому окрім поля Е на кордоні необхідно пам'ятати про прилеглих до кордону векторах Н. При цьому симетрії по Е і Н відрізняються один від одного тільки знаками переносимих значень.

При парній симетрії граничні поля Е і Н встановлюються не одночасно, а на відповідних півкроку за часом.

2.2.3 Прості умови поглинання (АВС)

Для умов поглинання значення векторів електричного поля на кордоні обчислюються на підставі відомих полів у прикордонних шарах. Причому беруться прикордонні поля не тільки на поточному кроці по часу, але і на попередніх кроках. Всі ці умови намагаються спрогнозувати полі на кордоні.

У літературі зустрічається опис цілої низки простих умов поглинання різних авторів. Всього їх набереться з десяток. Практично найчастіше застосовують умови Мура (Mur) і Ліао (Liao). Решта не застосовують або використовують дуже рідко через їх більш низькій ефективності (Trefethen-Halpern, Higdon) або незручності у використанні (Retarded Time - RT), або через незастосовність в декартових координатах (Bayliss-Turkel), або через тенденції до втрати стабільності (відноситься майже до всіх умов, але особливо - до умов, заснованим на вищих порядках точності кінцевих різниць).

Всі умови мають досить низький коефіцієнт відбиття від кордону, що становить порядку 0,1. 1%, але тільки при падінні хвилі на межу під прямим кутом. При падінні під гострим кутом коефіцієнт відбиття зростає аж до 100% при падінні по дотичній. Через це кордону необхідно розташовувати якомога далі від джерела електромагнітних хвиль, щоб хвилі приходили до кордону під якомога більшими кутами, бажано по нормалі до кордону.

Слід зазначити, що оцінка ефективності тих чи інших простих граничних умов різна у різних авторів. Наприклад, одне джерело пише, що умови Ліао на 20 dB ефективніше умов Мура 2-го порядку. Інший пише, що умови Ліао на 12 dB гірше умов Мура 1-го порядку. Мається на увазі коефіцієнт відбиття. І обидва доводять свої висновки графіками порівняльних розрахунків.

Істина, швидше за все, посередині: для кожної конкретної задачі є свої оптимальні граничні умови. Одне погано - заздалегідь ніколи не знаєш, які краще. Можливо, введення чисел подвійної точності радикально підвищить стабільність, але поки це недозволена розкіш, витрачає в два рази більше пам'яті і часу.

Отже, граничних умов багато. Тут розглянемо три варіанти найпростіших граничних умов: Мура 1-го порядку, Ліао 3-го порядку і RT. У таблиці S (велика) і s (мале) - різні змінні! S = (dT*c)/D - число Куранта, де D - крок за простором, dT-крок по часу, с-швидкість світла, а s - координата границі. s-1 - один осередок всередину рахункового обсягу, s-2 - два осередки від кордону і т.д. На пояснюючих малюнках від кордону (surface - зовнішньої поверхні рахункового обсягу) по горизонталі відкладена координата s, а по вертикалі - час у кроках рахунку (n - часовий індекс). Формула приведена в двох видах, що доповнюють один одного в просторі і часі. На малюнках порожній кружок - обчислюване значення, чорні кружки - необхідні значення. Відразу можна визначити, які змінні прикордонній області потрібно зберігати в додаткових масивах і скільки кроків по часу вони повинні зберігатися.

Таблиця 2.1 - Граничні умови

Видно, що найбільше масивів змінних потрібно зберігати для умов RT (5), найменше - для умов Мура (1). Для умов Ліао потрібно 3 додаткових масиву змінних.

2.3    Умови PМL

Примітка. Взагалі, дослівний переклад «Perfectly Matched Layer» виглядає як «ідеально узгоджений (поєднуваний) шар». Але, по-перше, по одному шару в цих граничних умовах не застосовують, а, по-друге, якщо один шар з чимось поєднується, то, безсумнівно, з іншим шаром. Тобто їх мінімум два, а значить, по-російськи, множина - «шари».

Умови PML вперше опубліковані в статті [J.P. Berenger, JOURNAL OF COMPUTATIONAL PHYSICS. 114, 185 (1994)], а потім отримали розвиток в публікаціях цього автора в 1995-1996 р. При виникненні інтересу до безумовно-стабільним алгоритмам були розроблені умови PML, адаптовані до нових алгоритмів.

Умови PML володіють низьким коефіцієнтом відображення (за деякими даними в мільйон разів менше, ніж у умов Мура), а також практичної незалежністю від кута падіння хвилі.

До недоліків умов PML слід віднести значно більший обсяг необхідної пам'яті, ніж для умов ABC і наявність нижньої граничної частоти, для зниження якої вимагається збільшення кількості шарів PML, а, отже, необхідної пам'яті. Як наслідок збільшення необхідного обсягу пам'яті відбуватися зниження швидкості обчислень.

За межами головного рахункового обсягу додаються додаткові осередки, навколишні рахунковий обсяг по периметру декількома шарами. Електричне і магнітне поле в цих осередках обчислюється майже так само, як і в рахунковому обсязі, але є відмінності.

По-перше, в рівняння в обов'язковому порядку вводяться електричні й магнітні втрати. У головному алгоритмі цих втрат може не бути взагалі. У рівняннях головного алгоритму, враховані тільки електричні втрати. Магнітні втрати вводяться так само, як і електричні, шляхом завдання «щільності магнітних струмів». Тоді рівняння Максвелла і виглядають так:

(H) = ∂D/∂t + J(E) = - ∂B/∂t + J^*                                            (2.9)

де:

= ε ε_o E= σ E= μ μ_o H                                                    (2.10)^* = σ^* H

Відміна від (2.1) і (2.2) тільки в появі J * - щільності «магнітних струмів» і σ * - «магнітної провідності». З введенням магнітних втрат рівняння (2.9) стали симетричними.

По-друге, запроваджується поділ векторів і їх роздільне обчислення. Кожен вектор декартовій сітки Йі в межах PML ділиться на два паралельних вектора (дві компоненти). Сума цих векторів є повний вектор:

= Exy + Exz; Ey = Eyx + Eyz; Hz = Hzx + Hzy і т.д.

Позначення розшифровуються так: Exy - вектор E в напрямку X, отриманий через сусідні вектори Hz, що лежать на прямій, паралельній осі Y. Зрозуміти це можна, подивившись на систему рівнянь Беренгера:

По-третє, теоретично, якщо виконується умова:

то на межі розділу двох середовищ швидкість електромагнітних хвиль не змінюється і відображення рано нулю. У той же час, оскільки σi і σi * не дорівнюють нулю, то відбувається поглинання електромагнітних хвиль в надрах PML.

На жаль, віддзеркалення все ж є:

Від першого шару PML;

Між шарами PML, оскільки для економії обчислювальних ресурсів втрати зростають від шару до шару (закон зміни втрат від першого шару до останнього називається «профілем втрат»);

Після останнього шару PML, оскільки там знаходиться PEC - межа.

Відображення від першого шару PML і між шарами PML викликано помилками звичайно-різницевої дискретизації, і, в першу чергу, тим, що вектори E і H (а, отже, σ_i та σ_i^*) не збігаються в просторі. Для зниження відбиття усередині PML необхідно обмежувати швидкість росту втрат деяким розумною межею.

Беренгер показав, що при кожному конкретному значенні σ_i і σ_i^*, внаслідок цифрового відображення, відбувається відображення нормальних (перпендикулярних) до кордону електромагнітних хвиль, частота яких нижче fc (частота відсічення). Чим більше σ_i і σ_i^*, тим вище f_c. Тому при проникненні хвилі в межі PML спочатку йде відбиття від першого шару з провідністю σ₀ тих частот, які нижче fc0. Потім іде відображення від півшару з провідністю σ₀^* частот нижче f_c0^*. Причому f_c0 < f_c0^*. І так далі - все більш і більш високі частоти відображаються, причому відбиття від σ_i і від σ_i^*.

Відображення від PEC кордону після останнього шару PML відбувається зазвичай вже для вельми ослабленою хвилі. Відбита хвиля на зворотному шляху продовжує послаблюватися. Але якщо шарів мало (звичайно <5), то відбита хвиля може бути істотною.

Для зменшення відбиття від першого шару значення σ₁ спеціально вибирається маленьким. Для зменшення відбиття між шарами профіль втрат вибирається з обмеженою швидкістю зростання втрат. Для зменшення впливу хвилі, відбитої від РЕС кордону, збільшується кількість шарів PML.

Як варіант, Беренгер пропонує наступний геометричний профіль втрат:

де g - коефіцієнт геометричній прогресії; ∆x - крок по простору; с - швидкість світла; N - номер PML-шару, відраховуючи від інтерфейсу рахункового регіону і кордони; r - відстань від межі; R (0) - коефіцієнт відбиття від першого шару. Рекомендоване значення R (0) = 0.01 (1%) Коефіцієнт прогресії g рекомендується брати 2.15. Це значення отримано Беренгер експериментально, хоча зустрічаються й інші рекомендації.

σ^*(r) виходить через σ (r) з використанням (2.13).

На низьких частотах спостерігається різке збільшення коефіцієнта відбиття від границь PML. Нижня гранична частота відсічення для відомого значення електричної провідності на кордоні знаходиться з виразу:

Для розрахунків відгуків від імпульсів - сходинок з постійної складової зворотна величина до (2.15) - 1/f_c - це максимальний час рахунку, при якому коефіцієнт відбиття від першого шару ще не перевищує заданого R (0).

Як можна помітити, у рівняннях для векторів Е має суми типу (Hzx^{n +1/2}_{(i +1/2, j+1/2,k)} + Hzy^n+1/2_{(i+1/2, j+1/2, k)}). Це і є повний вектор, в даному випадку Hz. Коли обчислюється поле на кордоні рахункового обсягу, то з рахункового обсягу береться повний вектор Hz, а з граничного PML-шару - сума Hzx + Hzy. Також і для інших векторів.

На зовнішній межі PML, як уже говорилося, тангенціальні вектори Е дорівнюють нулю (PEC).

Профіль втрат залежить тільки від координати, що веде від інтерфейсу «рахунковий обсяг» - PML вглиб PML. На будь грані прямокутного рахункового обсягу вглиб PML веде тільки одна координата. Припустимо, це координата X. Тоді σ_x(r) змінюється по заданому закону, а σ_y(r) = σ_z (r) = 0. На ребрах вглиб PML ведуть уже дві координати (одна з σ дорівнює нулю), а в кутах вглиб ведуть всі три координати і, отже, змінюються всі σ (рис. 2.1).

Рисунок 2.1 - Розрахунковий об’єм в оточені PML шарів

3.     
Програмна реалізація

3.1 Вибір мови програмування

Серед великої кількості доступних мов програмування перевага надається мовам:

·        простий та інтуїтивно зрозумілий синтаксис;

·        можливість розділення та зручної зміни структури;

·        принципи ООП;

·        кросплаформенність;

·        забезпеченість зручного інструментарію;

·        наявність бібліотек;

·        легкість розширення.

Даним вимогам підходять чимало мов: PHP, Ryby, Python, Java, C/C++ та інші.

Проте, хотілося б, щоб отримана програма могла з легкістю використовуватись на будь-яких ОС без перезбирання чи перекомпіляції. Тому мій вибір впав на Java.- об'єктно-орієнтована мова програмування, випущена компанією Sun Microsystems у 1995 році як основний компонент платформи Java. Синтаксис мови багато в чому походить від C та C++. У офіційній реалізації, Java програми компілюються у байткод, який при виконанні інтерпретується віртуальною машиною для конкретної платформи.Microsystems надає компілятор Java та віртуальну машину Java, які задовольняють специфікації Java Community Process, під ліцезією GNU General Public License.

Передусім, Java розроблялась як платформо-незалежна мова, тому вона має менше низькорівневих можливостей для роботи з апаратним забезпеченням. За необхідності таких дій java дозволяє викликати підпрограми, написані іншими мовами програмування.

Під «незалежністю від архітектури» мається на увазі те, що програма, написана на мові Java, працюватиме на будь-якій підтримуваній апаратній чи системній платформі без змін у початковому коді та перекомпіляції.

Цього можна досягти, компілюючи початковий Java код у байт-код, який являє собою спрощені машинні команди. Потім програму можна виконати на будь-якій платформі, що має встановлену віртуальну машину Java, яка інтерпретує байткод у код, пристосований до специфіки конкретної операційної системи і процесора. Зараз віртуальні машини Java існують для більшості процесорів і операційних систем.

Стандартні бібліотеки забезпечують загальний спосіб доступу до таких платформозалежних особливостей, як обробка графіки, багатопотоковість та роботу з мережами. У деяких версіях задля збільшення продуктивності JVM, байт-код можна компілювати у машинний код до або під час виконання програми.

Основна перевага використання байт-коду - це портативність. Тим не менш, додаткові витрати на інтерпретацію означають, що інтерпретовані програми будуть майже завжди працювати повільніше, ніж скомпільовані у машинний код, і саме тому Java одержала репутацію «повільної» мови. Проте, цей розрив суттєво скоротився після введення декількох методів оптимізації у сучасних реалізаціях JVM.

Швидкість офіційної віртуальної машини Java значно покращилася з моменту випуску ранніх версій, до того ж, деякі випробування показали, що продуктивність JIT компіляторів у порівнянні зі звичайними компіляторами у машинний код майже однакова.

3.2    Попередня структура програми

Програма буде складатись, як найменше, з двох головних модулів: розрахункового та графічного.

Розрахунковий модуль повинен підключати підмодулі, що будуть працювати з обрахунками, залежними від вхідних даних та вибраних опцій. Головна задача - організувати розрахунок. Інформація про доступність та наявність відповідних модулів знаходитиметься в конфігураційному файлі, що забезпечить розширення можливостей програми без її перезбирання та правки вихідного коду.

Важливим моментом є збереження швидкодії та використання пам’яті, для цього потрібно підібрати найбільш оптимізовані структури збереження даних (масив, список і т.д.). Також, для вирішення проміжних розрахунків, будуть використовуватись математичні бібліотеки, що значно може підвищити швидкодію.

Щодо графічного модуля він повинен надати користувачеві зручність та зрозумілість ідей програми, бути максимально простим та мати певні режими налаштувань (напр. звичайний, розширений). Сюди увійдуть інтерфейс та графічне представлення розрахунків. Для графічного представлення використовуватиметься бібліотека JavaFX, що є досить потужним інструментом у відображені графіки. Для інтерфейсу буде використовуватись поєднання бібліотек Swing та JavaFX.

Висновок

різниця граничний диференціальний рівняння

В даній роботі представлено огляд прогресивного методу рішення диференційних рівнянь, який використовується в задачах електродинаміки заснованого на дискретизації рівнянь Максвелла. Розглянуто теоретичні відомості рівнянь електромагнітних полів, визначення значень при поширені хвиль та проблеми отримання значень при граничних умовах.

Складено попередню структуру програми, яка матиме можливість розширюватись.

Література

1.      John B. Schneider Understanding the Finite-Difference Time-Domain Method. - 2011.

2.      Введение в метод FDTD [Електроний ресурс] - Режим доступу: http://zfdtd.narod.ru/method/intrfdtd.htm.

3.      Антены и техника СВЧ [Електроний ресурс] - Режим доступу: http://kepstr.eltech.ru/antenn.htm.

4.      Решение динамических задач акустики методом конечных разностей во временной области [Електроний ресурс] - Режим доступу: http://fpribor.ru/uploadedFiles/files/FDTD_1D.pdf.

5.      Численное моделирование методом конечных разностей во временной области компактных акустооптических фильтров на основе многоотражательного расширения пучка [Електроний ресурс] - Режим доступу: http://www.quantum-electron.ru/pdfrus/fullt/2007/4/13329.pdf.

6.      Метод конечних різностей во временной области [Електроний ресурс] - Режим доступу: http://ru.wikipedia.org/wiki/Метод_конечных_разностей_во_временной_области.

7.      Чисельное моделирование двухмерных фотонних кристалов [Електроний ресурс] - Режим доступу: http://jre.cplire.ru/jre/nov06/3/

8.      A. Taflove and K.R. Umashankar Radar cross section of general three-dimensional scatterers. - 1983 [Електроний ресурс] - Режим доступу: http://www.ece.northwestern.edu/ecefaculty/taflove/Paper10.pdf.

9.      D. Gallagher Photonics CAD Matures. - 2008 [Електроний ресурс] - Режим доступу: http://www.photond.com/files/docs/leos_newsletter_feb08_article.pdf.

10.    FDTD (Finite-Difference Time-Domain) [Електроний ресурс] - Режим доступу: http://fdtd.kintechlab.com/ru/fdtd.

11.    A. Deinega, I. Valuev, B. Potapkin, Y. Lozovik Minimizing light reflection from dielectric textured surfaces - 2011 [Електроний ресурс] - Режим доступу:://fdtd.kintechlab.com/_media/deinega_thesis.pdf.

12.    A. Deinega, S. John Solar power conversion efficiency in modulated silicon nanowire photonic crystals - 2012 [Електроний ресурс] - Режим доступу:://fdtd.kintechlab.com/_media/deinega_-_solar_power_conversion_efficiency_in_modulated_silicon_nanowire_photonic_crystals.pdf.

13.    Ярив А., Юх. П. Оптические волны в кристаллах - Москва «Мир», 1987. - 616 с.

15.    Литвинов О.С., Павлов К.Б., Горелик В.С. Электромагнитные волны и оптика - М.: МГТУ им. Н.Э. Баумана, 2002. - 763с