Сайт Александра Зеленина

Сайт Александра Зеленина

Введение в метод FDTD

Александр Зеленин.

Введение в метод 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 = εεoEx /∂t + σEx;

 

rot(E) Y = - μμoHy /∂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 (скорости света в вакууме).

файл:/zfdtd/method/intrfdtd.htm