Оценка погрешностей метода конечных разностей во временной области при моделировании протяженных проводников, разделенных одной ячейкой расчетной сетки Текст научной статьи по специальности «Электротехника, электронная техника, информационные технологии»

УДК 621.396.67

Д. В. Куклин

Оценка погрешностей метода конечных разностей во временной области при моделировании протяженных проводников, разделенных одной ячейкой расчетной сетки

Метод конечных разностей во временной области содержит вспомогательные способы для моделирования проводников, диаметр которых существенно меньше размера ячейки. Среди них существует группа методов, в которых моделирование проводников достигается при помощи корректировки параметров среды, окружающей проводник. Несмотря на то что в таком случае могут моделироваться только проводники, расположенные вдоль узлов расчетной сетки, методы данной группы не приводят к ошибке расчета при расположении проводников в области поглощающих граничных условий, что необходимо для моделирования проводников бесконечной длины. Выполнено сравнение нескольких методов моделирования проводников в случае, когда проводники разделены одной ячейкой расчетной сетки. Проведена оценка влияния корректировки магнитного поля вокруг промежутка между проводниками на результаты расчетов. На примере расчетов с дипольной антенной показано, что корректировка магнитного поля вокруг промежутка между проводниками не приводит к существенному увеличению точности расчетов (по крайней мере, в рассмотренных случаях). В то же время применение данной корректировки устраняет вычислительную неустойчивость. В работе предложено применять корректировку магнитного поля на концах проводника только в тех случаях, когда проводники разделены одной ячейкой, но не у каждого конца проводников, как это предлагалось ранее.

Ключевые слова: метод конечных разностей во временной области (FDTD), вычислительная неустойчивость, дипольная антенна, метод моделирования проводников, расчетные ошибки.

Один из основных способов моделирования протяженных проводников в методе конечных разностей во временной области [1] основан на корректировке параметров среды, окружающей проводник [2; 3]. При моделировании проводников, расположенных на расстоянии, равном одной ячейке расчетной сетки, могут возникать ошибки расчета (вплоть до возникновения вычислительной неустойчивости). Подобное расположение проводников встречается при проведении расчетов с заземлителями [4], антеннами и в прочих случаях [5; 6]. Как правило, между проводниками в таких случаях расположен источник тока или напряжения [4-6], либо проводники разделены промежутком с целью расчета напряжения между ними [5]. Поскольку данное расположение проводников встречается во многих расчетах, необходимо оценить уровень возникающих ошибок и определить возможные способы их устранения.

Материалы и методы

Все расчеты проводились при помощи метода конечных разностей во временной области (FDTD) [1]. Метод FDTD основан на решении двух из четырех уравнений Максвелла в дифференциальной форме:

E - напряженность электрического поля; H - напряженность магнитного поля; D - электрическое смещение; B - магнитная индукция; J - плотность тока проводимости; е - диэлектрическая проницаемость; р. - магнитная проницаемость.

Для устранения отражений волн от границ расчетной области в FDTD применяются т. н. "поглощающие граничные условия", являющиеся вспомогательным численным методом. В качестве поглощающих граничных условий в данной статье применяется convolutional perfectly matched layer (CPML) [1]. Толщина поглощающих граничных условий - 10 ячеек.

Условие Куранта (необходимое условие устойчивости численного решения) для метода FDTD с расчетной сеткой из кубических ячеек [1]

Введение

(1)

где

D = е E, B = ,

(2)

где dt - временной шаг; А - размер ячейки; с - скорость света в вакууме. Запишем формулу расчета временного интервала в виде

где £ - безразмерный коэффициент, который не должен превышать единицу для получения устойчивого решения.

Для моделирования проводников, диаметр которых существенно меньше шага расчетной сетки, необходимо применение дополнительных методов [1]. В одном из наиболее распространенных подходов проводники располагаются вдоль узлов, в которых рассчитывается электрическое поле. Однако после задания в модели удельной проводимости протяженных проводников вдоль узлов расчетной сетки необходимо точно задать диаметр проводника. Для точного моделирования диаметра проводников корректируются значения электрического и магнитного полей в узлах расчетной сетки, находящихся в непосредственной близости от проводника. Причем корректировка производится таким образом, чтобы характер электрического и магнитного полей соответствовал полям вокруг реальных проводников [1]. В методе FDTD корректировку электрического и магнитного полей удобно производить путем изменения параметров среды. В работах [2; 3] рассматривались методы моделирования проводников. В методе [2] для моделирования проводника радиусом г корректируются значения диэлектрической и магнитной проницаемости в соответствии с формулами

Значения параметров среды корректируются в конечно-разностных уравнениях метода FDTD, при помощи которых рассчитываются значения электрического и магнитного полей. Диэлектрическая проницаемость корректируется в узлах, радиальных к проводнику, магнитная проницаемость - в тангенциальных к проводнику узлах. В работе [3] предлагается модификация предыдущего метода, в соответствии с которой корректируются также компоненты магнитного поля, продольные по отношению к проводнику в тех случаях, когда радиус проводника меньше значения 0,208-А. Когда радиус проводника больше 0,208-А, корректируются продольные компоненты электрического поля. Это необходимо для устранения вычислительной неустойчивости при значениях коэффициента Б, близких к единице. Также в работе [3] предложено корректировать магнитное поле на расстоянии А/2 от концов проводника, если радиус проводника меньше значения 0,208-А.

Для проведения расчетов использовалась компьютерная программа, разработанная автором, обладающая основными возможностями метода FDTD и использующая дополнительные методы для моделирования протяженных проводников, источников, сред с частотной зависимостью диэлектрической проницаемости и др.

Результаты и обсуждение

Ошибки расчетов

Для демонстрации ошибок расчета проведены расчеты с дипольной антенной длиной 0,21 м (рис. 1). Дипольная антенна состоит из двух проводников, между которыми расположен источник тока. В расчетах использовался источник тока с внутренним сопротивлением 50 Ом. Диаметр проводников - 0,5 мм. Значение коэффициента Б, используемое в расчетах, - 0,588. Расчеты проводились до момента времени 0,5 мкс.

(4)

г' = г- т, ц ' = ц / т,

при расчетной сетке из кубических ячеек значение т находится как

п

т =-.

(5)

(6)

источник токя

проводники

у'

W

X

поглощающие граничные условия

Рис. 1. Расчетная модель с дипольной антенной Fig. 1. The calculation model with a dipole antenna

Рассчитывалось значение возвратных потерь антенны (return loss, RLin) для пятидесятиомного кабеля

V (©) - 50/(©)

RLm =-201og10 |Sn| = -20 log.

V (©) + 50 / (©)

(7)

где У(ю) и 1(ю) - входные напряжение и ток. Значения У(ю) и 1(ю) получены после расчетов У(() и /(/) и применения к ним преобразования Фурье.

Форма тока источника задавалась функцией следующего вида [7]:

at)

n 141

T

7 exP

+ 0.5т,

(8)

где п = 10; T = 10-8; т = 3 10-9; п = 1; ¡тах = 1 А. Ток задается после его выражения через плотность тока проводимости в уравнениях Максвелла.

Корректируемые компоненты электрического и магнитного полей представлены на рис. 2 на примере центральной части дипольной антенны, для которой проводились расчеты. Компоненты И2, Ну, Ez, Ey корректируются согласно методу [2], за исключением компонентов, которые находятся вокруг источника тока. Компоненты Ey на рисунке не показаны, все компоненты корректируются в близлежащих к проводнику ячейках исходя из симметрии. Нх - компоненты магнитного поля - предложено корректировать в работе [3]. Также в этой работе предложено корректировать компоненты магнитного поля вокруг оси проводника, отстоящие на А/2 от концов проводника; в данном случае они включают в себя компоненты вокруг источника тока (дополнительно корректируемые компоненты на рис. 2).

Рис. 2. Корректируемые компоненты магнитного поля Fig. 2. Correction of the magnetic field

n

t

т

Напряжение Щ) рассчитывается как интеграл напряженности электрического поля вдоль определенного

пути:

V (Т) = | Е (Т )Ш. (9)

В данном случае оно рассчитывается как произведение напряженности электрического поля между проводниками на расстояние между проводниками: E•А. Ток /(/) рассчитывается как интеграл напряженности магнитного поля вдоль контура вокруг источника тока:

I (Т) = | Н (Т )Ш. (10)

Поскольку в связи с особенностями метода расчеты выполнены во временной области, а параметр ЯЬЫ следует вычислять в частотной области, то к найденному току и напряжению, как упоминалось выше, необходимо применить преобразование Фурье и далее для каждой частоты рассматриваемого диапазона частот рассчитать значение по формуле (7).

В расчетах использованы сетки с тремя разными размерами ячеек: 10 мм; 5,12 мм; 2,593 мм. Для ячейки 10 мм проведены расчеты с корректировкой магнитного поля только у тех концов проводников, где расположен источник тока (рис. 2). Результаты расчетов представлены на рис. 3-5.

Рис. 3. Результаты расчета RLin для ячейки размером 10 мм: 1 - метод, предложенный Railton и др.; 2 - метод, предложенный Taniguchi и др.; 3 - метод, предложенный Taniguchi и др. с корректировкой

магнитного поля только у источника тока Fig. 3. The calculation results of the RLin for the cell with the size equal to 10 mm: 1 - the method proposed by Railton et al.; 2 - the method proposed by Taniguchi et al.; 3 - the method proposed by Taniguchi et al. with the correction of the magnetic field only around the current source

w

4

J (X

0 -2 -4

-6 - 8

- 10 - 12

- 14

- 16

ill;

\ : ;

\ ! f !

; ч // ........: А...... .....// ....................

X ! \\ i

\ i

i X V 1 "1 V2-------- 1 1 1 1

0,5 0,55

0,6 0,65 0,7 Частота (ГГц)

0,75 0,8

Рис. 4. Результаты расчета RLin для ячейки размером 5,12 мм: 1 - метод, предложенный Railton и др.; 2 - метод, предложенный Taniguchi и др. Fig. 4. The calculation results of the RLin for the cell with the size equal to 5.12 mm: 1 - the method proposed by Railton et al.; 2 - the method proposed by Taniguchi et al.

0

-2

-4

W п -6

3 - 8

л

- 10

- 12

- 14

- 16

---J 1111

\ / 1 --

1 V 2 ------- .....

0,5 0,55

0,6 0,65 0,7 Частота (ГГц)

0,75 0,8

Рис. 5. Результаты расчета RLin для ячейки размером 2,593 мм: 1 - метод, предложенный Railton и др.; 2 - метод, предложенный Taniguchi и др. Fig. 5. The calculation results of the RLin for the cell with the size equal to 2.593 mm: 1 - the method proposed by Railton et al.; 2 - the method proposed by Taniguchi et al.

Из рисунков видно, что корректировка магнитного поля вокруг источника заметно повлияла на результаты расчетов при крупном размере ячейки. При уменьшении размера ячейки данное влияние существенно уменьшается. При сравнении результатов на рис. 3 с более точными результатами на рис. 5

(за счет меньшего размера ячейки сетки) можно отметить, что в отличие от [8] корректировка магнитного поля вокруг источника существенно не увеличила точность расчетов (нельзя исключить однако, что этот вывод может зависеть от рассчитываемого параметра антенны). Наличие корректировки магнитного поля на концах проводников, не присоединенных к источнику тока, не повлияло на результаты расчетов.

Вычислительная неустойчивость

Следует отметить, что отсутствие корректировки магнитного поля вокруг источника может приводить к вычислительной неустойчивости. Ниже представлены расчеты интеграла напряженности магнитного поля вдоль контура вокруг источника (что равно полному току через поверхность, ограниченную данным контуром).

В первом примере проведены расчеты с небольшим значением коэффициента S. Использована прежняя модель, но с диаметром проводников 0,5 10-3 мм. Проводники моделировались при помощи метода [2]. Параметры функции (8): n = 10, T = 10-6, т = 3• 10-7, п = 1. В качестве источника использовался идеальный источник тока. Размер ячейки - 10 мм, коэффициент S - 0,294. Конечно, данные параметры импульса и диаметр нереалистичны для антенны длиной 0,21 м, но они применяются лишь для демонстрации вычислительной неустойчивости при использовании модели, близкой к той, что была использована в предыдущих расчетах. Результат расчета без корректировки магнитного поля представлен на рис. 6. Из рисунка видно, что результаты расчета ошибочны.

°'8 -1-!-!-!-1-

0,6

<

м" 0,4 о Н

0,2 0

0 0,2 0,4 0,6 0,8 1 1,2 Время (мкс)

Рис. 6. Результат расчета тока без корректировки магнитного поля вокруг источника тока Fig. 6. Result of the current calculation. The magnetic field around the source is not corrected

Далее проведен расчет с применением корректировки компонентов магнитного поля. Корректируемые компоненты магнитного поля показаны на рис. 2. Результаты представлены на рис. 7. Можно видеть, что такой большой ошибки, как в предыдущем случае, нет. Причем поскольку проводники в данных двух случаях моделировались при помощи метода [2], продольные по отношению к проводнику компоненты магнитного поля не корректировались. Так как явная вычислительная неустойчивость устранена (корректировкой компонентов магнитного поля вокруг промежутка между проводниками), то можно предположить, что причина неустойчивости на рис. 3 отличается от причины неустойчивости, описанной в [3].

0,8 -|-|-|-|-1-

<

н

0 Z_i_i_i_i_i_

0 0,2 0,4 0,6 0,8 1 1,2

Время (мкс)

Рис. 7. Результат расчета тока с корректировкой магнитного поля вокруг источника тока Fig. 7. Result of the current calculation. The magnetic field around the source is corrected

В следующем примере неустойчивость возникает в самом начале расчетов. В данном случае проводники моделировались при помощи метода [3]. Параметры функции тока те же, что и в предыдущем случае: диаметр проводников - 0,5 мм; размер ячейки - 10 мм; коэффициент S - 0,98. В результате на рис. 8

видна явная ошибка. Конечно, устранения этой ошибки можно было бы добиться путем уменьшения значения коэффициента но это привело бы к увеличению времени расчета (пропорционально тому, во сколько раз уменьшено значение 5). Из результата расчета с применением корректировки магнитного поля на рис. 9 можно видеть, что вычислительная неустойчивость отсутствует.

Рис. 8. Результат расчета тока без корректировки магнитного поля вокруг источника тока Fig. 8. Result of the current calculation. The magnetic field around the source is not corrected

Рис. 9. Результат расчета тока с корректировкой магнитного поля вокруг источника тока Fig. 9. Result of the current calculation. The magnetic field around the source is corrected

Также проведены расчеты для случая, когда диаметры двух проводников, разделенных промежутком, различны. Расчеты выполнены на прежней модели, но с разными диаметрами двух проводников - 2 мм и 0,008 мм. Проведено два расчета, в которых магнитное поле корректировалось в соответствии с одним из проводников. Результаты (рис. 10) демонстрируют отсутствие вычислительной неустойчивости в обоих случаях: если магнитное поле вокруг источника корректируется в соответствии 1) с проводником большего диаметра и 2) с проводником меньшего диаметра.

Рис. 10. Результат расчета тока с корректировкой магнитного поля вокруг источника тока (и проводниками разного диаметра): 1 - магнитное поле корректируется в соответствии с проводником диаметром 2 мм; 2 - в соответствии с проводником диаметром 0,008 мм Fig. 10. Result of the current calculation for the wires with a different diameter: 1 - the magnetic field around the source is corrected according to the wire with the 2 mm diameter; 2 - according to the wire with the 0.008 mm diameter

Корректировать магнитное поле следует и в тех случаях, когда в промежутке между проводниками нет источника, например, когда рассчитывается разность потенциалов между проводниками.

Стоит отметить, что после многочисленных расчетов с проводниками рассмотренная вычислительная неустойчивость была выявлена только в случаях, когда проводники расположены на расстоянии, равном одной ячейке. Возможно, что в корректировке магнитного поля на концах проводников в остальных случаях нет необходимости, либо такие случаи редки.

Заключение

Для устранения вычислительной неустойчивости при расположении проводников на расстоянии, равном одной ячейке расчетной сетки, вокруг разделяющего проводники промежутка необходимо корректировать магнитное поле точно так же, как и вокруг проводников. При этом нет необходимости корректировать магнитное поле на расстоянии А/2 от концов проводника во всех случаях, когда радиус проводника меньше значения 0,208-А, как это предложено в [3].

Данную корректировку следует применять также и при использовании метода [2], т. е. если продольные по отношению к проводнику компоненты магнитного поля не скорректированы.

Применение этой корректировки не только устраняет вычислительную неустойчивость, но и влияет на результаты расчетов в целом (что не приводит к существенному повышению точности в случае расчета возвратных потерь). Поэтому корректировку предлагается использовать преимущественно для устранения неустойчивости.

Библиографический список

1. Taflove A., Hagness S. C. Computational electrodynamics: the finite-difference time-domain method. 3rd ed. Boston ; London : Artech House, 2005. 1006 p.

2. Railton C. J., Paul D. L., Craddock I. J., Hilton G. S. The treatment of geometrically small structures in FDTD by the modification of assigned material parameters // IEEE Transactions on Antennas and Propagation. 2005. V. 53, Iss. 12. P. 4129-4136. DOI: https://doi.org/10.1109/tap.2005.860008.

3. Taniguchi Y., Baba Y., Nagaoka N., Ametani A. An improved arbitrary-radius-wire representation for FDTD electromagnetic and surge calculations // International Conference on Power Systems Transients (IPST2009). 3-6 June 2009. Kyoto, Japan, 2009. URL: http://www.ipstconf.org/papers/ Proc_IPST2009/09IPST003.pdf.

4. Kuklin D. Choosing configurations of transmission line tower grounding by back flashover probability value // Frontiers in Energy. 2016. V. 10, Iss. 2. P. 213-226. DOI: https://doi.org/10.1007/s11708-016-0398-6.

5. Куклин Д. В. Оценка параметров измерительного устройства электрических характеристик грунта // Труды Кольского научного центра РАН. Энергетика. 2018. № 8 (9), вып. 17. С. 60-67.

6. Куклин Д. В. Коррекция ошибок, вызванных применением ступенчатой аппроксимации проводника при использовании метода конечных разностей во временной области // Вестник МГТУ. 2013. Т. 16, № 4. С. 728-733.

7. Heidler F., Cvetic J. A class of analytical functions to study the lightning effects associated with the current front // International Transactions on Electrical Energy System. 2002. V. 12, Iss. 2. P. 141-150. DOI: https://doi.org/10.1002/etep.4450120209.

8. Watanabe S., Taki M. An improved FD TD model for the feeding gap of a thin-wire antenna // IEEE Microwave and Guided Wave Letters. 1998. V. 8, Iss. 4. P. 152-154.

References

1. Taflove A., Hagness S. C. Computational electrodynamics: the finite-difference time-domain method. 3rd ed. Boston ; London : Artech House, 2005. 1006 p.

2. Railton C. J., Paul D. L., Craddock I. J., Hilton G. S. The treatment of geometrically small structures in FDTD by the modification of assigned material parameters // IEEE Transactions on Antennas and Propagation. 2005. V. 53, Iss. 12. P. 4129-4136. DOI: https://doi.org/10.1109/tap.2005.860008.

3. Taniguchi Y., Baba Y., Nagaoka N., Ametani A. An improved arbitrary-radius-wire representation for FDTD electromagnetic and surge calculations // International Conference on Power Systems Transients (IPST2009). 3-6 June 2009. Kyoto, Japan, 2009. URL: http://www.ipstconf.org/papers/ Proc_IPST2009/09IPST003.pdf.

4. Kuklin D. Choosing configurations of transmission line tower grounding by back flashover probability value // Frontiers in Energy. 2016. V. 10, Iss. 2. P. 213-226. DOI: https://doi.org/10.1007/s11708-016-0398-6.

5. Kuklin D. V. Otsenka parametrov izmeritel'nogo ustroystva elektricheskikh kharakteristik grunta [Estimation of parameters of measurement device for electrical soil] // Trudy Kol'skogo nauchnogo tsentra RAN. Energetika. 2018. N 8 (9), vyp. 17. P. 60-67.

6. Kuklin D. V. Korrektsiya oshibok, vyzvannykh primeneniyem stupenchatoy approksimatsii provodnika pri ispol'zovanii metoda konechnykh raznostey vo vremennoy oblasti [Correction of errors caused by the staircase approximation of wires in the FDTD method] // Vestnik MGTU. 2013. V. 16, N 4. P. 728-733.

7. Heidler F., Cvetic J. A class of analytical functions to study the lightning effects associated with the current front // International Transactions on Electrical Energy System. 2002. V. 12, Iss. 2. P. 141-150. DOI: https://doi.org/10.1002/etep.4450120209.

8. Watanabe S., Taki M. An improved FD TD model for the feeding gap of a thin-wire antenna // IEEE Microwave and Guided Wave Letters. 1998. V. 8, Iss. 4. P. 152-154.

Сведения об авторе

Куклин Дмитрий Владимирович - ул. Ферсмана, 14, г. Апатиты, Мурманская обл., Россия, 184209; Центр физико-технических проблем энергетики Севера - филиал Федерального исследовательского центра КНЦ РАН, науч. сотрудник; e-mail: kuklindima@gmail.com

Kuklin D. V. - 14, Fersmana Str., Apatity, Murmansk region, Russia, 184209; Northern Energetics Research Centre - Branch of the Federal Research Centre KSC RAS, Researcher; e-mail: kuklindima@gmail.com

D. V. Kuklin

Calculation of errors in the finite difference time domain method when modeling thin wires separated by a single cell of the computational grid

The finite difference time domain method contains additional techniques for modeling thin wires whose diameter is substantially smaller than the cell size. Among them, there is a group of techniques where thin wires are modeled by correcting the medium surrounding the wire. Although in this case only the wires located along the electric field nodes can be modeled, techniques of this group do not cause calculation errors if the wires are located in the absorbing boundary conditions necessary to model thin wires with infinite length. Several techniques have been compared for the case when two wires are separated by one-cell gap. Influence of magnetic field modification around the gap between wires on calculation results has been examined. By calculations with the dipole antenna, it has been shown that correction of the magnetic field around the gap does not increase calculation accuracy (at least for the considered cases). However, this correction eliminates numerical instability. In this paper, it has been proposed to apply a magnetic field correction at the ends of the conductor only in cases where the conductors are separated by one cell, but not at each end of the conductors, as previously suggested.

Key words: finite difference time domain (FDTD) method, numerical instability, dipole antenna, wire modeling technique, calculation errors.