Александр Зеленин.
Введение в метод FDTD.
Аббревиатура FDTD расшифровывается как "finite-difference time-domain", а в русскоязычной литературе иногда выглядит как КРВО - "конечные разности во временной области", что является переводом с английского. В принципе этот метод - понятие чисто математическое и обозначает один из многочисленных методов решения дифференциальных уравнений, но среди тех, кто занимается решением задач электротехники, аббревиатура FDTD в настоящее время является синонимом решения вихревых дифференциальных уравнений Максвелла.
В 1966 г. Йе (Yee) разработал технику, реализующую явную конечно - разностную схему второго порядка для решения вихревых уравнений Максвелла в пространстве и времени.
Исходными являются уравнения Максвелла в дифференциальной форме:
|
rot(H) = ∂D/∂t + J; rot(E) = - ∂B/∂t; |
(1) |
а также
|
D = ε εo E; J = σ E; B = μ μo H; |
(2) |
Здесь E - вектор напряженности электрического поля (В/м), Н - вектор напряженности магнитного поля (А/м), ε, μ - относительные диэлектрическая и магнитная проницаемости (без размерности), εo - диэлектрическая постоянная (Ф/м), μo - магнитная постоянная (Гн/м), B - вектор магнитной индукции (Тл), D - вектор электрического смещения (Кл/м2), J - вектор плотности тока (А/м2), σ - электрическая проводимость (См/м), и t - время в секундах.
εo = 107/(4pc2), где с - скорость света в вакууме (2,997925010 x 108 м/с).
μo = 4p/107.
Оба уравнения (1) содержат пространственные и временные производные.
Для решения уравнения (1) следует выразить в декартовых координатах векторы Е и Н:
|
Е = 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; |
(3) |
где Ex, Ey, Ez, Hx, Hy, Hz - проекции векторов на координатные оси, а X, Y, Z - единичные векторы.
Остальные величины в (1) - D, B, J - выразим через E и H. Величины E и H для нас будут основными.
Примечание: существуют и другие подходы, когда в уравнениях (1) вначале оставляют D и/или B, но в конце концов всё равно выражаются вектора Е и Н . Также следует указать, что уравнения (1) записаны не полностью. Например, в них не учитываются сторонние токи.
Yee (1966) предложил пространственную сетку для конечно-разностной аппроксимации, в которую поместил вектора Ex, Ey, Ez, Hx, Hy, Hz. Фрагмент сетки Yee показан на (рис.1).
|
Рис. 1. |
Все компоненты (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 - шаг по времени.
Поля E и H вычисляются со сдвигом на полшага по времени. Обозначения, введенные Yee, следующие: En - значение поля E на только что вычисленном шаге; En+1 - значение поля E на вычисляемом сейчас шаге по времени. Hn-1/2 - значение поля H на только что вычисленном шаге; Hn+1/2 - значение поля на вычисляемом сейчас полушаге по времени. Из этих обозначений следует, что процедура вычислений начинается с поля Hn+1/2, потому что в момент t=0 (n=0) установлены начальные условия по всему счетному объему: все значения полей E и H равны нулю. Хотя в принципе это лишь наиболее распространенная условность. Можно считать, что пространственная сетка проходит через вектора H, что процедура счета начинается с поля E.
Теперь, когда введены основные обозначения, покажем вывод выражений, пригодных для расчетов с помощью компьютера и которым уже 40 лет.
Поставим (3) и (2) в (1). Получим:
|
rot(H) X = εεo∂Ex /∂t + σEx;
rot(E) Y = - μμo∂Hy /∂t; |
(5) |
Применяя конечно-разностную аппроксимацию, преобразуем (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) - Hy n-1/2(i+1/2,j,k+1/2))/∆t;
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; |
(6) |
Подставляя (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);
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);
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); |
(8) |
Аналогичные выражения можно получить для остальных четырех компонент ячейки Yee.
Из выражений (7) и (8) видно, что значения μ, ε и σ задаются для каждого их векторов ячейки и могут быть различными в разных направлениях. Т.е. при необходимости можно задать анизотропию материалов для Е и/или Н полей.
Выражения (7) и (8) являются достаточными для многих решаемых задач, но для расчетов сосредоточенных элементов (источников напряжения, индуктивностей, транзисторов и т.п.), а также для расчетов материалов с нелинейными свойствами требуется их модификация.
В заключение следует упомянуть, что явные конечно-разностные схемы требуют специальных условий для устойчивой работы. Для метода FDTD это условие имеет вид:
∆t ≤ 1/(v √((1/∆x2)+ (1/∆y2)+ (1/∆z2))),
где v - максимальная скорость электромагнитных волн в счетном объеме, а выражение (1/∆x2)+ (1/∆y2)+ (1/∆z2) находится под знаком квадратного корня.
Обычно v = c (скорости света в вакууме).
