🗊Презентация Алгоритм FDTD. Введение в метод FDTD

Алгоритм FTDT СМвТФ Егоров Н.Н.

Введение в метод FDTD FDTD - "finite-difference time-domain", в русскоязычной литературе КРВО - "конечные разности во временной области" метод FDTD - один из многочисленных методов решения дифференциальных уравнений, но среди тех, кто занимается решением задач электродинамики, FDTD является синонимом решения вихревых дифференциальных уравнений Максвелла.

В 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.

Исходными являются уравнения Максвелла в дифференциальной форме rot(H) = ∂D/∂t + J; div(B)=0; (1) rot(E) = - ∂B/∂t; div(D)=ρ; D = ε εo E; J = σ E; (2) B = μ μo H Два уравнения (1) содержат пространственные и временные производные.

Для решения уравнения (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) записаны не полностью. Например, в них не учитываются сторонние токи.

Yee (1966) предложил пространственную сетку для конечно-разностной аппроксимации, в которую поместил вектора Ex, Ey, Ez, Hx, Hy, Hz. Фрагмент сетки Yee показан на (рис.1).

Все компоненты (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 - шаг по времени.

Поля E и H вычисляются со сдвигом на полшага по времени. Обозначения, введенные Yee, следующие: Поля 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. Теперь, когда введены основные обозначения, покажем вывод выражений, пригодных для расчетов с помощью компьютера и которым уже 50 лет.

Поставим (3) и (2) в (1). Получим: rot(H) X = εεo∂Ex /∂t + σEx; (5) 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;

Подставляя (6) в (5) и решая получившиеся выражения относительно Hyn+1/2(i+1/2,j,k+1/2) и Exn+1(i+1/2,j,k) получим: Подставляя (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); (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) видно, что значения μ, ε и σ задаются для каждого их векторов ячейки и могут быть различными в разных направлениях. Т.е. при необходимости можно задать анизотропию материалов для Е и/или Н полей.

Выражения (7) и (8) являются достаточными для многих решаемых задач, но для расчетов сосредоточенных элементов (источников напряжения, индуктивностей, транзисторов и т.п.), а также для расчетов материалов с нелинейными свойствами требуется их модификация. Выражения (7) и (8) являются достаточными для многих решаемых задач, но для расчетов сосредоточенных элементов (источников напряжения, индуктивностей, транзисторов и т.п.), а также для расчетов материалов с нелинейными свойствами требуется их модификация. Следует упомянуть, что явные конечно-разностные схемы требуют специальных условий для устойчивой работы. Для метода FDTD это условие имеет вид: где v - максимальная скорость электромагнитных волн в счетном объеме. Обычно v = c (скорости света в вакууме).

Базовый алгоритм FDTD Выше кратко описан вывод конечно-разностных уравнений в декартовых координатах для двух компонент сетки Yee из шести. На рис. 2. приводится полная система для всех векторов сетки Yeе

Примечания. 1. Полуцелые индексы, которые применены в нотации Yee, часто заменяются на целые путем уменьшения на пол-индекса, например, (j+1/2) заменено на (j), а (j-1/2) заменено на (j-1). Это сделано для удобства программирования, т.к. не бывает полуцелых индексов массивов. 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).

Шесть уравнений базового алгоритма очень просты в реализации. Для этого необходимо создать шесть массивов Ex, Ey, Ez, Hx, Hy и Hz. А также массивы электрофизических характеристик ε σ μ, которые назовем Eps, Sig и Mu. Этих массивов три штуки, если электрофизические характеристики задаются для ячейки Yee в целом. Здесь возможны разные варианты. Если ε σ и μ планируется задавать для каждого вектора свои, то потребуются массивы Epsx, Epsy, Epsz, Sigx и т.д. Если не планируется применение магнитных материалов, то массив Mu можно вообще не создавать. Если планируется проводить расчеты только для идеальных проводников (PEC) в свободном пространстве, то массивы электрофизических констант не нужны вовсе. В этом случае векторы Е, которые принадлежат PEC-объекту, попросту обнуляются, а все остальные считаются принадлежащими свободному пространству. Шесть уравнений базового алгоритма очень просты в реализации. Для этого необходимо создать шесть массивов Ex, Ey, Ez, Hx, Hy и Hz. А также массивы электрофизических характеристик ε σ μ, которые назовем Eps, Sig и Mu. Этих массивов три штуки, если электрофизические характеристики задаются для ячейки Yee в целом. Здесь возможны разные варианты. Если ε σ и μ планируется задавать для каждого вектора свои, то потребуются массивы Epsx, Epsy, Epsz, Sigx и т.д. Если не планируется применение магнитных материалов, то массив Mu можно вообще не создавать. Если планируется проводить расчеты только для идеальных проводников (PEC) в свободном пространстве, то массивы электрофизических констант не нужны вовсе. В этом случае векторы Е, которые принадлежат PEC-объекту, попросту обнуляются, а все остальные считаются принадлежащими свободному пространству. Чем больше возможностей мы хотим реализовать в базовом алгоритме, тем больший объем памяти требуется для переменных и, соответственно, больше затраты времени при расчетах. Раньше значительное ускорение вычислений можно было получить создав еще ряд массивов и поместив в них предварительно вычисленные коэффициенты С, С1 и С2 из уравнений (7) и (8). Но сейчас в гигагерцовых процессорах скорость вычислений арифметики с плавающей точкой выросла больше, чем скорость выборки дополнительных элементов массива из памяти, поэтому смысл в предварительном вычислении этих коэффициентов снизился. В программе ZFDTD пожертвовали скоростью во имя экономии памяти, хотя в ранних версиях программы было наоборот.

Вычисление полей Е и Н ведется рекурсивно: для вычисления значений на текущем шаге по времени используются значения на предыдущем. При этом результат сразу помещается в прежний массив "затирая" старые значения. Но в некоторых случаях все-таки требуется иметь значения полей на предыдущем шаге по времени одновременно с текущим шагом. Например, для того, чтобы вычислить значение поля Е на полушаге по времени как среднее арифметическое значения на двух соседних шагах. В этом случае волей-неволей приходится использовать еще несколько массивов. При вычислении энергии и других величин, производных от полей Е и Н, также потребуются дополнительные массивы для хранения результатов. Вычисление полей Е и Н ведется рекурсивно: для вычисления значений на текущем шаге по времени используются значения на предыдущем. При этом результат сразу помещается в прежний массив "затирая" старые значения. Но в некоторых случаях все-таки требуется иметь значения полей на предыдущем шаге по времени одновременно с текущим шагом. Например, для того, чтобы вычислить значение поля Е на полушаге по времени как среднее арифметическое значения на двух соседних шагах. В этом случае волей-неволей приходится использовать еще несколько массивов. При вычислении энергии и других величин, производных от полей Е и Н, также потребуются дополнительные массивы для хранения результатов. В общем, по отношению к оперативной памяти алгоритм FDTD довольно прожорлив, но для больших задач этот алгоритм, по литературным данным, требует меньше памяти, чем интегральные методы. Ход вычислений рассмотрим на примере, написанном на языке Pascal

Две процедуры вычисляют поля E и H: // ОПРЕДЕЛЕНИЕ МАГНИТНОГО ПОЛЯ H Procedure TimeStepForH; 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

// ОПРЕДЕЛЕНИЕ ЭЛЕКТРИЧЕСКОГО ПОЛЯ E // ОПРЕДЕЛЕНИЕ ЭЛЕКТРИЧЕСКОГО ПОЛЯ E Procedure TimeStepForE; var I,J,K: Integer; begin for I:=1 to XSize-1 do for J:=2 to YSize-1 do for K:=2 to ZSize-1 do 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

Массивы AE[I,J,K], BE[I,J,K] и AН[I,J,K] - это предварительно вычисленные коэффициенты С1, С2 и С соответственно из (7) и (8). XSize, YSize и ZSize - размеры счетного пространства. А коэффициенты TX, TY и TZ равны 1/∆X, 1/∆Y и 1/∆Z соответственно. Как видно, здесь применены предварительно вычисленные массивы коэффициентов С1, С2 и С. Массивы AE[I,J,K], BE[I,J,K] и AН[I,J,K] - это предварительно вычисленные коэффициенты С1, С2 и С соответственно из (7) и (8). XSize, YSize и ZSize - размеры счетного пространства. А коэффициенты TX, TY и TZ равны 1/∆X, 1/∆Y и 1/∆Z соответственно. Как видно, здесь применены предварительно вычисленные массивы коэффициентов С1, С2 и С. Приведенный код вычисления полей рассчитан на такое представление объектов, когда считается, что электрофизические характеристики пространства в ячейке (i,j,k) задаются одинаковыми для всех трех Нx,y,z(i,j,k) и трех Ex,y,z(i,j,k) векторов. При таком представлении точность вычислений на границах раздела сред снижается до первого порядка. Обычно это проявляется в снижении на 25..40% вычисленной резонансной частоты, такого же увеличения вычисленных токов в длинных структурах, снижении в разы скорости затухания свободных колебаний и в других ошибках.

Приведенные процедуры вызываются в цикле по времени: // ГЛАВНЫЙ ЦИКЛ ПО ВРЕМЕНИ T:=0; //Исходное значение времени for I:=0 to <заданное количество шагов по времени> do begin TimeStepForH; //Вычисление H по всему объему BordersSymmetryH; //Выполнение граничных условий для симметрии по Н TimeStepForE; //Вычисление E по всему объему BordersABC; // Граничные УСЛОВИЯ ПОГЛОЩЕНИЯ BordersSymmetryE;//Выполнение граничных условий для симметрии по Е BordersРЕС;//Выполнение граничных условий РЕС T:=T+dT; // инкремент времени на один шаг end; //ГЛАВНЫЙ ЦИКЛ ПО ВРЕМЕНИ это самое простое место в программе

Граничные условия для FDTD В счетном объеме каждый вектор Е или Н вычисляется через 4 соседних вектора. Так происходит по всему объему. Но на границах самые последние векторы Е имеют: на гранях параллелепипеда счетного объема только три соседних вектора Н из четырех необходимых; на ребрах - два. Поэтому точно вычислить поле Е на границах невозможно. Проблема вычисления граничных полей решается различными способами. 1. Условия PEC – (photo electrochemical cell) идеальный проводник. Условия PEC таковы, что граничные вектора Е никогда не вычисляются, а, значит, всегда равны нулю. Как известно, поле Е всегда равно нулю в идеальном проводнике, поэтому такие границы ведут себя как идеальный проводник: электромагнитные волны 100 % отражаются обратно в счетный объем.

2. Условия симметрии. 2. Условия симметрии. В некоторых случаях поле Е или поле Н может быть симметричным относительно некоторой плоскости. Тогда в этой плоскости можно задать условие симметрии и тем самым вдвое уменьшить счетный объем. При этом рядом с данной плоскостью симметрии будет проходить граница счетного объема с условиями симметрии. Симметрия может быть четной или нечетной. При нечетной симметрии плоскость симметрии проходит внутри счетного объема параллельно грани на расстоянии пол-ячейки от грани. Условия нечетной симметрии для симметрии по Е получаются простым переносом значений ближайших к границе векторов Е на саму границу, а для симметрии по Н получаются так же, но при этом вектор Е меняет свой знак. Ниже приведены примеры нечетной симметрии для Е и Н.

Пример: нечетная симметрия по Е: // ПЛОСКОСТЬ Y=1/2 симметрична по E 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]; Пример: нечетная симметрия по Н: // ПЛОСКОСТЬ 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];

Условия четной симметрии несколько сложнее. Плоскость симметрии проходит на расстоянии целой ячейки от границы, поэтому кроме поля Е на границе необходимо помнить о прилегающих к границе векторах Н. При этом симметрии по Е и Н отличаются друг от друга только знаками переносимых значений. Пример: четная симметрия по Н: // ПЛОСКОСТЬ 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; При четной симметрии граничные поля Е и Н устанавливаются не одновременно, а на соответствующих полушагах по времени.

Простые условия поглощения (АВС) Для условий поглощения значения векторов электрического поля на границе вычисляются на основании известных полей в приграничных слоях. Причем берутся приграничные поля не только на текущем шаге по времени, но и на предыдущих шагах. Все эти условия пытаются спрогнозировать поле на границе. В литературе встречается описание целого ряда простых условий поглощения разных авторов. Всего их наберется с десяток. Практически чаще всего применяют условия Мура (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 - временной индекс). Формула приведена в двух видах, дополняющих друг друга в пространстве и времени. На рисунках пустой кружок - вычисляемое значение, черные кружки - требуемые значения. Сразу можно определить, какие переменные приграничной области нужно хранить в дополнительных массивах и сколько шагов по времени они должны храниться. Итак, граничных условий много. Здесь рассмотрим три варианта простейших граничных условий: Мура 1-го порядка, Лиао 3-го порядка и RT. В таблице S (большое) и s (малое) - разные переменные! S = (dT * c)/ D - число Куранта, где D - шаг по пространству, dT- шаг по времени, с- скорость света, а s - координата границы. s-1 - одна ячейка внутрь счетного объема, s-2 - две ячейки от границы и т.д. На поясняющих рисунках от границы (surface - внешней поверхности счетного объема) по горизонтали отложена координата s, а по вертикали - время в шагах счета (n - временной индекс). Формула приведена в двух видах, дополняющих друг друга в пространстве и времени. На рисунках пустой кружок - вычисляемое значение, черные кружки - требуемые значения. Сразу можно определить, какие переменные приграничной области нужно хранить в дополнительных массивах и сколько шагов по времени они должны храниться.

Видно, что больше всего массивов переменных нужно хранить для условий RT (5), меньше всего - для условий Мура (1). Для условий Лиао нужно 3 дополнительных массива переменных. Видно, что больше всего массивов переменных нужно хранить для условий RT (5), меньше всего - для условий Мура (1). Для условий Лиао нужно 3 дополнительных массива переменных. А теперь пример реализации условий Мура 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;

//ВЕРХНЯЯ И НИЖНЯЯ ПЛОСКОСТИ //ВЕРХНЯЯ И НИЖНЯЯ ПЛОСКОСТИ // 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

// Ex для двух плоскостей XxY [ближней и дальней] // Ex для двух плоскостей XxY [ближней и дальней] for I:=1 to NX-1 do for J:=2 to NY-1 do begin 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;

// сохранение приграничных полей для MursFirstABC // сохранение приграничных полей для 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-го порядка.

4. Условия PМL - идеально сочетающиеся слои "Perfectly Matched Layer" выглядит как "идеально согласованный (сочетающийся) слой“ Условия PML обладают низким коэффициентом отражения (по некоторым данным в миллион раз меньше, чем у условий Мура), а также практической независимостью от угла падения волны. К недостаткам условий PML следует отнести значительно больший объем требуемой памяти, чем для условий ABC и наличие нижней граничной частоты, для снижения которой требуется увеличения количества слоев PML, а, следовательно, требуемой памяти. Как следствие увеличения требуемого объема памяти происходить снижение скорости вычислений. За пределами главного счетного объема добавляются дополнительные ячейки, окружающие счетный объем по периметру несколькими слоями. Электрическое и магнитное поле в этих ячейках вычисляется почти так же, как и в счетном объеме, но есть отличия.

Во-первых, в уравнения в обязательном порядке вводятся электрические и магнитные потери. В главном алгоритме этих потерь может не быть вообще. В уравнениях главного алгоритма (рис. 2) учтены только электрические потери. Магнитные потери вводятся так же, как и электрические, путем задания "плотности магнитных токов". Тогда уравнения Максвелла (1) и (2) выглядят так: Во-первых, в уравнения в обязательном порядке вводятся электрические и магнитные потери. В главном алгоритме этих потерь может не быть вообще. В уравнениях главного алгоритма (рис. 2) учтены только электрические потери. Магнитные потери вводятся так же, как и электрические, путем задания "плотности магнитных токов". Тогда уравнения Максвелла (1) и (2) выглядят так: 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)

(11)

Примечание Схему Беренгера называют «раздельной», т.к. вводится разделение векторов. Есть также предложения других авторов применять «нераздельную» схему. Если кто-нибудь проверял работу «нераздельной» схемы, поделитесь впечатлениями. Я не доверяю автору «нераздельной» схемы из-за целого ряда некорректных статей за его авторством, а также авторством его учеников и соратников. Сложилось впечатление, что некий клан печет статьи как пирожки не придавая внимания содержанию – лишь бы побольше опубликовать. Буду рад, если мое впечатление ошибочно.

В третьих, теоретически, если выполняется условие В третьих, теоретически, если выполняется условие σi/ εo = σi*/ μo, (12) то на границе раздела двух сред скорость электромагнитных волн не изменяется и отражение рано нулю. В то же время, поскольку σi и σi* не равны нулю, то происходит поглощение электромагнитных волн в недрах PML. К сожалению, отражение все же есть: - от первого слоя 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*.

Отражение от PEC границы после последнего слоя PML происходит обычно уже для весьма ослабленной волны. Отраженная волна на обратном пути продолжает ослабляться. Но если слоев мало (обычно < 5), то отраженная волна может быть существенной. Отражение от PEC границы после последнего слоя PML происходит обычно уже для весьма ослабленной волны. Отраженная волна на обратном пути продолжает ослабляться. Но если слоев мало (обычно < 5), то отраженная волна может быть существенной. Для уменьшения отражения от первого слоя значения σ1 специально выбирается маленьким. Для уменьшения отражения между слоями профиль потерь выбирается с ограниченной скоростью роста потерь. Для уменьшения влияния волны, отраженной от РЕС границы, увеличивается количество слоев PML. Как вариант, Беренгер предлагает следующий геометрический профиль потерь: (13) где g - коэффициент геометрической прогрессии; Dx - шаг по пространству; с - скорость света; N - номер PML -слоя, считая от интерфейса счетного региона и границы; r - расстояние от границы; R(0) - коэффициент отражения от первого слоя. Рекомендуемое значениеR(0) =0.01 (1 %) Коэффициент прогрессии g рекомендуется брать 2,15. Это значение получено Беренгером экспериментально, хотя встречаются и другие рекомендации. σ*(r) получается через σ(r) с использованием (12).

На низких частотах наблюдается резкое увеличение коэффициента отражения от границ PML. Нижняя граничная частота отсечки для известного значения электрической проводимости на границе находится из выражения: (14) На низких частотах наблюдается резкое увеличение коэффициента отражения от границ PML. Нижняя граничная частота отсечки для известного значения электрической проводимости на границе находится из выражения: (14) Для расчетов откликов от импульсов - ступенек с постоянной составляющей обратная величина к (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))

Как можно заметить, в уравнениях для векторов Е имеются суммы типа (Hzxn+1/2(i+1/2,j+1/2,k) + Hzyn+1/2 (i+1/2,j+1/2,k)). Это и есть полный вектор, в данном случае Hz. Когда вычисляется поле на границе счетного объема, то из счетного объема берется полный векторHz, а из граничного PML-слоя – сумма Hzx+ Hzy. Также и для других векторов. Как можно заметить, в уравнениях для векторов Е имеются суммы типа (Hzxn+1/2(i+1/2,j+1/2,k) + Hzyn+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 ведут уже две координаты (одна из σ равна нулю), а в углах вглубь ведут все три координаты и, следовательно, изменяются все σ (рис. 1).

Рис. 1. Счетный объем в окружении слоев PML

Расширения базового алгоритма Существует множество дополнений к базовому алгоритму, расширяющих его возможности. Здесь представлены два вида расширений: · Сосредоточенные элементы · Полярные диэлектрики Сосредоточенные элементы В англоязычной литературе сосредоточенные элементы называются «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 может быть дифференциальной, интегральной, линейной или нелинейной функцией электрического потенциала

Проведя вывод дискретного уравнения так, как это было во введении, но с учетом тока через сосредоточенный элемент, получим: Проведя вывод дискретного уравнения так, как это было во введении, но с учетом тока через сосредоточенный элемент, получим: (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), Узнаете, что это за уравнение?

Резистор Ток резистора находится по известной формуле I=U/R. В нашем случае U=∆Z(En+1+En)/2, т.е. берется среднее арифметическое напряженности Е поля на соседних шагах по времени, чтобы получить напряженность на шаге (n+1/2). Это приближение необходимо для придания устойчивости вычислительной схемы. Имеем: Подставляя последнее уравнение в (1) и решая его относительно En+1 получим:

Резистивный источник напряжения Это резистор, в котором встроен источник напряжения. Очень удобный и абсолютно физичный элемент. Позволяет моделировать генератор с заданным выходным сопротивлением. Отличается от резистора только наличием дополнительного тока, равного напряжению источника Vs, деленному на сопротивление: Окончательно имеем:

Конденсатор Ток через конденсатор I=C(du/dt), где С – емкость. Производную напряжения по времени выразим через разность напряженности поля на соседних шагах по времени: Дальнейший вывод несложен, как и предыдущие:

Индуктивность Ток в индуктивности, как известно, равен I=(1/L)∫Udt, т.е. величина интегральная. В дискретной форме ток запишем как: Отсюда: Полученная формула работает очень хорошо и устойчиво. Для вычислений требуется хранение суммы поля Ez в данной ячейке на всех предыдущих шагах, что не представляет особой трудности. Итак, приведено описание четырех сосредоточенных элементов.

Расширения базового алгоритма Полярные диэлектрики Наличие относительной диэлектрической проницаемости большей, чем единица, связано с поляризацией молекул диэлектрика. При приложении внешнего электрического поля к диэлектрику молекулы полярного диэлектрика поворачиваются, а неполярного - деформируются. Здесь рассматриваются полярные диэлектрики. Их молекулы всегда представляют собой диполи. Например, молекулы воды. Наличие механически поворачивающихся молекул приводит к двум основным эффектам: 1. Всегда найдется такая высокая частота, при которой молекулы не успевают повернуться. Из-за этого диэлектрическая проницаемость падает. 2. При повороте молекул часть их кинетической энергии растрачивается на соударения с соседними молекулами и происходит нагрев диэлектрика. С ростом частоты этот процесс становится все заметнее и, наконец, говорят о том, что в диэлектрике растут потери на поляризацию или просто потери. С точки зрения стороннего наблюдателя, нагрев диэлектрика происходит так же, как и нагрев проводника с током за счет его сопротивления. Поэтому говорят о том, что в диэлектрике имеется некая проводимость потерь. Непосвященный в эти процессы будет удивлен, но факт остается фактом: дистиллированная вода на частоте в пятьсот мегагерц ведет себя как проводник с проводимостью ~0,07 См/м, но постоянный ток вообще не проводит. Белки, жиры, углеводы и вода, насыщенная солями, из которых состоят живые организмы, - все это в основном полярные диэлектрики. Сейчас модно рассчитывать эффекты взаимодействия электромагнитного излучения с человеком. Для синусоидального сигнала в расчетах можно брать электрофизические характеристики для конкретной частоты. Но для широкополосных воздействий учет частотных зависимостей совершенно необходим.

Применение метода 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 частотно-зависимыми материалами.

В [2] В [2] 2. Luebbers R. et all.: "A frequency-depended finite-difference time-domain formulation for dispersive materials." IEEE Transactions on Electromagnetic Compatibility, Vol. 32, No. 3, Aug. 1990, pp. 222-227. приводится частотно-зависимая формулировка (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 для биологических тканей с учетом проводимости среды.

Учет проводимости биологических тканей в уравнениях (FD2)TD Предположим, что биологические материалы линейные и изотропные, магнитная проницаемость является константой, не зависящей от частоты. Для этого случая в [2] приводится вывод уравнений без учета проводимости. Добавляя проводимость (это J в (1), см. ниже), повторим вывод уравнений. Исходными являются уравнения Максвелла: (1) (2) где - удельная проводимость. Из уравнения (2) следует выражение для H, которое полностью совпадает с классической формулировкой [1], поэтому вывод формулировки для Н не приводится.

Вектор смещения из (1) во временной области записывается как [2, 3]: (3) Вектор смещения из (1) во временной области записывается как [2, 3]: (3) где - диэлектрическая постоянная, - относительная диэлектрическая проницаемость на "бесконечной" частоте, - электрическая восприимчивость. Используя обозначения, принятые в [1], запишем и для каждой компоненты векторов и в (3) можем записать: (4) Уравнение (1) с учетом того, что x=ix, y=jy и z=kz, в конечно-разностной дискретной форме в декартовых координатах для компоненты Dy имеет вид: (5)

Принимая во внимание, что D(t) и E(t) равны нулю при t<0, можем записать (4) в виде: Принимая во внимание, что D(t) и E(t) равны нулю при t<0, можем записать (4) в виде: (6) Выражение для D на следующем временном шаге имеет вид: (7) Подставляя выражения (6) и (7) в (5) получаем: (8)

Для удобства введем обозначения: Для удобства введем обозначения: (9) Подставляя (9) в (8) и решая последнее уравнение относительно получаем: (10) Выражения для компонент Ex и Ez могут быть получены аналогичным образом.

Если относительная диэлектрическая проницаемость не зависит от частоты, то m=0 для всех m и 0=0. В этом случае выражение (10) совпадает с [1]. Если относительная диэлектрическая проницаемость не зависит от частоты, то m=0 для всех m и 0=0. В этом случае выражение (10) совпадает с [1]. Для полярных диэлектриков комплексная диэлектрическая проницаемость определяется как (формула Дебая): ()=+(), (11) где - комплексная диэлектрическая проницаемость, зависящая от частоты. Фурье- преобразование функции () имеет вид: где - функция, зависящая от шага по времени. Так как t=mt имеем:

Выражение (10) содержит член Выражение (10) содержит член который, на первый взгляд, для вычисления требует хранения большого количества переменных. Но, поскольку функция () носит экспоненциальный характер, то, как показано в [2], определяя, что можем записать: (12) Таким образом, для каждой компоненты и требуется только одна переменная.

Вычисление электромагнитного поля в дальней зоне с использованием метода FDTD и интеграла Кирхгофа Решение уравнений Максвелла методом FDTD дает возможность вычисления электромагнитного внутри некоторой ограниченной области пространства. Ограничение на размеры области, в которой проводятся вычисления, накладываются ограниченной вычислительной мощность применяемых ЭВМ, а именно тремя факторами: ограниченной скоростью выполнения арифметических операций, ограниченным объемом оперативной памяти и ограниченной скоростью обмена данными между процессором с памятью. Однако существует ряд задач, в которых необходим расчет электромагнитных полей на больших расстояниях от некоторого объекта, излучающего или рассевающего электромагнитное поле. В их числе расчет диаграммы направленности антенн, определение эффективной поверхности отражения радиолокационных целей, решение задач дифракции и. др. В данной статье предложен эффективный способ вычисления полей в дальней зоне с использованием результатов вычислений ближнего поля методом FDTD. Под дальним полем здесь подразумевается поле за пределами вычислительного объема, а ближнее поле здесь означает поле, вычисляемое непосредственно методом FDTD.

Преобразование ближнего поля в дальнее поле с использованием интеграла Кирхгофа Существует несколько методов преобразования ближнего поля в дальнее поле. Все они включают интегрирование по замкнутой поверхности, которая охватывает излучающий или рассеивающий объект. А в остальном имеются значительные отличия. Одни методы интегрируют, используя поверхностные эквивалентные токи (электрические и магнитные), другие непосредственно используют поля Е и Н. Одни методы применяются в частотной области, другие работают во временной области. Например, для вычисления дальнего поля часто применяют интегрирование по элементам, на которые разбивается замкнутая поверхность. Каждый элемент поверхности является элементарным электрическим и магнитным диполем одновременно. Если при этом ближнее поле вычисляется методом FDTD, выполняются следующие операции: - Все компоненты E и H полей на поверхности, по которой происходит интегрирование, сохраняются на всех шагах по времени. Т.е. к концу вычислений известна временная форма поля в каждой точке поверхности. - Поля Е и Н на поверхности преобразуются в эквивалентные электрические и магнитные токи. - Осуществляется преобразование токов на поверхности в частотную область. - Вычисляются составляющие Eφ, Eθ, Еr и соответствующие составляющие поля H для всех частот и от каждого элемента Гюйгенса на поверхности. Элементы Гюйгенса в данном случае - это обычные ячейки FDTD. С каждой ячейкой сопоставляется своя собственная система полярных координат. - Значения Eφ, Eθ, Еr и соответствующих составляющих поля H преобразуются в прямоугольные координаты (Ex, Ey, Еz) и складываются в точке наблюдения. Только на этом шаге происходит собственно интегрирование по поверхности. - Производится обратное преобразование во временную область. Теперь искомое поле найдено.

Существует ряд модификаций приведенного метода расчета дальнего поля, но при применении метода FDTD более логичным выглядит вычисление поля в дальней зоне с использованием только временной области и без преобразования полей в эквивалентные токи. Кроме того, желательно обойтись без преобразования в другие системы координат. Такой метод существует – это вычисления с использованием поверхностного интеграла Кирхгофа. Существует ряд модификаций приведенного метода расчета дальнего поля, но при применении метода FDTD более логичным выглядит вычисление поля в дальней зоне с использованием только временной области и без преобразования полей в эквивалентные токи. Кроме того, желательно обойтись без преобразования в другие системы координат. Такой метод существует – это вычисления с использованием поверхностного интеграла Кирхгофа. Способ применения поверхностного интеграла Кирхгофа совместно с методом FDTD приводится в [3]. Однако в [3] допущены три грубые ошибки в формулах. Кроме того, в порядок самого интегрирования введена путаница. Пришлось заново вывести формулы, воспользовавшись идеей Рамахи (Ramahi, автор [3]). Результат получился хорошим, поэтому, несмотря на досадные ошибки в статье, хочется выразить автору свою благодарность. Интеграл Кирхгофа связывает поле внутри ограниченного объема с полем и его производными на поверхности, ограничивающей объем. Эта формула выведена в середине 19 века немецким физиком Г. Кирхгофом, и во временной области имеет вид [3]: (6) где p, p’ – точка наблюдения и точка на поверхности соответственно; ñ – единичный вектор нормали к поверхности; Ψ – скалярная функция, которая может быть любой из шести компонент поля; R – расстояние от точки наблюдения до точки на поверхности; - вектор p ‑ p’; с – скорость света; A’ – площадь элемента поверхности. В формуле (6) ret означает, что интегрирование осуществляется с учетом запаздывающего времени t’ = t ‑ R/c. Вектор нормали направлен внутрь замкнутого объема 3. O. M. Ramahi, “Near- and far-field calculations in FDTD simulations using. Kirchhoff surface integral representation”, IEEE Transactions on Antennas and Propagation, vol. 45, no 5, may 1997.

Формула (6) выражает принцип Гюйгенса, согласно которому каждая точка на волновом фронте служит фиктивным источником воображаемой сферической волны. Каждый участок поверхности dA’ излучает волну, которая в точку наблюдения p приходит с задержкой R/c. При этом на каждом шаге FDTD по времени на поверхности интегрирования возникает совокупность фиктивных источников, поле от которых придет в точку наблюдения с разным запаздыванием, поскольку расстояние R для всех точек различно. Это означает, что на одном временном шаге FDTD из (6) получаются вклады участков dA’ в разные временные участки выходного сигнала в точке наблюдения. Формула (6) выражает принцип Гюйгенса, согласно которому каждая точка на волновом фронте служит фиктивным источником воображаемой сферической волны. Каждый участок поверхности dA’ излучает волну, которая в точку наблюдения p приходит с задержкой R/c. При этом на каждом шаге FDTD по времени на поверхности интегрирования возникает совокупность фиктивных источников, поле от которых придет в точку наблюдения с разным запаздыванием, поскольку расстояние R для всех точек различно. Это означает, что на одном временном шаге FDTD из (6) получаются вклады участков dA’ в разные временные участки выходного сигнала в точке наблюдения. Шаг по времени при вычислении интеграла (6) тесно связан с шагом по времени FDTD и равен ему. Выходная последовательность в точке наблюдения имеет такой же шаг по времени. Однако задержка R/c может не быть кратной шагу по времени. Поэтому получаемое время задержки округляется до ближайшего времени, кратного шагу FDTD. Возникающая при этом ошибка незначительна, т.к. шаг по времени в классическом FDTD мал по сравнению с периодом колебаний вычисляемого сигнала. Чтобы привести (6) к виду, удобному для применения совместно с алгоритмом FDTD, необходимо выполнить ряд преобразований. При этом учтем, что все участки поверхности, по которой ведется интегрирование, в алгоритме FDTD всегда перпендикулярны одной из осей координат.

Предположим, что поверхность dA’ перпендикулярна оси Z и точка p’ имеет координаты (i, j, ko) (рис. 2). Предположим, что поверхность dA’ перпендикулярна оси Z и точка p’ имеет координаты (i, j, ko) (рис. 2). Применим преобразование подынтегральных слагаемых (6) в конечные разности: (7)

Для упрощения дальнейшей записи введем обозначения: Для упрощения дальнейшей записи введем обозначения: (8) и применим стандартные для FDTD обозначения для дискретных значений времени: Ψ(t’) = Ψ(n+1), Ψ(t’-Δt) = Ψ(n), Ψ(t’+ Δt) = Ψ(n+2). Кроме того, вспомним, что время в точке наблюдения t=t’+R/c. Округляя t до ближайшего целого шага по времени обозначим его как tn*. Выражение (6) с учетом введенных обозначений (7) и (8) запишется для одной площадки dA в виде: (9) где Ai,j = Δx Δy. В (9) нет смысла записывать сумму по всем участкам поверхности, т.к. tn* в общем случае для каждого участка имеет разное значение, поскольку различно расстояние R. Временная последовательность в точке наблюдения p получается суммированием значений Ψ(p, tn*)на каждом шаге по времени по всем участкам поверхности с учетом времени запаздывания.

В (9) значение функции Ψ(p, tn*) вычисляется с использованием значений Ψ на трех шагах по времени. Но это не означает, что необходимо всегда помнить значения двух предыдущих шагов. Сгруппируем слагаемые (9) по временным шагам и запишем в виде: В (9) значение функции Ψ(p, tn*) вычисляется с использованием значений Ψ на трех шагах по времени. Но это не означает, что необходимо всегда помнить значения двух предыдущих шагов. Сгруппируем слагаемые (9) по временным шагам и запишем в виде: (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).

Так что не требуется хранить значения полей на поверхности интегрирования. Но для ускорения вычислений можно заранее вычислить значения расстояний и косинусы углов от всех участков поверхности до точки наблюдения. Так что не требуется хранить значения полей на поверхности интегрирования. Но для ускорения вычислений можно заранее вычислить значения расстояний и косинусы углов от всех участков поверхности до точки наблюдения. Вместо функции Ψ можно свободно подставлять Ex, Ey, Ez, Hx, Hy или Hz. Все компоненты ячейки Yee находятся в разных местах пространства. Но расстояния и косинусы углов вполне можно вычислить общие для всех компонент. Главное, чтобы они принадлежали одной ячейке. Тогда при вычислении дальнего поля рассчитанные компоненты также окажутся сдвинутыми в пространстве, что удобно использовать при тестировании программы, когда поле одновременно, в одной и той же ячейке, вычисляется как непосредственно алгоритмом FDTD, так и решением поверхностного интеграла Кирхгофа. Здесь приведен вывод формулы только для одной поверхности. Для остальных пяти поверхностей вывод формул аналогичен. Интегрирование следует вести строго по замкнутой поверхности. Границы интегрирования можно приблизить вплотную к объекту, заданному в расчетах FDTD. При вычислении полей на довольно больших расстояниях из уравнения (9) можно исключить слагаемое, убывающие как 1/R2, а угол направления на точку наблюдения вычислить один на всю грань поверхности интегрирования. Этот случай соответствует случаю дальней зоны в ее обычном понимании и позволяет ускорить расчеты и уменьшить требуемую для расчетов память.