Расчет переходного сопротивления заземлителей с помощью метода конечных разностей во временной области

Куклин Д.В. Расчет переходного сопротивления заземлителей...

УДК 621.311

Расчет переходного сопротивления заземлителей с помощью метода конечных разностей во временной области

Д.В. Куклин

Центр физико-технических проблем энергетики Севера КНЦРАН, лаборатория высоковольтной электроэнергетики и технологии

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

Abstract. A short review of widespread methods for calculation of transient grounding characteristics and more detail review of the Finite Difference Time Domain method (FDTD) have been given in the paper. Comparison between FDTD and other methods and main features of application FDTD for transient grounding characteristics calculations have been made. An example of transient resistance calculation of metallic rod using FDTD has been represented.

Ключевые слова: метод конечных разностей во временной области, переходное сопротивление, заземление Key words: finite difference time domain method, transient resistance, grounding

1. Введение

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

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

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

2. Методы расчетов

2.1. Обзор методов

Одним из распространенных методов, используемых для расчета переходных характеристик заземлителей сложной конфигурации, является метод создания эквивалентной цепи с сосредоточенными параметрами (Circuit approach) (Liu, 2004).

При использовании метода необходимо проделать следующие этапы:

• разделение заземлителя на множество конечных элементов;

• создание эквивалентной цепи с сосредоточенными параметрами и вычисление их значений (собственной и взаимной индуктивности, емкости, проводимости и внутреннего сопротивления);

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

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

114

Вестник МГТУ, том 16, №1, 2013 г.

стр.114-122

задержку распространения волн (Liu, 2004).

В методе моментов (Method of Moments, MOM) (Liu, 2004) анализируемая структура разделяется на металлические пластины и проводники (в некоторых алгоритмах пластины не применяются; в них сплошные объекты заменяются их каркасными моделями). После того как структура определена, все проводники разделяются на отрезки проводников (они должны быть достаточно короткими по сравнению с длиной волны вследствие допущения, согласно которому ток не меняется при переходе через сегмент), а пластины разделяются на небольшие по сравнению с длиной волны участки. После этого для получившейся структуры создается система линейных уравнений. Через решение этой системы уравнений вычисляются токи для каждого отрезка провода и участка пластины. После нахождения токов можно узнать электрическое поле в любой точке пространства через векторную сумму решений уравнений для каждого отрезка/участка.

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

В методе конечных элементов (Finite Element Method, FEM) (Liu, 2004) вычисления производятся по всему объему. Рассчитываемая область разделяется на небольшие элементы, обычно треугольной или тетраэдрической формы, и представляется как сетка, состоящая из конечных элементов. В этом методе используется допущение, согласно которому поля постоянны внутри каждого элемента сетки. Поле в каждом элементе аппроксимируется полиномом небольшого порядка с неизвестными коэффициентами. Эти аппроксимирующие функции подставляются в вариационное выражение, и получающаяся в результате система уравнений решается для определения неизвестных коэффициентов. После того как коэффициенты найдены, поля внутри каждого элемента известны.

Метод длинной линии (Transmission Line Approach) является первым методом, который был применен для расчета переходных процессов в заземлителях. Однако разработка этого метода не была столь быстрой, как для других методов, и поначалу он мог применяться только для заземлителей простой формы. Он исходит из теории длинных линий и основывается на решении телеграфных уравнений. При современных вычислительных возможностях этот метод получил дальнейшее развитие; совместно с другими методами он может быть применен для заземлителей более сложной формы (Liu, 2004), однако нуждается в экспериментальной проверке.

Метод конечных разностей во временной области (Finite Difference Time Domain, FDTD) (Yee, 1966; Taflove, Hagness, 2000) основан на центрально-разностной дискретизации уравнений Максвелла во времени и пространстве.

2.2. Сравнение метода FDTD с другими методами

Метод FDTD имеет следующие преимущества перед остальными методами:

• простота в понимании;

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

• усложнение формы объектов не ведет к увеличению занимаемой памяти и времени расчета. Объем занимаемой памяти и время расчета могут измениться только при изменении алгоритма расчета (например, при учете частотной зависимости физических свойств, использовании различных сеток);

• отсутствие необходимости задавать граничные условия. Их учет происходит непосредственно при расчете.

Вместе с тем, данный метод имеет свои недостатки:

• кроме самих исследуемых объектов, необходимо моделировать среду, которая их окружает, поэтому объем памяти, занимаемой программой, и время выполнения могут быть слишком большими;

• при использовании ортогональной сетки расчеты, связанные с объектами цилиндрической, сферической или другой неортогональной формы, могут быть неточными. А повышение точности может потребовать существенного усложнения алгоритма расчета;

• сложно моделировать произвольно ориентированные тонкие проводники.

2.3. Метод FDTD с ортогональной сеткой

Уравнения Максвелла для изотропной среды:

115

Куклин Д.В. Расчет переходного сопротивления заземлителей...

Vx H = J +

dD

dt ’

Vx E = -

дБ

dt ’

где

D = ёЕ, B = Щ.

В декартовой системе координат уравнения (1), с учетом материальных соотношений (2), эквивалентны следующей системе уравнений (здесь не учитывается ток переноса):

д ту дУ д j

Обозначим точку пространства как

тогда для функции пространства и времени

F

Запишем уравнение (3а), используя центрально-разностную аппроксимацию во времени и пространстве:

1

dEx 1 "f dHz dHy )-S?x ,

dt e dy dz dt m Уdz dy J

dEy 1 Г

dt e у dz dx dt m У dx dz J

dEz 1 "f dHy dHx' j-sEz ,

dt e У dx dy dt m у ^ dx z

e

ЕХП

At

Ay

Hn-1/2

Az

-s

Вместо значений поля Ex, Ey, Ez на шаге n-1/2 в (Taflove, Hagness, 2000) предлагается использование следующего приближения на примере Ex:

Exn-1/2(i +1/2, j, k) =

Enx (i +1/2, j,k) - Exn-1(i +1/2, j,k) 2

Соответствующим образом записываются уравнения (3б) и (3в). Для (3г) получаем:

Hxn+1/2(i, j +1/2, k +1/2) - Hxn-1/2(i, j +1/2, k +1/2) =_1_______

At _m

Г Eny (i, j +1/2,k +1) - Eny (i, j +1/2,k) Ez”(i, j +1,k +1/2) - Ez”(i, j,k +1/2)^

Az

Ay

Аналогично записываются уравнения (3д) и (3е).

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

En(i +1/2, j,k) = 2 +1/2,.Лk)-s(i +1/2,/,k)'At • Exn-1 (i +1/2, j,k) +

2 e(i +1/2, j,k) + s(i +1/2, j,k)• At

2At

2 e(i +1/2, j, k) + s(i +1/2, j, k )• At

( Hn~1/2(i +1/2, j +1/2, k) - Hzn~1/2(i +1/2, j -1/2, k) ^

Ay

Hn~1/2

Az

116

Вестник МГТУ, том 16, №1, 2013 г.

стр.114-122

е; (/,. ,+1/2, k) =2e((-1+',2-k 1 -s('; j+1/2 k )-Al - еУ//,. /+1/2, k)+

2-e(i, j +1/2, k) + s(i, j +1/2, k)- At

f Hn-1/2/;

2At

2-e(i, j +1/2, k) + s(i, j +1/2, k)- А

(i, j +1/2, k +1/2) - H;-1!2(i, j +1/2, k -1/2) Az

H;-1/2(i +1/2, j +1/2, k) - H;-1/2(i -1/2, j +1/2, k) Ax

E;(i, j,k +1/2) = 2- ^-Ak +1/2) g(i»Ak +1/2)- At - е;-1 (i, j,k +1/2) +

2 - e(i, j,k +1/2) + s(i, j,k +1/2) - A 2At

f H;-1/2(i +1/2, j, k +1/2) - H;-1/2(i -1/2, j, k +1/2) A Ax

2- e(i, j, k +1/2) + s(i, j, k +1/2)- At

H;-1/2

Ay

Hf1/2(i, j +1/2, k +1/2) = H;-1/2(i, j +1/2, k +1/2) +

At

Mi, j +1/2, k +1/2)

re;(i,j +1/2,k +1) -e;(i,j +1/2,k) e;(i,j +1,k+1/2) -e;(i,j,k +1/2)^

Az

H”+1/2(i +1/2, j, k +1/2) = H;-1/2(i +1/2, j, k +1/2) +

Ay

At

Mi +1/2, j, k +1/2)

Г E;(i +1, j,k +1/2) -E;(i, j,k +1/2) E;(i +1/2, j,k +1) -Enx (i +1/2, j,k)Л

Ax

Az

H;+1/2(i +1/2, j +1/2, k) = H;-1/2(i +1/2, j +1/2, k) +

At

^

e;

Ay

Ax

На рис. 1 изображено взаимное расположение векторов электрического и магнитного поля небольшого участка сетки с индексами полей i, j, k и некоторыми соседними точками. Ax, Ay, Az -размеры ячейки FDTD, которые часто принимают равными между собой.

Из рис. 1 видно, что в каждой конкретной точке пространства рассчитывается только одна составляющая электрического или магнитного поля. Следует отметить, что непосредственно при вычислениях индексы i + 1/2, j + 1/2, k + 1/2 невозможны, поэтому они "округляются" до i, j, k, а вместо i - 1/2, j - 1/2, k - 1/2 используются i - 1, j - 1, k - 1.

Время учитывается неявным образом непосредственно при расчете, который удобно совершать в цикле. Одна итерация цикла равна промежутку времени At, и она включает следующую последовательность:

1. Расчет Ex, Ey, Ez по всему объему с использованием формул (9а-9в). Во время расчета предыдущие сохраняемые в массиве значения поля в каждой точке заменяются новыми.

2. Задание источников электрического поля (если имеются) через установку значений электрического поля в необходимых точках.

3. Расчет Hx, Hy, Hz по всему объему с использованием формул (9г-9е). Во время расчета предыдущие значения поля в каждой точке заменяются новыми.

4. Задание источников магнитного поля (если имеются) через установку значения магнитного поля в необходимых точках.

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

2.4. Граничные условия для имитации бесконечной области вычислений

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

117

Куклин Д.В. Расчет переходного сопротивления заземлителей...

Рис. 1. Относительное расположение векторов электрического и магнитного полей

Поскольку в методе FDTD используется преимущественно центрально-разностная аппроксимация, то вычисление поля в конкретной точке требует знания полей в соседних точках. В случае с ортогональной сеткой необходимо знать поля в четырех точках, окружающих данную. Но для расчета полей в точках, находящихся на границах области, информация о значениях полей в некоторых соседних точках отсутствует, т.к. они находятся за границей области вычислений. Отсутствие расчета полей в точках, находящихся на границах области, выражается в виде отражений волн от границ и приводит к большим ошибкам вычислений. Следовательно, необходимо применять специальные меры для того, чтобы избежать отражений от границ. Для этого были созданы так называемые поглощающие граничные условия (Absorbing Boundary Conditions, ABC). Хотя, безусловно, это не граничные условия в обычном смысле этого термина, их задача - не учитывать поведение электромагнитных волн на границе между двумя средами, а препятствовать отражению волн от границ области вычислений.

Существует несколько граничных условий для имитации бесконечной области вычислений, и наиболее часто используемые - это граничные условия Мура, Ляо и PML (Perfectly Matched Layer). Но в последнее время чаще всего используется PML, т.к. это граничное условие обладает на порядки меньшими коэффициентами отражений волн от границ, чем граничные условия Мура и Ляо. Единственное преимущество граничных условий Мура и Ляо в том, что они требуют меньших компьютерных ресурсов. Но PML позволяет располагать границы области гораздо ближе к исследуемым объектам, и за счет этого объем вычисляемой области становится небольшим, поэтому условия Мура и Ляо могут и не давать выигрыша в скорости и используемой памяти. Вместо PML часто используют его разновидность - UPML (Uniaxial Perfectly Matched Layer) (Gedney, 1996). Граничное условие UPML по эффективности идентично PML, но проще в использовании, а также ближе к реальному веществу - его можно сравнить с анизотропным поглощающим электромагнитные волны веществом, окружающим исследуемую область.

2.5. Выбор размера шага сетки и временного интервала вычислений

Параметры Ах, Ду и Д выбираются из соображений точности, а At - стабильности расчета. Величины, определяющие размеры сетки, выбираются таким образом, чтобы они составляли небольшую долю от минимальной необходимой в расчете длины волны (часто берется 1/10 длины волны), или они должны быть достаточными для точного моделирования объектов. Для предотвращения неустойчивости численного решения временной шаг вычисляется по формуле (условие Куранта):

At <

1 1 1 11

c • - + - —- +

1 (Ax)2 (Ay)2 (Az)2

1

где с - скорость света в вакууме.

В случае если Дх = Ду = Az, формула (10) упрощается до

At < Дх / V3c.

118

Вестник МГТУ, том 16, №1, 2013 г.

стр.114-122

2.6. Моделирование тонкого проводника

Для моделирования вещества его свойства (е, а, у) просто назначаются в тех точках пространства, в которых присутствует это вещество. В том случае, если моделируемые объекты имеют "удобную" форму, аппроксимация получается точной. В противном случае расчет будет неточным.

Возможны следующие варианты решения этой проблемы.

1. Выбор удобной сетки, с помощью которой можно точно аппроксимировать форму нужного объекта. Например, это может быть ортогональная сетка с переменным шагом или неструктурированная сетка.

2. Использование методов, позволяющих аппроксимировать "неудобные" объекты с помощью ортогональной сетки. Это достигается за счет того, что поле в отдельных точках рассчитывается с учетом формы объекта и его свойств.

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

Один из наиболее важных объектов - тонкий проводник. Но размеры сетки часто оказываются гораздо больше диаметра проводника, поэтому было создано несколько методов, позволяющих моделировать проводники, диаметр которых меньше шага сетки (Taflove, Hagness, 2000; Railton et al., 2005). Существует два подхода для моделирования тонких проводников, располагающихся вдоль узлов расчетной сетки:

1. Коррекция значений магнитного поля вокруг проводника в близлежащих к проводнику точках

(Taflove, Hagness, 2000).

2. Коррекция значений е, у (и иногда а) для окружающих проводник близлежащих точек, в которых рассчитываются электрическое и магнитное поля (Railton et al., 2005).

Метод, описанный в (Railton et al., 2005), позволяет моделировать тонкий проводник в том числе и в проводящих средах, поэтому он подходит для расчета характеристик заземлителей. Рассмотрим его подробнее.

Рис. 2 иллюстрирует тонкий проводник с радиусом г, ось которого совпадает с компонентой поля Ez(i,j, к). В этом примере радиус проводника меньше Ах и Ау.

Рис. 2. Тонкий проводник, ось которого совпадает с Ez компонентами поля

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

*

ex ~ ex

e = e

°y °y

*

*

fix =Mx

*

my =my

2tan ‘

2 tan-1

ln

2 tan-1(Ax / Ay) Ay ’

ln

2 tan-1(Ay / Ax) Ax ’

119

Куклин Д.В. Расчет переходного сопротивления заземлителей...

где a - радиус проводника; е, ц - первоначальные свойства вещества, окружающего проводник; е , ц скорректированные свойства вещества.

В случае с Ат = Ay формулы (12а-12г) упрощаются:

p/2

*

e = e

1и(Дг /

a)

* 1и

m = m—-------

p/2

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

Ex(i, j, k), Ey(i, j, k), Ex(i, j, k + 1), Ey(i, j, к + 1), Hx(i, j, k), Hy(i, j, k), Hx(i, j - 1, k), Hy(i - 1, j, k), Ex(i - 1, j, k), Ey(i, j - 1, k) (последние два не изображены на рисунке).

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

2.7. Пример расчета переходного сопротивления заземлителя

Среди методов численного решения задач электродинамики метод FDTD оказался довольно эффективным для расчета переходных сопротивлений заземлителей. Точность этого метода для заземлителей различной конфигурации проверена экспериментально и является удовлетворительной (Tanabe, 2001).

Рассчитаем с помощью метода FDTD переходное импульсное сопротивление металлического стержня. Выберем стержень квадратного сечения, т.к. его удобно моделировать в ортогональной сетке. Сечение стержня - 0,5 х 0,5 м2, длина - 2 м. Шаг сетки - 0,25 м, причем Ах, Ау и Az равны между собой. Временной шаг выбираем согласно условию стабильности расчетов по формуле (11). Выбирая максимальное значение At, получаем At = 481,5 пс.

На рис. 3 изображен вид сверху на моделируемую область.

Удаленный токовый электрод

Рис. 3. Вид на расчетную область сверху

Токовый контур образован металлическим стержнем квадратного сечения, источником напряжения, резистором сопротивлением 750 Ом, горизонтальным проводом длиной 20 м и удаленным токовым электродом. Источник напряжения задается через электрическое поле. Форма импульса напряжения определяется следующими выражениями:

nDt < 1,5T

Vs (nDt) = [e~ainM - e

f

-a-nDt 1

“2

An

Vs (nDt) =

)• sin 2

nDt > 1,5T

f

V

Ao

-anDt

-e

-a-nDt

)

t0 =

где aj = 0,69315/Tt, a2 = 2,55843/Tf,

= Цу a2) A = -V0

A0 e

a1 - a2

- e

Юп =

3T,

ж

n

120

Вестник МГТУ, том 16, №1, 2013 г.

стр.114-122

Время фронта Tf равно 0,05 мкс, длина импульса Tt - 500 мкс, амплитуда Vmax - 300 В. На рис. 4 показана форма импульса напряжения.

Потенциальный контур образован проводом длиной 50 м и удаленным потенциальным электродом.

Заземляющий электрод смоделирован как идеальный проводник с нулевой удельной электрической проводимостью. Проводники токового и потенциального контуров, а также удаленные электроды смоделированы как тонкие проводники путем приравнивания к нулю электрического поля вдоль оси проводников. Их эквивалентный радиус примерно равен 0,208Дг = 0,052 м.

Относительная диэлектрическая проницаемость земли принята равной 10. Удельное электрическое сопротивление грунта принято равным 200 Омм.

На рис. 5 изображены yz и xz плоскости рассчитываемой области.

Рис. 5. yz и xz плоскости

Размер области - 80 м х 50 м х 50 м.

Ток, стекающий с заземлителя, вычислялся через интеграл напряженности магнитного поля вдоль замкнутого контура, охватывающего токовый контур в точке, находящейся над вершиной заземляющего электрода:

I(nAt) = | H(nAt)dl.

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

V (nAt) = EZ (nAt) -Az.

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

Метод FDTD требует больших вычислительных ресурсов. Но, используя параллельные вычисления (Yu et al., 2006) в многопроцессорных системах, время расчета может быть существенно снижено, либо может быть увеличена моделируемая область. Хотя при увеличении области вычислений время расчета растет быстрее. Также нужно учитывать, что при большом количестве точек сетки этот

121

Куклин Д.В. Расчет переходного сопротивления заземлителей...

метод требует большого объема оперативной памяти компьютера, но объем занимаемой памяти с увеличением области расчетов растет линейно.

« -

Я

Я ч

о О

и 15 10 5 0

R /

/Л

I

/

\

V

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

Время, мкс

Рис. 6. Результаты вычислений

1

0,9

0,8

0,7

0,6

0,5 < g

0,4 H

0,3

0,2

0,1

0

Большие импульсные токи, создаваемые разрядом молнии, приводят к процессам искрообразования в почве, которые влияют на значение сопротивления заземления. Эти процессы подробно не изучены по причине их сложности. Однако, используя существующие модели процессов искрообразования, их учет возможен, в том числе и в расчетах с помощью метода FDTD (Santos et al., 2009).

3. Заключение

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

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

Учитывая простоту и гибкость метода FDTD, по мнению автора, имеет смысл создание гибридных методов, в которых неоднородности сред учитывались бы методом FDTD, а расчет токов в проводниках - каким-либо другим методом.

Литература

Gedney S.D. An anisotropic perfectly matched layer-absorbing medium for the truncation of FDTD lattices.

IEEE Transactions on Antennas and Propagation, v. 44, p.1630-1639, 1996.

Liu Y. Transient response of grounding systems caused by lightning: Modeling and experiments. Ph.D. Thesis, Uppsala University, Uppsala, p.19-50, 2004.

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, v. 53, p.4129-4136, 2005.

Santos T., Oliveira R., Sobrinho S.S., Almeida J.F. Soil ionization in different types of grounding grids simulated by FDTD method. Microwave and Optoelectronics Conference (IMOC), p.127-132, 2009. Taflove A., Hagness S.C. Computational electrodynamics: The finite-difference time-domain method. Artech House, 852 p., 2000.

Tanabe K. Novel method for analyzing the transient behavior of grounding systems based on the finite-difference time-domain method. Power Engineering Society Winter Meeting, v. 3, p.1128-1132, 2001.

Yee K.S. Numerical solution of initial boundary value problems involving Maxwell’s equation in isotropic media. IEEE Transactions on Antennas and Propagation, v. 14, p.302-307, 1966.

Yu W., Mittra R., Su T., Liu Y., Yang X. Parallel finite-difference time-domain method. Artech House, 262 p., 2006.

122