Содержание
- 2. Введение в метод FDTD FDTD - "finite-difference time-domain", в русскоязычной литературе КРВО - "конечные разности во
- 3. В 1966 г. Йе (Yee) разработал технику, реализующую явную конечно - разностную схему второго порядка для
- 4. Исходными являются уравнения Максвелла в дифференциальной форме rot(H) = ∂D/∂t + J; div(B)=0; (1) rot(E) =
- 5. Для решения уравнения (1) следует выразить в декартовых координатах векторы Е и Н: Е = Ex(t,x,y,z)X+
- 6. Yee (1966) предложил пространственную сетку для конечно-разностной аппроксимации, в которую поместил вектора Ex, Ey, Ez, Hx,
- 7. Все компоненты (Ex, Ey, Ez, Hx, Hy, Hz) находятся в разных местах, т.е. разнесены в пространстве.
- 8. Поля E и H вычисляются со сдвигом на полшага по времени. Обозначения, введенные Yee, следующие: En
- 9. Поставим (3) и (2) в (1). Получим: rot(H) X = εεo∂Ex /∂t + σEx; (5) rot(E)
- 10. Подставляя (6) в (5) и решая получившиеся выражения относительно Hyn+1/2(i+1/2,j,k+1/2) и Exn+1(i+1/2,j,k) получим: Hyn+1/2(i+1/2,j,k+1/2)=Hyn-1/2(i+1/2,j,k+1/2)+CHy(i+1/2,j,k+1/2)((Ezn(i+1,j,k+1/2)-Ezn(i,j,k+1/2))/∆x- (Exn(i+1/2,j,k+1)-Exn(i+1/2,j,k))/ ∆z);
- 11. 16.11.2016 Егоров Н.Н.
- 12. Базовый алгоритм FDTD Выше кратко описан вывод конечно-разностных уравнений в декартовых координатах для двух компонент сетки
- 13. Примечания. 1. Полуцелые индексы, которые применены в нотации Yee, часто заменяются на целые путем уменьшения на
- 14. Шесть уравнений базового алгоритма очень просты в реализации. Для этого необходимо создать шесть массивов Ex, Ey,
- 15. Вычисление полей Е и Н ведется рекурсивно: для вычисления значений на текущем шаге по времени используются
- 16. Две процедуры вычисляют поля E и H: // ОПРЕДЕЛЕНИЕ МАГНИТНОГО ПОЛЯ H Procedure TimeStepForH; var I,J,K:
- 17. // ОПРЕДЕЛЕНИЕ ЭЛЕКТРИЧЕСКОГО ПОЛЯ E Procedure TimeStepForE; var I,J,K: Integer; begin for I:=1 to XSize-1 do
- 18. Массивы AE[I,J,K], BE[I,J,K] и AН[I,J,K] - это предварительно вычисленные коэффициенты С1, С2 и С соответственно из
- 19. Приведенные процедуры вызываются в цикле по времени: // ГЛАВНЫЙ ЦИКЛ ПО ВРЕМЕНИ T:=0; //Исходное значение времени
- 20. Граничные условия для FDTD В счетном объеме каждый вектор Е или Н вычисляется через 4 соседних
- 21. 2. Условия симметрии. В некоторых случаях поле Е или поле Н может быть симметричным относительно некоторой
- 22. Пример: нечетная симметрия по Е: // ПЛОСКОСТЬ Y=1/2 симметрична по E For I:=1 to NX-1 do
- 23. Условия четной симметрии несколько сложнее. Плоскость симметрии проходит на расстоянии целой ячейки от границы, поэтому кроме
- 24. Простые условия поглощения (АВС) Для условий поглощения значения векторов электрического поля на границе вычисляются на основании
- 25. Итак, граничных условий много. Здесь рассмотрим три варианта простейших граничных условий: Мура 1-го порядка, Лиао 3-го
- 26. 16.11.2016 Егоров Н.Н.
- 27. Видно, что больше всего массивов переменных нужно хранить для условий RT (5), меньше всего - для
- 28. //ВЕРХНЯЯ И НИЖНЯЯ ПЛОСКОСТИ // Ex для двух плоскостей XxZ [верхней и нижней] for K:=2 to
- 29. // Ex для двух плоскостей XxY [ближней и дальней] for I:=1 to NX-1 do for J:=2
- 30. // сохранение приграничных полей для MursFirstABC Procedure GetBorderValues; var I,J,K: Integer; begin // слева и справа
- 31. 4. Условия PМL - идеально сочетающиеся слои "Perfectly Matched Layer" выглядит как "идеально согласованный (сочетающийся) слой“
- 32. Во-первых, в уравнения в обязательном порядке вводятся электрические и магнитные потери. В главном алгоритме этих потерь
- 33. (11) 16.11.2016 Егоров Н.Н.
- 34. Примечание Схему Беренгера называют «раздельной», т.к. вводится разделение векторов. Есть также предложения других авторов применять «нераздельную»
- 35. В третьих, теоретически, если выполняется условие σi/ εo = σi*/ μo, (12) то на границе раздела
- 36. Отражение от PEC границы после последнего слоя PML происходит обычно уже для весьма ослабленной волны. Отраженная
- 37. На низких частотах наблюдается резкое увеличение коэффициента отражения от границ PML. Нижняя граничная частота отсечки для
- 38. Как можно заметить, в уравнениях для векторов Е имеются суммы типа (Hzxn+1/2(i+1/2,j+1/2,k) + Hzyn+1/2 (i+1/2,j+1/2,k)). Это
- 39. Рис. 1. Счетный объем в окружении слоев PML 16.11.2016 Егоров Н.Н.
- 40. Расширения базового алгоритма Существует множество дополнений к базовому алгоритму, расширяющих его возможности. Здесь представлены два вида
- 41. Проведя вывод дискретного уравнения так, как это было во введении, но с учетом тока через сосредоточенный
- 42. Резистор Ток резистора находится по известной формуле I=U/R. В нашем случае U=∆Z(En+1+En)/2, т.е. берется среднее арифметическое
- 43. Резистивный источник напряжения Это резистор, в котором встроен источник напряжения. Очень удобный и абсолютно физичный элемент.
- 44. Конденсатор Ток через конденсатор I=C(du/dt), где С – емкость. Производную напряжения по времени выразим через разность
- 45. Индуктивность Ток в индуктивности, как известно, равен I=(1/L)∫Udt, т.е. величина интегральная. В дискретной форме ток запишем
- 46. Расширения базового алгоритма Полярные диэлектрики Наличие относительной диэлектрической проницаемости большей, чем единица, связано с поляризацией молекул
- 47. Применение метода FDTD для расчетов распространения нестационарных электромагнитных полей в условиях частотной дисперсии В классической постановке
- 48. В [2] 2. Luebbers R. et all.: "A frequency-depended finite-difference time-domain formulation for dispersive materials." IEEE
- 49. Учет проводимости биологических тканей в уравнениях (FD2)TD Предположим, что биологические материалы линейные и изотропные, магнитная проницаемость
- 50. Вектор смещения из (1) во временной области записывается как [2, 3]: (3) где - диэлектрическая постоянная,
- 51. Принимая во внимание, что D(t) и E(t) равны нулю при t (6) Выражение для D на
- 52. Для удобства введем обозначения: (9) Подставляя (9) в (8) и решая последнее уравнение относительно получаем: (10)
- 53. Если относительная диэлектрическая проницаемость не зависит от частоты, то χm=0 для всех m и χ0=0. В
- 54. Выражение (10) содержит член который, на первый взгляд, для вычисления требует хранения большого количества переменных. Но,
- 55. Вычисление электромагнитного поля в дальней зоне с использованием метода FDTD и интеграла Кирхгофа Решение уравнений Максвелла
- 56. Преобразование ближнего поля в дальнее поле с использованием интеграла Кирхгофа Существует несколько методов преобразования ближнего поля
- 57. 16.11.2016 Егоров Н.Н.
- 58. Формула (6) выражает принцип Гюйгенса, согласно которому каждая точка на волновом фронте служит фиктивным источником воображаемой
- 59. Предположим, что поверхность dA’ перпендикулярна оси Z и точка p’ имеет координаты (i, j, ko) (рис.
- 60. Для упрощения дальнейшей записи введем обозначения: (8) и применим стандартные для FDTD обозначения для дискретных значений
- 61. В (9) значение функции Ψ(p, tn*) вычисляется с использованием значений Ψ на трех шагах по времени.
- 63. Скачать презентацию
Введение в метод FDTD
FDTD - "finite-difference time-domain",
в русскоязычной литературе КРВО
Введение в метод FDTD
FDTD - "finite-difference time-domain",
в русскоязычной литературе КРВО
метод FDTD - один из многочисленных методов решения дифференциальных уравнений, но среди тех, кто занимается решением задач электродинамики, FDTD является синонимом решения вихревых дифференциальных уравнений Максвелла.
16.11.2016
Егоров Н.Н.
В 1966 г. Йе (Yee) разработал технику, реализующую явную конечно -
В 1966 г. Йе (Yee) разработал технику, реализующую явную конечно -
K. S. Yee, “Numerical Solution of Initial Boundary Value Problems Involving Maxwell’s Equations in Isotropic Media,” IEEE Transactions on Antennas and Propagation, vol. 14, pp. 302-307, May 1966.
16.11.2016
Егоров Н.Н.
Исходными являются уравнения Максвелла в дифференциальной форме
rot(H) = ∂D/∂t + J;
Исходными являются уравнения Максвелла в дифференциальной форме
rot(H) = ∂D/∂t + J;
rot(E) = - ∂B/∂t; div(D)=ρ;
D = ε εo E;
J = σ E; (2)
B = μ μo H
Два уравнения (1) содержат пространственные и временные производные.
16.11.2016
Егоров Н.Н.
Для решения уравнения (1) следует выразить в декартовых координатах векторы Е
Для решения уравнения (1) следует выразить в декартовых координатах векторы Е
Е = Ex(t,x,y,z)X+ Ey(t,x,y,z)Y+ Ez(t,x,y,z)Z; (3)
H = Hx(t,x,y,z)X+ Hy(t,x,y,z)Y+ Hz(t,x,y,z)Z;
где Ex, Ey, Ez, Hx, Hy, Hz - проекции векторов на координатные оси, а X, Y, Z - единичные векторы.
Остальные величины в (1) - D, B, J - выразим через E и H. Величины E и H для нас будут основными.
Примечание: существуют и другие подходы, когда в уравнениях (1) вначале оставляют D и/или B, но в конце концов всё равно выражаются вектора Е и Н. Также следует указать, что уравнения (1) записаны не полностью. Например, в них не учитываются сторонние токи.
16.11.2016
Егоров Н.Н.
Yee (1966) предложил пространственную сетку для конечно-разностной аппроксимации, в которую поместил
Yee (1966) предложил пространственную сетку для конечно-разностной аппроксимации, в которую поместил
16.11.2016
Егоров Н.Н.
Все компоненты (Ex, Ey, Ez, Hx, Hy, Hz) находятся в разных
Все компоненты (Ex, Ey, Ez, Hx, Hy, Hz) находятся в разных
Е-компоненты находятся посередине ребер,
Н - компоненты - по центру граней.
Все компоненты независимы друг от друга, т.е. каждой из них можно присвоить свои уникальные электрические (для Е) и магнитные (для Н) параметры.
Пространственные координаты каждого вектора x, y и z выражаются в номерах ячеек i, j и k соответственно, время t выражается в шагах n по времени:
x = i∆x; y = j∆y; z= k∆z; t= n∆t; (4)
где ∆x, ∆y, ∆z - размеры пространственной ячейки, ∆t - шаг по времени.
16.11.2016
Егоров Н.Н.
Поля E и H вычисляются со сдвигом на полшага по времени.
Поля E и H вычисляются со сдвигом на полшага по времени.
En - значение поля E на только что вычисленном шаге;
En+1 - значение поля E на вычисляемом сейчас шаге по времени.
Hn-1/2 - значение поля H на только что вычисленном шаге;
Hn+1/2 - значение поля на вычисляемом сейчас полушаге по времени.
Из этих обозначений следует, что процедура вычислений начинается с поля Hn+1/2, потому что в момент t=0 (n=0) установлены начальные условия по всему счетному объему: все значения полей E и H равны нулю. Хотя в принципе это лишь наиболее распространенная условность. Можно считать, что пространственная сетка проходит через вектора H, что процедура счета начинается с поля E.
Теперь, когда введены основные обозначения, покажем вывод выражений, пригодных для расчетов с помощью компьютера и которым уже 50 лет.
16.11.2016
Егоров Н.Н.
Поставим (3) и (2) в (1). Получим:
rot(H) X = εεo∂Ex /∂t
Поставим (3) и (2) в (1). Получим:
rot(H) X = εεo∂Ex /∂t
rot(E) Y = - μμo∂Hy /∂t;
Применяя конечно-разностную аппроксимацию, преобразуем (5) в выражения для шагов n и n+1, учитывая (4), получим:
σExn+1/2 ≈ σ(i+1/2,j,k)(Exn(i+1/2,j,k)+ Exn+1(i+1/2,j,k))/2;
εεo∂Exn+1/2 /∂t ≈ ε(i+1/2,j,k)εo(Exn+1(i+1/2,j,k) - Exn(i+1/2,j,k))/∆t;
μμo∂Hyn/∂t≈μ(i+1/2,j,k+1/2)μo(Hy n+1/2(i+1/2,j,k+1/2)-Hyn-1/2(i+1/2,j,k+1/2))/∆t; (6)
rot(Hn+1/2)X≈(Hzn+1/2(i+1/2,j+1/2,k)-Hzn+1/2(i+1/2,j-1/2,k))/∆y-(Hyn+1/2(i+1/2,j,k+1/2)-
Hyn+1/2(i+1/2,j,k-1/2))/∆z;
rot(En)Y≈(Exn(i+1/2,j,k+1)-Exn(i+1/2,j,k))/∆z-(Ezn(i+1,j,k+1/2)-Ezn(i,j,k+1/2))/∆x;
16.11.2016
Егоров Н.Н.
Подставляя (6) в (5) и решая получившиеся выражения относительно Hyn+1/2(i+1/2,j,k+1/2) и
Подставляя (6) в (5) и решая получившиеся выражения относительно Hyn+1/2(i+1/2,j,k+1/2) и
Hyn+1/2(i+1/2,j,k+1/2)=Hyn-1/2(i+1/2,j,k+1/2)+CHy(i+1/2,j,k+1/2)((Ezn(i+1,j,k+1/2)-Ezn(i,j,k+1/2))/∆x-
(Exn(i+1/2,j,k+1)-Exn(i+1/2,j,k))/ ∆z);
CHy(i+1/2,j,k+1/2) = ∆t/ (μ(i+1/2,j,k+1/2) μo); (7)
Exn+1(i+1/2,j,k)=C1Ex(i+1/2,j,k)Exn(i+1/2,j,k)+C2Ex(i+1/2,j,k)(Hzn+1/2(i+1/2,j+1/2,k)-
Hzn+1/2(i+1/2,j-1/2,k))/∆y-(Hyn+1/2(i+1/2,j,k+1/2)-Hyn+1/2(i+1/2,j,k-1/2))/∆z); (8)
C1Ex(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);
C2Ex(i+1/2,j,k)=∆t/(ε(i+1/2,j,k)εo+0,5σ(i+1/2,j,k)∆t).
Аналогичные выражения можно получить для остальных четырех компонент ячейки Yee.
Из выражений (7) и (8) видно, что значения μ, ε и σ задаются для каждого их векторов ячейки и могут быть различными в разных направлениях. Т.е. при необходимости можно задать анизотропию материалов для Е и/или Н полей.
16.11.2016
Егоров Н.Н.
16.11.2016
Егоров Н.Н.
16.11.2016
Егоров Н.Н.
Базовый алгоритм FDTD
Выше кратко описан вывод конечно-разностных уравнений в декартовых координатах
Базовый алгоритм FDTD
Выше кратко описан вывод конечно-разностных уравнений в декартовых координатах
На рис. 2. приводится полная система для всех векторов сетки Yeе
16.11.2016
Егоров Н.Н.
Примечания.
1. Полуцелые индексы, которые применены в нотации Yee, часто заменяются на
Примечания.
1. Полуцелые индексы, которые применены в нотации Yee, часто заменяются на
2. Переменные ε, σ и μ повсюду должны быть представлены как ε(i,j,k) σ(i,j,k) μ(i,j,k), но на рисунке индексы (i,j,k) опущены для экономии места.
3. ε и μ здесь - это абсолютные проницаемости, т.е. сразу готовое произведение εεo и μμo, как это представлено ранее (см. введение). Путаница в этих обозначениях - широко распространенное явление в литературе, поэтому надо держать ухо востро.
4. Коэффициенты С1 и С2 из (8), на первый взгляд, не такие, как на этом рисунке. На самом деле это то же самое. Попробуйте разделить числитель и знаменатель коэффициентов С1 и С2 (8) на (ε(i+1/2,j,k)εo).
16.11.2016
Егоров Н.Н.
Шесть уравнений базового алгоритма очень просты в реализации. Для этого необходимо
Шесть уравнений базового алгоритма очень просты в реализации. Для этого необходимо
Чем больше возможностей мы хотим реализовать в базовом алгоритме, тем больший объем памяти требуется для переменных и, соответственно, больше затраты времени при расчетах.
Раньше значительное ускорение вычислений можно было получить создав еще ряд массивов и поместив в них предварительно вычисленные коэффициенты С, С1 и С2 из уравнений (7) и (8). Но сейчас в гигагерцовых процессорах скорость вычислений арифметики с плавающей точкой выросла больше, чем скорость выборки дополнительных элементов массива из памяти, поэтому смысл в предварительном вычислении этих коэффициентов снизился.
В программе ZFDTD пожертвовали скоростью во имя экономии памяти, хотя в ранних версиях программы было наоборот.
16.11.2016
Егоров Н.Н.
Вычисление полей Е и Н ведется рекурсивно: для вычисления значений на
Вычисление полей Е и Н ведется рекурсивно: для вычисления значений на
В общем, по отношению к оперативной памяти алгоритм FDTD довольно прожорлив, но для больших задач этот алгоритм, по литературным данным, требует меньше памяти, чем интегральные методы.
Ход вычислений рассмотрим на примере, написанном на языке Pascal
16.11.2016
Егоров Н.Н.
Две процедуры вычисляют поля E и H:
// ОПРЕДЕЛЕНИЕ МАГНИТНОГО ПОЛЯ H
Procedure
Две процедуры вычисляют поля E и H:
// ОПРЕДЕЛЕНИЕ МАГНИТНОГО ПОЛЯ H
Procedure
var I,J,K: Integer;
begin
For I:=1 to XSize do
for J:=1 to YSize-1 do
for K:=1 to ZSize-1 do
HX[I,J,K]:=HX[I,J,K]-((EZ[I,J+1,K]-EZ[I,J,K])*TY-(EY[I,J,K+1]-EY[I,J,K])*TZ)*AH[I,J,K];
For I:=1 to XSize-1 do
for J:=1 to YSize do
for K:=1 to ZSize-1 do
HY[I,J,K]:=HY[I,J,K]-((EX[I,J,K+1]-EX[I,J,K])*TZ-(EZ[I+1,J,K]-EZ[I,J,K])*TX)*AH[I,J,K];
For I:=1 to XSize-1 do
for J:=1 to YSize-1 do
for K:=1 to ZSize do
HZ[I,J,K]:=HZ[I,J,K]-((EY[I+1,J,K]-EY[I,J,K])*TX-(EX[I,J+1,K]-EX[I,J,K])*TY)*AH[I,J,K];
end;
// КОНЕЦ ОПРЕДЕЛЕНИЯ МАГНИТНОГО ПОЛЯ H
16.11.2016
Егоров Н.Н.
// ОПРЕДЕЛЕНИЕ ЭЛЕКТРИЧЕСКОГО ПОЛЯ E
Procedure TimeStepForE;
var I,J,K: Integer;
begin
for I:=1 to
// ОПРЕДЕЛЕНИЕ ЭЛЕКТРИЧЕСКОГО ПОЛЯ E
Procedure TimeStepForE;
var I,J,K: Integer;
begin
for I:=1 to
EX[I,J,K]:=EX[I,J,K]*AE[I,J,K]+((HZ[I,J,K]-HZ[I,J-1,K])*TY-(HY[I,J,K]-HY[I,J,K-1])*TZ)*BE[I,J,K];
for I:=2 to XSize-1 do for J:=1 to YSize-1 do for K:=2 to ZSize-1 do
EY[I,J,K]:=EY[I,J,K]*AE[I,J,K]+((HX[I,J,K]-HX[I,J,K-1])*TZ- (HZ[I,J,K]-HZ[I-1,J,K])*TX)*BE[I,J,K]
for I:=2 to XSize-1 do for J:=2 to YSize-1 do for K:=1 to ZSize-1 do
EZ[I,J,K]:=EZ[I,J,K]*AE[I,J,K]+((HY[I,J,K]-HY[I-1,J,K])*TX- (HX[I,J,K]-HX[I,J-1,K])*TY)*BE[I,J,K];
end;
// КОНЕЦ ОПРЕДЕЛЕНИЯ ЭЛЕКТРИЧЕСКОГО ПОЛЯ E
16.11.2016
Егоров Н.Н.
Массивы AE[I,J,K], BE[I,J,K] и AН[I,J,K] - это предварительно вычисленные коэффициенты С1,
Массивы AE[I,J,K], BE[I,J,K] и AН[I,J,K] - это предварительно вычисленные коэффициенты С1,
Приведенный код вычисления полей рассчитан на такое представление объектов, когда считается, что электрофизические характеристики пространства в ячейке (i,j,k) задаются одинаковыми для всех трех Нx,y,z(i,j,k) и трех Ex,y,z(i,j,k) векторов. При таком представлении точность вычислений на границах раздела сред снижается до первого порядка. Обычно это проявляется в снижении на 25..40% вычисленной резонансной частоты, такого же увеличения вычисленных токов в длинных структурах, снижении в разы скорости затухания свободных колебаний и в других ошибках.
16.11.2016
Егоров Н.Н.
Приведенные процедуры вызываются в цикле по времени:
// ГЛАВНЫЙ ЦИКЛ ПО ВРЕМЕНИ
T:=0;
Приведенные процедуры вызываются в цикле по времени:
// ГЛАВНЫЙ ЦИКЛ ПО ВРЕМЕНИ
T:=0;
for I:=0 to <заданное количество шагов по времени> do
begin
TimeStepForH; //Вычисление H по всему объему
BordersSymmetryH; //Выполнение граничных условий для симметрии по Н
TimeStepForE; //Вычисление E по всему объему
BordersABC; // Граничные УСЛОВИЯ ПОГЛОЩЕНИЯ
BordersSymmetryE;//Выполнение граничных условий для симметрии по Е
BordersРЕС;//Выполнение граничных условий РЕС
T:=T+dT; // инкремент времени на один шаг
end;
//ГЛАВНЫЙ ЦИКЛ ПО ВРЕМЕНИ
это самое простое место в программе
16.11.2016
Егоров Н.Н.
Граничные условия для FDTD
В счетном объеме каждый вектор Е или Н
Граничные условия для FDTD
В счетном объеме каждый вектор Е или Н
Проблема вычисления граничных полей решается различными способами.
1. Условия PEC – (photo electrochemical cell) идеальный проводник.
Условия PEC таковы, что граничные вектора Е никогда не вычисляются, а, значит, всегда равны нулю. Как известно, поле Е всегда равно нулю в идеальном проводнике, поэтому такие границы ведут себя как идеальный проводник: электромагнитные волны 100 % отражаются обратно в счетный объем.
16.11.2016
Егоров Н.Н.
2. Условия симметрии.
В некоторых случаях поле Е или поле Н может
2. Условия симметрии.
В некоторых случаях поле Е или поле Н может
Симметрия может быть четной или нечетной.
При нечетной симметрии плоскость симметрии проходит внутри счетного объема параллельно грани на расстоянии пол-ячейки от грани. Условия нечетной симметрии для симметрии по Е получаются простым переносом значений ближайших к границе векторов Е на саму границу, а для симметрии по Н получаются так же, но при этом вектор Е меняет свой знак. Ниже приведены примеры нечетной симметрии для Е и Н.
16.11.2016
Егоров Н.Н.
Пример: нечетная симметрия по Е:
// ПЛОСКОСТЬ Y=1/2 симметрична по E
Пример: нечетная симметрия по Е:
// ПЛОСКОСТЬ Y=1/2 симметрична по E
For K:=1 to NZ do
EX[I,1,K]:=EX[I,2,K];
For I:=1 to NX do
For K:=1 to NZ-1 do
EZ[I,1,K]:=EZ[I,2,K];
Пример: нечетная симметрия по Н:
// ПЛОСКОСТЬ Y=1/2 симметрична по H
For I:=1 to NX-1 do
For K:=1 to NZ do
EX[I,1,K]:=-EX[I,2,K];
For I:=1 to NX do
For K:=1 to NZ-1 do
EZ[I,1,K]:=-EZ[I,2,K];
16.11.2016
Егоров Н.Н.
Условия четной симметрии несколько сложнее. Плоскость симметрии проходит на расстоянии целой
Условия четной симметрии несколько сложнее. Плоскость симметрии проходит на расстоянии целой
Пример: четная симметрия по Н:
// ПЛОСКОСТЬ Y=1 симметрична по H, вычисление поля Е
begin
For K:=1 to NZ do For I:=1 to NX-1 do EX[I,1,K]:=-EX[I,3,K]; For K:=1 to NZ-1 do For I:=1 to NX do EZ[I,1,K]:=-EZ[I,3,K];
end;
// ПЛОСКОСТЬ Y=1 симметрична по H, вычисление поля Н
begin
For K:=1 to NZ-1 do For I:=1 to NX do HX[I,1,K]:=HX[I,2,K];
For K:=1 to NZ do For I:=1 to NX-1 do HZ[I,1,K]:=HZ[I,2,K];
end;
Пример: четная симметрия по Е:
// ПЛОСКОСТЬ Y=1 симметрична по Е, вычисление поля Е
begin
For K:=1 to NZ do For I:=1 to NX-1 do EX[I,1,K]:=EX[I,3,K];
For K:=1 to NZ-1 do For I:=1 to NX do EZ[I,1,K]:=EZ[I,3,K];
end;
// ПЛОСКОСТЬ Y=1 симметрична по Е, вычисление поля Н
begin
For K:=1 to NZ-1 do For I:=1 to NX do HX[I,1,K]:=-HX[I,2,K];
For K:=1 to NZ do For I:=1 to NX-1 do HZ[I,1,K]:=-HZ[I,2,K];
end;
При четной симметрии граничные поля Е и Н устанавливаются не одновременно, а на соответствующих полушагах по времени.
16.11.2016
Егоров Н.Н.
Простые условия поглощения (АВС)
Для условий поглощения значения векторов электрического поля на
Простые условия поглощения (АВС)
Для условий поглощения значения векторов электрического поля на
В литературе встречается описание целого ряда простых условий поглощения разных авторов. Всего их наберется с десяток. Практически чаще всего применяют условия Мура (Mur) и Лиао (Liao). Остальные не применяют или используют очень редко из-за их более низкой эффективности (Trefethen-Halpern, Higdon) или неудобства в использовании (Retarded Time - RT), или из-за неприменимости в декартовых координатах (Bayliss-Turkel), или из-за тенденции к потере стабильности (относится почти ко всем условиям, но особо - к условиям, основанным на высших порядках точности конечных разностей).
Все условия имеют довольно низкий коэффициент отражения от границы, составляющий порядка 0,1..1 %, но только при падении волны на границу под прямым углом. При падении под острым углом коэффициент отражения растет вплоть до 100 % при падении по касательной. Из-за этого границы необходимо располагать как можно дальше от источника электромагнитных волн, чтобы волны приходили к границе под как можно большими углами, желательно по нормали к границе.
Следует отметить, что оценка эффективности тех или иных простых граничных условий различна у разных авторов. Например, один источник пишет, что условия Лиао на 20 dB эффективнее условий Мура 2-го порядка. Другой пишет, что условия Лиао на 12 dB хуже условий Мура 1-го порядка. Имеется ввиду коэффициент отражения. И оба доказывают свои выводы графиками сравнительных расчетов.
Истина, скорее всего, посередине: для каждой конкретной задачи имеются свои оптимальные граничные условия. Одно плохо - заранее никогда не знаешь, КАКИЕ лучше.
Примечание. Возможно, введение чисел двойной точности радикально повысит стабильность, но пока это непозволительная роскошь, расходующая в два раза больше памяти и времени.
16.11.2016
Егоров Н.Н.
Итак, граничных условий много. Здесь рассмотрим три варианта простейших граничных условий:
Итак, граничных условий много. Здесь рассмотрим три варианта простейших граничных условий:
16.11.2016
Егоров Н.Н.
16.11.2016
Егоров Н.Н.
16.11.2016
Егоров Н.Н.
Видно, что больше всего массивов переменных нужно хранить для условий RT
Видно, что больше всего массивов переменных нужно хранить для условий RT
А теперь пример реализации условий Мура 1-го порядка.
Procedure MursFirstABC;
var I,J,K: Integer;
begin
// УСЛОВИЯ ПОГЛОЩЕНИЯ. Массивы вида EZX1, EZXN и т.п. - массивы сохраненных на предыдущем шаге по времени значений поля Е в приграничной области. Процедура сохранения приводится ниже.
// коэффициенты Tx, Ty и Tz равны c*dT/dX, c*dT/dY и c*dT/dZ соответственно. c-скорость света, dT- шаг по времени; dX, dY, dZ - шаги по пространству.
//ЛЕВАЯ И ПРАВАЯ ПЛОСКОСТИ
// Ez для двух плоскостей YxZ (левой и правой)
for J:=2 to NY-1 do for K:=1 to NZ-1 do begin
EZ[1,J,K]:=EZX1[J,K]+(Tx-1)/(Tx+1)*(EZ[2,J,K]-EZ[1,J,K]);
EZ[NX,J,K]:=EZXN[J,K]+(Tx-1)/(Tx+1)*(EZ[NX-1,J,K]-EZ[NX,J,K]);
end;
// Ey для двух плоскостей YxZ [левой и правой]
for J:=1 to NY-1 do for K:=2 to NZ-1 do begin
EY[1,J,K]:=EYX1[J,K]+(Tx-1)/(Tx+1)* (EY[2,J,K]-EY[1,J,K]);
EY[NX,J,K]:=EYXN[J,K]+(Tx-1)/(Tx+1)* (EY[NX-1,J,K] -EY[NX,J,K]);
end;
16.11.2016
Егоров Н.Н.
//ВЕРХНЯЯ И НИЖНЯЯ ПЛОСКОСТИ
// Ex для двух плоскостей XxZ [верхней
//ВЕРХНЯЯ И НИЖНЯЯ ПЛОСКОСТИ
// Ex для двух плоскостей XxZ [верхней
for K:=2 to NZ-1 do for I:=1 to NX-1 do begin
EX[I,1,K]:= EXY1[I,K]+(Ty-1)/(Ty+1)*(EX[I,2,K]-EX[I,1,K]);
EX[I,NY,K]:= EXYN[I,K]+(Ty-1)/(Ty+1)*(EX[I,NY-1,K]-EX[I,NY,K]);
end;
// Ez для двух плоскостей XxZ [верхней и нижней]
for K:=1 to NZ-1 do for I:=2 to NX-1 do begin
EZ[I,1,K]:= EZY1[I,K]+(Ty-1)/(Ty+1)*(EZ[I,2,K]-EZ[I,1,K]);
EZ[I,NY,K]:= EZYN[I,K]+(Ty-1)/(Ty+1)*(EZ[I,NY-1,K]-EZ[I,NY,K]);
end;
//БЛИЖНЯЯ И ДАЛЬНЯЯ ПЛОСКОСТИ
// Ey для двух плоскостей XxY [ближней и дальней]
for I:=2 to NX-1 do for J:=1 to NY-1 do begin
EY[I,J,NZ]:=EYZN[I,J]+(Tz-1)/(Tz+1)*(EY[I,J,NZ-1]-EY[I,J,NZ]);
EY[I,J,1]:=EYZ1[I,J]+(Tz-1)/(Tz+1)*(EY[I,J,2]-EY[I,J,1]);
end
16.11.2016
Егоров Н.Н.
// Ex для двух плоскостей XxY [ближней и дальней]
for I:=1
// Ex для двух плоскостей XxY [ближней и дальней]
for I:=1
EX[I,J,NZ]:= EXZN[I,J]+(Tz-1)/(Tz+1)*(EX[I,J,NZ-1]-EX[I,J,NZ]);
EX[I,J,1]:= EXZ1[I,J]+(Tz-1)/(Tz+1)*(EX[I,J,2]-EX[I,J,1]);
end;
// РЕБРА, параллельные оси Z
for K:=1 to NZ-1 do begin
EZ[1,1,K]:=EZX1[1,K]+(Tx-1)/(Tx+1)*(EZ[2,1,K]-EZ[1,1,K]);
EZ[1,NY,K]:=EZX1[NY,K]+(Tx-1)/(Tx+1)*(EZ[2,NY,K]-EZ[1,NY,K]);
EZ[NX,1,K]:=EZXN[1,K]+(Tx-1)/(Tx+1)*(EZ[NX-1,1,K]-EZ[NX,1,K]);
EZ[NX,NY,K]:=EZXN[NY,K]+(Tx-1)/(Tx+1)*(EZ[NX-1,NY,K]-EZ[NX,NY,K]);
end;
// РЕБРА, параллельные оси Y
for J:=1 to NY-1 do begin
EY[1,J,1]:=EYX1[J,1]+(Tx-1)/(Tx+1)*(EY[2,J,1]-EY[1,J,1]);
EY[1,J,NZ]:=EYX1[J,NZ]+(Tx-1)/(Tx+1)*(EY[2,J,NZ]-EY[1,J,NZ]);
EY[NX,J,1]:=EYXN[J,1]+(Tx-1)/(Tx+1)*(EY[NX-1,J,1]-EY[NX,J,1]);
EY[NX,J,NZ]:=EYXN[J,NZ]+(Tx-1)/(Tx+1)*(EY[NX-1,J,NZ]-EY[NX,J,NZ]);
end;
// РЕБРА, параллельные оси X
for I:=1 to NX-1 do begin
EX[I,1,1]:=EXY1[I,1]+ (Ty-1)/(Ty+1)*(EX[I,2,1]-EX[I,1,1]);
EX[I,1,NZ]:=EXYN[I,NZ]+(Ty-1)/(Ty+1)*(EX[I,2,NZ]-EX[I,1,NZ]);
EX[I,NY,1]:=EXY1[I,1]+(Ty-1)/(Ty+1)*(EX[I,NY-1,1]-EX[I,NY,1]);
EX[I,NY,NZ]:=EXYN[I,NZ]+(Ty-1)/(Ty+1)*(EX[I,NY-1,NZ]-EX[I,NY,NZ]);
end;
end;
16.11.2016
Егоров Н.Н.
// сохранение приграничных полей для MursFirstABC
Procedure GetBorderValues;
var I,J,K: Integer;
begin
//
// сохранение приграничных полей для MursFirstABC
Procedure GetBorderValues;
var I,J,K: Integer;
begin
//
for J:=1 to NY do for k:=1 to NZ-1 do begin
EZX1[J,K]:=EZ[2,J,K];
EZXN[J,K]:=EZ[NX-1,J,K]; end;
for J:=1 to NY-1 do for k:=1 to NZ do begin
EYX1[J,K]:=EY[2,J,K];
EYXN[J,K]:=EY[NX-1,J,K]; end;
// верх и низ
for I:=1 to NX do for k:=1 to NZ-1 do begin
EZY1[I,K]:=EZ[I,2,K];
EZYN[I,K]:=EZ[I,NY-1,K]; end;
for I:=1 to NX-1 do for k:=1 to NZ do begin
EXY1[I,K]:=EX[I,2,K];
EXYN[I,K]:=EX[I,NY-1,K]; end;
// ближняя и дальняя
for I:=1 to NX do for J:=1 to NY-1 do begin
EYZ1[I,J]:=EY[I,J,2];
EYZN[I,J]:=EY[I,J,NZ-1]; end;
for I:=1 to NX-1 do for J:=1 to NY do begin
EXZ1[I,J]:=EX[I,J,2];
EXZN[I,J]:=EX[I,J,NZ-1]; end;
end;
Примечание: в литературе встречается и несколько иная форма записи условий Мура 1-го порядка. К приведенной здесь формуле добавляется еще один член, который сам Мур относит к условиям второго порядка, а другие авторы утверждают, что этот член должен быть и в формуле 1-го порядка.
16.11.2016
Егоров Н.Н.
4. Условия PМL - идеально сочетающиеся слои
"Perfectly Matched Layer" выглядит как
4. Условия PМL - идеально сочетающиеся слои
"Perfectly Matched Layer" выглядит как
Условия PML обладают низким коэффициентом отражения (по некоторым данным в миллион раз меньше, чем у условий Мура), а также практической независимостью от угла падения волны.
К недостаткам условий PML следует отнести значительно больший объем требуемой памяти, чем для условий ABC и наличие нижней граничной частоты, для снижения которой требуется увеличения количества слоев PML, а, следовательно, требуемой памяти. Как следствие увеличения требуемого объема памяти происходить снижение скорости вычислений.
За пределами главного счетного объема добавляются дополнительные ячейки, окружающие счетный объем по периметру несколькими слоями. Электрическое и магнитное поле в этих ячейках вычисляется почти так же, как и в счетном объеме, но есть отличия.
16.11.2016
Егоров Н.Н.
Во-первых, в уравнения в обязательном порядке вводятся электрические и магнитные потери.
Во-первых, в уравнения в обязательном порядке вводятся электрические и магнитные потери.
rot(H) = ∂D/∂t + J;
rot(E) = - ∂B/∂t + J*; (9)
где D = ε εo E;
J = σ E;
B = μ μo H; (10)
J* = σ* H;
Отличие от (1) и (2) только в появлении J* - плотности "магнитных токов" и σ* - "магнитной проводимости". С введением магнитных потерь уравнения (9) стали симметричными.
Во-вторых, вводится разделение векторов и их раздельное вычисление. Каждый вектор декартовой сетки Yee в границах PML делится на два параллельных вектора (две компоненты). Сумма этих векторов есть полный вектор: Eх=Exy+Exz; Ey=Eyx+Eyz; Hz=Hzx+Hzy и т.д. Обозначения расшифровываются так: Exy - вектор E в направлении X, полученный через соседние векторы Hz, лежащие на прямой, параллельной оси Y. Понять это можно, посмотрев на систему уравнений Беренгера:
(11)
16.11.2016
Егоров Н.Н.
(11)
16.11.2016
Егоров Н.Н.
(11)
16.11.2016
Егоров Н.Н.
Примечание
Схему Беренгера называют «раздельной», т.к. вводится разделение векторов. Есть также предложения
Примечание
Схему Беренгера называют «раздельной», т.к. вводится разделение векторов. Есть также предложения
16.11.2016
Егоров Н.Н.
В третьих, теоретически, если выполняется условие
σi/ εo = σi*/ μo, (12)
то на
В третьих, теоретически, если выполняется условие
σi/ εo = σi*/ μo, (12)
то на
К сожалению, отражение все же есть:
- от первого слоя PML;
- между слоями PML, поскольку для экономии вычислительных ресурсов потери растут от слоя к слою (закон изменения потерь от первого слоя к последнему называется "профилем потерь");
- после последнего слоя PML, поскольку там находится PEC - граница.
Отражение от первого слоя PML и между слоями PML вызвано ошибками конечно-разностной дискретизации, и, в первую очередь, тем, что векторы E и H (а, следовательно, σi и σi*) не совпадают в пространстве. Для снижения отражения внутри PML необходимо ограничивать скорость роста потерь некоторым разумным пределом.
Беренгер показал, что при каждом конкретном значении σi и σi*, вследствие цифрового отражения, происходит отражение нормальных (перпендикулярных) к границе электромагнитных волн, частота которых ниже fc (частота отсечки). Чем больше σi и σi*, тем выше fc. Поэтому при проникновении волны в границы PML вначале идет отражение от первого слоя с проводимостью σ0 тех частот, которые ниже fc0. Затем идет отражение от полуслоя с проводимостью σ0* частот ниже fc0*. Причем fc0 < fc0*. И так далее – все более и более высокие частоты отражаются, причем отражение от σi и от σi*.
16.11.2016
Егоров Н.Н.
Отражение от PEC границы после последнего слоя PML происходит обычно уже
Отражение от PEC границы после последнего слоя PML происходит обычно уже
Для уменьшения отражения от первого слоя значения σ1 специально выбирается маленьким. Для уменьшения отражения между слоями профиль потерь выбирается с ограниченной скоростью роста потерь. Для уменьшения влияния волны, отраженной от РЕС границы, увеличивается количество слоев PML.
Как вариант, Беренгер предлагает следующий геометрический профиль потерь: (13)
где g - коэффициент геометрической прогрессии; Dx - шаг по пространству; с - скорость света; N - номер PML -слоя, считая от интерфейса счетного региона и границы; r - расстояние от границы; R(0) - коэффициент отражения от первого слоя. Рекомендуемое значениеR(0) =0.01 (1 %) Коэффициент прогрессии g рекомендуется брать 2,15. Это значение получено Беренгером экспериментально, хотя встречаются и другие рекомендации.
σ*(r) получается через σ(r) с использованием (12).
16.11.2016
Егоров Н.Н.
На низких частотах наблюдается резкое увеличение коэффициента отражения от границ PML.
На низких частотах наблюдается резкое увеличение коэффициента отражения от границ PML.
Для расчетов откликов от импульсов - ступенек с постоянной составляющей обратная величина к (14) - 1/fc - это максимальное время счета, при котором коэффициент отражения от первого слоя еще не превышает заданного R(0).
На самом деле, fc полученное по (14), - сильно заниженное значение. Реально происходят отражения и от других слоев, которые, как указано выше, имеют более высокую частоту отсечки. Практически fc нужно брать в 5 раз больше, что показывается в разделе сравнение граничных условий.
Существуют и другие профили потерь. В статье Беренгера "PML for the FDTD Solution of Wave-Structure Interaction Problems" (IEEE trans. on antennas and prop., vol. 44, N1, Jan 1996) есть богатый материал для размышлений о выборе профиля потерь и количества слоев в зависимости от разных факторов.
Дискретизация уравнений (11) происходит так же, как и дискретизация уравнений главного алгоритма. Для двух векторов Ex:
Exyn+1(i+1/2,j,k)=(1-0,5σy(r)∆t)/(1+0,5σy(r)∆t)Exyn(i+1/2,j,k)+∆t/(1+0,5 σy(r)∆t)/∆y (Hzxn+1/2(i+1/2,j+1/2,k)+Hzyn+1/2(i+1/2,j+1/2,k)-Hzxn+1/2(i+1/2,j-1/2,k)-Hzyn+1/2(i+1/2,j-1/2,k))
Exzn+1(i+1/2,j,k)=(1-0,5σz(r)∆t)/(1+0,5σz(r)∆t)Exzn(i+1/2,j,k)+∆t/(1+0,5σz(r)∆t)/∆z (Hyxn+1/2(i+1/2,j,k+1/2)+Hyzn+1/2(i+1/2,j,k+1/2)-Hyxn+1/2(i+1/2,j,k-1/2)-Hyzn+1/2(i+1/2,j,k-1/2))
Для остальных векторов Е все аналогично. Для векторов Н уравнения аналогичные, но вместо σ(r) берется σ*(r).
Hxyn+1/2(i,j+1/2,k+1/2)=(1-0,5σ*y(r)∆t)/(1+0,5σ*y(r)∆t)Hxyn-1/2(i,j+1/2,k+1/2)-∆t/(1+0,5σ*y(r)∆t)/∆y (Ezxn(i,j+1,k+1/2)+Ezyn(i,j+1,k+1/2)-Ezxn(i,j,k+1/2) - Ezyn (i,j,k+1/2))
16.11.2016
Егоров Н.Н.
Как можно заметить, в уравнениях для векторов Е имеются суммы типа
Как можно заметить, в уравнениях для векторов Е имеются суммы типа
На внешней границе PML, как уже говорилось, тангенциальные векторы Е равны нулю (PEC).
Профиль потерь зависит только от координаты, ведущей от интерфейса "счетный объем" - PML вглубь PML. На любой грани прямоугольного счетного объема вглубь PML ведет только одна координата. Допустим, это координата X. Тогда σx(r) меняется по заданному закону, а σy(r) = σz(r) = 0. На ребрах вглубь PML ведут уже две координаты (одна из σ равна нулю), а в углах вглубь ведут все три координаты и, следовательно, изменяются все σ (рис. 1).
16.11.2016
Егоров Н.Н.
Рис. 1. Счетный объем в окружении слоев PML
16.11.2016
Егоров Н.Н.
Рис. 1. Счетный объем в окружении слоев PML
16.11.2016
Егоров Н.Н.
Расширения базового алгоритма
Существует множество дополнений к базовому алгоритму, расширяющих его возможности.
Здесь
Расширения базового алгоритма
Существует множество дополнений к базовому алгоритму, расширяющих его возможности.
Здесь
· Сосредоточенные элементы
· Полярные диэлектрики
Сосредоточенные элементы
В англоязычной литературе сосредоточенные элементы называются «lumped», «lumped circuit elements».
М. Piket-May, A. Taflove J. Baron «FD-TD Modeling of Digital Signal Propagation in 3-D Circuits With Passive and Active Loads», IEEE Trans. On Mirowave Theory and Techiques, 1994, vol.42, №8. В этой статье приведены уравнения резисторов, резистивных источников напряжения, конденсаторов, индуктивностей, диодов и биполярных транзисторов. В более поздних публикациях разные авторы совершенствуют модели, добавляют новые элементы (например, источники тока, полевые транзисторы).
Отсутствие сосредоточенных элементов сужает область применения метода FDTD. При их отсутствии в расчетных задачах можно конструировать резисторы – ячейки с такой проводимостью, чтобы их сопротивление было равно заданному, и конденсаторы – три ячейки вдоль одной линии, две из них являются обкладками, а центральная – диэлектриком с такой диэлектрической проницаемостью, чтобы емкость была равна заданной. Я использовал этот способ задания элементов, пока не потребовалось задать индуктивность. Большую индуктивность задавать из ячеек практически нереально. И здесь на помощь пришла вышеупомянутая статья. Приведу кратко суть статьи с собственными поясненями.
Добавим к базовой формулировке вихревого уравнения для H (см. систему (1 из Введения)), которое выглядит как
(здесь и далее форма записи взята из статьи) ток сосредоточенного элемента JL:
Предположим, что элемент находится в свободном пространстве и ориентирован по оси Z. Тогда плотность тока выразится через ток IL как
IL может быть дифференциальной, интегральной, линейной или нелинейной функцией электрического потенциала
16.11.2016
Егоров Н.Н.
Проведя вывод дискретного уравнения так, как это было во введении, но
Проведя вывод дискретного уравнения так, как это было во введении, но
(1)
Примечание. Среднее слагаемое правой части уравнения – это общепринятый способ краткой записи разностной схемы. В данном случае это
(∆t/εo(Hxn+1/2(i,j+1/2,k+1/2)-Hxn+1/2(i,j-1/2,k+1/2))/∆y-(Hyn+1/2(i+1/2,j,k+1/2)-Hyn+1/2(i-1/2,j,k+1/2))/∆x),
Узнаете, что это за уравнение?
16.11.2016
Егоров Н.Н.
Резистор
Ток резистора находится по известной формуле I=U/R. В нашем случае U=∆Z(En+1+En)/2,
Резистор
Ток резистора находится по известной формуле I=U/R. В нашем случае U=∆Z(En+1+En)/2,
Подставляя последнее уравнение в (1) и решая его относительно En+1 получим:
16.11.2016
Егоров Н.Н.
Резистивный источник напряжения
Это резистор, в котором встроен источник напряжения. Очень удобный
Резистивный источник напряжения
Это резистор, в котором встроен источник напряжения. Очень удобный
Окончательно имеем:
16.11.2016
Егоров Н.Н.
Конденсатор
Ток через конденсатор I=C(du/dt), где С – емкость. Производную напряжения по
Конденсатор
Ток через конденсатор I=C(du/dt), где С – емкость. Производную напряжения по
Дальнейший вывод несложен, как и предыдущие:
16.11.2016
Егоров Н.Н.
Индуктивность
Ток в индуктивности, как известно, равен I=(1/L)∫Udt, т.е. величина интегральная. В
Индуктивность
Ток в индуктивности, как известно, равен I=(1/L)∫Udt, т.е. величина интегральная. В
Отсюда:
Полученная формула работает очень хорошо и устойчиво. Для вычислений требуется хранение суммы поля Ez в данной ячейке на всех предыдущих шагах, что не представляет особой трудности.
Итак, приведено описание четырех сосредоточенных элементов.
16.11.2016
Егоров Н.Н.
Расширения базового алгоритма
Полярные диэлектрики
Наличие относительной диэлектрической проницаемости большей, чем единица, связано
Расширения базового алгоритма
Полярные диэлектрики
Наличие относительной диэлектрической проницаемости большей, чем единица, связано
Здесь рассматриваются полярные диэлектрики. Их молекулы всегда представляют собой диполи. Например, молекулы воды. Наличие механически поворачивающихся молекул приводит к двум основным эффектам:
1. Всегда найдется такая высокая частота, при которой молекулы не успевают повернуться. Из-за этого диэлектрическая проницаемость падает.
2. При повороте молекул часть их кинетической энергии растрачивается на соударения с соседними молекулами и происходит нагрев диэлектрика. С ростом частоты этот процесс становится все заметнее и, наконец, говорят о том, что в диэлектрике растут потери на поляризацию или просто потери. С точки зрения стороннего наблюдателя, нагрев диэлектрика происходит так же, как и нагрев проводника с током за счет его сопротивления. Поэтому говорят о том, что в диэлектрике имеется некая проводимость потерь.
Непосвященный в эти процессы будет удивлен, но факт остается фактом: дистиллированная вода на частоте в пятьсот мегагерц ведет себя как проводник с проводимостью ~0,07 См/м, но постоянный ток вообще не проводит.
Белки, жиры, углеводы и вода, насыщенная солями, из которых состоят живые организмы, - все это в основном полярные диэлектрики. Сейчас модно рассчитывать эффекты взаимодействия электромагнитного излучения с человеком. Для синусоидального сигнала в расчетах можно брать электрофизические характеристики для конкретной частоты. Но для широкополосных воздействий учет частотных зависимостей совершенно необходим.
16.11.2016
Егоров Н.Н.
Применение метода FDTD для расчетов распространения нестационарных электромагнитных полей в условиях
Применение метода FDTD для расчетов распространения нестационарных электромагнитных полей в условиях
В классической постановке [1]
1. Yee, K. S.: "Numerical solution of initial boundary value problems involving Maxwell’s equations in isotropic media." IEEE Transactions on Antennas and Propagation, Vol. 14, No. 3, 1966, pp. 302-307.
метод FDTD требует, чтобы значения проводимости и диэлектрической проницаемости материалов не зависели от частоты. В то же время, для многих материалов эти параметры имеют ярко выраженную частотную зависимость. К частотно-зависимым материалам относятся, в частности биологические ткани и вода. Для биологических тканей в диапазоне частот от 40 МГц до 400 МГц диэлектрическая проницаемость уменьшается в 2-3 раза, проводимость возрастает также в 2-3 раза.
Если вычисления по методу FDTD проводятся для одной частоты, то можно задать значения проводимости и диэлектрической проницаемости материалов для конкретной частоты и этого обычно достаточно. Проблемы возникают тогда, когда необходимо провести численное моделирование взаимодействия широкополосных электромагнитных полей (например, коротких импульсов) c частотно-зависимыми материалами.
16.11.2016
Егоров Н.Н.
В [2]
2. Luebbers R. et all.: "A frequency-depended finite-difference time-domain
В [2]
2. Luebbers R. et all.: "A frequency-depended finite-difference time-domain
приводится частотно-зависимая формулировка (FD)2TD (frequency-dependent FDTD). Эта формулировка приведена с упрощениями, она не учитывает проводимости диэлектрика. Между тем, биологические ткани имеют существенную проводимость, которую следует учесть.
В [3]
3. Sullivan D. M.: "A frequency-depended FDTD Method for Biological Applications." IEEE Transactions on Microwave Theory and Techniques, Vol. 40, No. 3, March 1992, pp. 532-539.
приводит новую формулировку (FD2)TD, в которой учитывается проводимость биологических тканей. К несчастью, эта формулировка в статье приведена с ошибками.
В этой статье приводится новая формулировка (FD2)TD для биологических тканей с учетом проводимости среды.
16.11.2016
Егоров Н.Н.
Учет проводимости биологических тканей в уравнениях (FD2)TD
Предположим, что биологические материалы линейные
Учет проводимости биологических тканей в уравнениях (FD2)TD
Предположим, что биологические материалы линейные
Исходными являются уравнения Максвелла:
(1)
(2)
где - удельная проводимость. Из уравнения (2) следует выражение для H, которое полностью совпадает с классической формулировкой [1], поэтому вывод формулировки для Н не приводится.
16.11.2016
Егоров Н.Н.
Вектор смещения из (1) во временной области записывается как [2, 3]: (3)
где
Вектор смещения из (1) во временной области записывается как [2, 3]: (3)
где
Используя обозначения, принятые в [1], запишем и для каждой компоненты векторов и в (3) можем записать:
(4)
Уравнение (1) с учетом того, что x=iΔx, y=jΔy и z=kΔz, в конечно-разностной дискретной форме в декартовых координатах для компоненты Dy имеет вид:
(5)
16.11.2016
Егоров Н.Н.
Принимая во внимание, что D(t) и E(t) равны нулю при t<0,
Принимая во внимание, что D(t) и E(t) равны нулю при t<0,
(6)
Выражение для D на следующем временном шаге имеет вид:
(7)
Подставляя выражения (6) и (7) в (5) получаем:
(8)
16.11.2016
Егоров Н.Н.
Для удобства введем обозначения:
(9)
Подставляя (9) в (8) и решая последнее уравнение
Для удобства введем обозначения:
(9)
Подставляя (9) в (8) и решая последнее уравнение
получаем:
(10)
Выражения для компонент Ex и Ez могут быть получены аналогичным образом.
16.11.2016
Егоров Н.Н.
Если относительная диэлектрическая проницаемость не зависит от частоты, то χm=0 для
Если относительная диэлектрическая проницаемость не зависит от частоты, то χm=0 для
Для полярных диэлектриков комплексная диэлектрическая проницаемость определяется как (формула Дебая):
ε(ω)=ε∞+χ(ω), (11)
где - комплексная диэлектрическая проницаемость, зависящая от частоты.
Фурье- преобразование функции χ(ω) имеет вид:
где - функция, зависящая от шага по времени.
Так как t=mΔt имеем:
16.11.2016
Егоров Н.Н.
Выражение (10) содержит член
который, на первый взгляд, для вычисления требует хранения
Выражение (10) содержит член
который, на первый взгляд, для вычисления требует хранения
можем записать: (12)
Таким образом, для каждой компоненты
и требуется только одна переменная.
16.11.2016
Егоров Н.Н.
Вычисление электромагнитного поля в дальней зоне с использованием метода FDTD и
Вычисление электромагнитного поля в дальней зоне с использованием метода FDTD и
Решение уравнений Максвелла методом FDTD дает возможность вычисления электромагнитного внутри некоторой ограниченной области пространства. Ограничение на размеры области, в которой проводятся вычисления, накладываются ограниченной вычислительной мощность применяемых ЭВМ, а именно тремя факторами: ограниченной скоростью выполнения арифметических операций, ограниченным объемом оперативной памяти и ограниченной скоростью обмена данными между процессором с памятью. Однако существует ряд задач, в которых необходим расчет электромагнитных полей на больших расстояниях от некоторого объекта, излучающего или рассевающего электромагнитное поле. В их числе расчет диаграммы направленности антенн, определение эффективной поверхности отражения радиолокационных целей, решение задач дифракции и. др. В данной статье предложен эффективный способ вычисления полей в дальней зоне с использованием результатов вычислений ближнего поля методом FDTD. Под дальним полем здесь подразумевается поле за пределами вычислительного объема, а ближнее поле здесь означает поле, вычисляемое непосредственно методом FDTD.
16.11.2016
Егоров Н.Н.
Преобразование ближнего поля в дальнее поле с использованием интеграла Кирхгофа
Существует несколько
Преобразование ближнего поля в дальнее поле с использованием интеграла Кирхгофа
Существует несколько
Например, для вычисления дальнего поля часто применяют интегрирование по элементам, на которые разбивается замкнутая поверхность. Каждый элемент поверхности является элементарным электрическим и магнитным диполем одновременно. Если при этом ближнее поле вычисляется методом FDTD, выполняются следующие операции:
- Все компоненты E и H полей на поверхности, по которой происходит интегрирование, сохраняются на всех шагах по времени. Т.е. к концу вычислений известна временная форма поля в каждой точке поверхности.
- Поля Е и Н на поверхности преобразуются в эквивалентные электрические и магнитные токи.
- Осуществляется преобразование токов на поверхности в частотную область.
- Вычисляются составляющие Eφ, Eθ, Еr и соответствующие составляющие поля H для всех частот и от каждого элемента Гюйгенса на поверхности. Элементы Гюйгенса в данном случае - это обычные ячейки FDTD. С каждой ячейкой сопоставляется своя собственная система полярных координат.
- Значения Eφ, Eθ, Еr и соответствующих составляющих поля H преобразуются в прямоугольные координаты (Ex, Ey, Еz) и складываются в точке наблюдения. Только на этом шаге происходит собственно интегрирование по поверхности.
- Производится обратное преобразование во временную область. Теперь искомое поле найдено.
16.11.2016
Егоров Н.Н.
16.11.2016
Егоров Н.Н.
16.11.2016
Егоров Н.Н.
Формула (6) выражает принцип Гюйгенса, согласно которому каждая точка на волновом
Формула (6) выражает принцип Гюйгенса, согласно которому каждая точка на волновом
Шаг по времени при вычислении интеграла (6) тесно связан с шагом по времени FDTD и равен ему. Выходная последовательность в точке наблюдения имеет такой же шаг по времени. Однако задержка R/c может не быть кратной шагу по времени. Поэтому получаемое время задержки округляется до ближайшего времени, кратного шагу FDTD. Возникающая при этом ошибка незначительна, т.к. шаг по времени в классическом FDTD мал по сравнению с периодом колебаний вычисляемого сигнала.
Чтобы привести (6) к виду, удобному для применения совместно с алгоритмом FDTD, необходимо выполнить ряд преобразований. При этом учтем, что все участки поверхности, по которой ведется интегрирование, в алгоритме FDTD всегда перпендикулярны одной из осей координат.
16.11.2016
Егоров Н.Н.
Предположим, что поверхность dA’ перпендикулярна оси Z и точка p’ имеет
Предположим, что поверхность dA’ перпендикулярна оси Z и точка p’ имеет
Применим преобразование подынтегральных слагаемых (6) в конечные разности:
(7)
16.11.2016
Егоров Н.Н.
Для упрощения дальнейшей записи введем обозначения:
(8)
и применим стандартные для FDTD обозначения
Для упрощения дальнейшей записи введем обозначения:
(8)
и применим стандартные для FDTD обозначения
Выражение (6) с учетом введенных обозначений (7) и (8) запишется для одной площадки dA в виде:
(9)
где Ai,j = Δx Δy.
В (9) нет смысла записывать сумму по всем участкам поверхности, т.к. tn* в общем случае для каждого участка имеет разное значение, поскольку различно расстояние R. Временная последовательность в точке наблюдения p получается суммированием значений Ψ(p, tn*)на каждом шаге по времени по всем участкам поверхности с учетом времени запаздывания.
16.11.2016
Егоров Н.Н.
В (9) значение функции Ψ(p, tn*) вычисляется с использованием значений Ψ
В (9) значение функции Ψ(p, tn*) вычисляется с использованием значений Ψ
(10)
где
(11)
В приведенных (10) и (11) присутствуют шаги n, n+1 и n+2. Но можно вычислить все эти значения для одного шага, допустим, текущего шага n+1. Тогда, если считать, что вычисляется шаг n+1, то F1(n) добавляется в tn+1*- шаг вычисляемой последовательности. F2(n+1) добавляется в предыдущую временную точку tn*, а F3(n+2) добавляется еще на один временной шаг раньше в вычисляемой последовательности (tn-1*). В этом случае F3(n+2) = ‑ F1(n).
16.11.2016
Егоров Н.Н.