1 Московский энергетический институт
 Скачать 3.23 Mb.
|
|
Глава 2. ЧИСЛЕННЫЕ МЕТОДЫ РАСЧЕТА Все макроскопические электромагнитные явления, наблюдаемые на практике, могут быть математически описаны полным набором уравнений Максвелла. Метод конечного интегрирования ( Finite Integration Technique, FIT ) представляет собой последовательную схему дискретизации уравнений Максвелла в интегральной форме [1,2]. Получаемые матричные уравнения дискретизированных полей могут быть использованы для численного моделирования на современных компьютерах. Кроме того, алгебраические свойства этой теории дискретного электромагнитного поля позволяют аналитически и алгебраически использовать законы сохранения энергии и заряда для дискретной формулировки и дают стабильный алгоритм расчета численным методом во временной области. Формулировка метода конечного интегрирования - это общий подход, и поэтому может быть применен во всем частотном диапазоне, от постоянного тока до высоких частот. Исторически, этот универсальный метод уже был в пакете программ MAFIA , разработка которого было начата компанией CST более 20 лет назад. На основании этого опыта, в 1997 году была начата разработка программ семейства STUDIO . Была выполнена модернизация интерфейса и трехмерной визуализации и модифицировано ядро расчета. Однако основные улучшения коснулись стратегии разбиения на ячейки, реализации техники идеальной аппроксимации границы ( Perfect Boundary Approximation, PBA ), расширенной техникой тонкого листа ( Thin Sheet Technique, ТSТ ) и иерархической схемой нанесения подсетки ( Multilevel Subgridding Scheme, MSS ). Метод конечного интегрирования доказал свою эффективность в моделировании электромагнитных явлений. Его можно рассматривать как обобщение и дальнейшее развитие предложенного в 1966 году Кейном Йи ( Kane Yee ) метода конечных разностей во временной области ( Finite Difference Time Domain, FDTD ) [3], описываемого в п. 2.1.1. 2.1. Метод конечных разностей во временной области Метод конечных разностей во временной области ( FDTD ) используется для решения системы уравнений Максвелла во временной области. Несмотря на сложность этих уравнений и наличие в них частных 48 производных, их можно преобразовать в форму, удобную для численных расчетов. Для изотропной среды в отсутствие магнитных токов основные уравнения Максвелла могут быть записаны в следующей форме: 0 = × ∇ + ∂ ∂ E t B r r , (2.1) J H t D r r r = × ∇ − ∂ ∂ , (2.2) H B r r μ = , (2.3) E D r r ε = . (2.4) где вектор электрического тока J r , диэлектрическая ε и магнитная μ проницаемости среды считаются известными функциями пространства и времени, а ∇ - векторный дифференциальный оператор. В декартовой системе координат имеем ⎟⎟ ⎠ ⎞ ⎜⎜ ⎝ ⎛ ∂ ∂ ∂ ∂ ∂ ∂ = ∇ z y x , , Отметим, что для произвольного вектора ( ) z y x F F F F , , = r справедливо ⎟⎟ ⎠ ⎞ ⎜⎜ ⎝ ⎛ ∂ ∂ − ∂ ∂ + ⎟ ⎠ ⎞ ⎜ ⎝ ⎛ ∂ ∂ + ∂ ∂ − + ⎟⎟ ⎠ ⎞ ⎜⎜ ⎝ ⎛ ∂ ∂ − ∂ ∂ = ∂ ∂ ∂ ∂ ∂ ∂ = ≡ × ∇ y F x F z z F x F y z F y F x F F F z y x z y x F F x y x z y z z y x 0 0 0 0 0 0 rot r r r r r r r r , где 0 xr , 0 yr , 0 zr - единичные орты вдоль координатных осей. Тогда из (2.1) и (2.2) для декартовой системы координат получаем z E y E t B y z x ∂ ∂ − ∂ ∂ = ∂ ∂ − , (2.5) x E z E t B z x y ∂ ∂ − ∂ ∂ = ∂ ∂ − , (2.6) x E y E t B y x z ∂ ∂ − ∂ ∂ = ∂ ∂ , (2.7) x y z x J z H y H t D − ∂ ∂ − ∂ ∂ = ∂ ∂ , (2.8) 49 y z x y J x H z H t D − ∂ ∂ − ∂ ∂ = ∂ ∂ , (2.9) z x y z J y H x H t D − ∂ ∂ − ∂ ∂ = ∂ ∂ . (2.10) Дискретизируем пространство задачи и время. Обозначим точку пространства как ( ) ( ) z k y j x i k j i Δ ⋅ Δ ⋅ Δ ⋅ = , , , , , а для любой функции пространства и времени примем такую запись: ( ) ( ) k j i F t n z k y j x i F n , , ; , , = Δ Δ Δ Δ Применив конечно-разностную аппроксимацию к уравнениям (2.5) и (2.8), можем получить ( ) ( ) ( ) ( ) ( ) ( ) y k j i E k j i E z k j i E k j i E t k j i B k j i B n z n z n y n y n x n x Δ + − + + − Δ + − + + = = Δ + + − + + − + 2 1 2 1 2 1 2 1 2 1 2 1 2 1 2 1 , , , 1 , , , 1 , , , , , , 2 1 2 1 (2.11) и ( ) ( ) ( ) ( ) ( ) ( ) ( ) , , , , , , , , , , , , , , 2 1 2 1 2 1 2 1 2 1 2 1 2 1 2 1 2 1 2 1 1 2 1 2 1 2 1 2 1 2 1 2 1 k j i J z k j i H k j i H z k j i H k j i H t k j i D k j i D n x n z n y n z n z n x n x + − Δ − + − + + − − Δ − + − + + = = Δ + − + − − − − − − (2.12) Выражения, соответствующие уравнениям (2.6), (2.9) и (2.7),(2.10) получаются аналогично. При этом значения компонент поля на каждом шаге находятся по значениям на предыдущих шагах. Именно такой итерационный процесс положен в основу расчета методом FDTD . На рис. 2.1 показан кубик пространства (или т.н. элементарная ячейка Йи [3]) с компонентами E и H на гранях, иллюстрирующая связь компонент поля в конечно-разностных уравнениях Максвелла. Отметим, что последние содержат обычные операции сложения/вычитания и умножения/деления, элементарно реализуемые на компьютере. Для проведения расчетов требуется разбить пространство задачи на элементарные ячейки и установить начальные значения для всех компонент полей, которые определяются условиями возбуждения. 50 Рис. 2.1. Представление компонент поля в кубической ячейке Йи Размер пространственной сетки разбиения должен быть таким, чтобы вблизи одного элемента разбиения электромагнитное поле не претерпевало значительных изменений. Это означает, что для получения адекватных результатов линейные размеры сетки должны составлять доли длины волны. Для получения устойчивости результатов вычислений необходимо удовлетворять соотношению между пространственным приращением и временным приращением ∆ t . В случае пространства задачи с непостоянными ε и μ , строгий критерий стабильности получить непросто. Для постоянных же ε и μ стабильные вычисления требуют ( ) ( ) ( ) εμ t t c z y x Δ = Δ > Δ + Δ + Δ 2 2 2 , (2.13) где с – это скорость света в среде. Если с max является максимальной скоростью распространения света в рассматриваемой области задачи, необходимо выбирать ( ) ( ) ( ) t c z y x Δ > Δ + Δ + Δ max 2 2 2 . (2.14) Это требование вводит ограничение на ∆ t для выбранных ∆ x , ∆ y и ∆ z Положив для примера с м 10 3 8 0 max ⋅ = = c c и м 10 мм 10 2 − = = Δ = Δ = Δ z y x , из (2.14) можем получить следующее ограничение на шаг ∆ t : 51 ( ) ( ) ( ) пс 60 с 10 6 с м 10 3 м 10 3 11 8 2 max 2 2 2 = ⋅ ≈ ⋅ ⋅ = Δ + Δ + Δ < Δ − − c z y x t Это означает, что временной процесс, который длится 1 мсек, требует 17 тыс. шагов. 2.2. МЕТОД КОНЕЧНОГО ИНТЕГРИРОВАНИЯ Метод конечного интегрирования (МКИ), предложенный Томасом Вейландом в 1977 году [1], представляет собой дискретную формулировку уравнений Максвелла в интегральной форме, удобную для реализации на компьютерах и позволяющую моделировать реальные электромагнитные задачи со сложной геометрией. Первый шаг формализации МКИ состоит в ограничении электромагнитной задачи, которая обычно представляет собой задачу с открытыми границами, ограниченной областью 3 R ∈ Ω , содержащей пространственную область задачи. Следующий шаг заключается в разбиении расчетной области Ω на конечное число ячеек V i , таких как тетраэдральные (четырехгранные) или гексагональные (шестигранные) при условии, что все ячейки точно прилегают друг к другу, то есть пересечение двух различных ячеек либо отсутствует, либо должно быть двухмерным многоугольником, общей одномерной гранью обеих ячеек или точкой. Это разбиение дает конечную группу ячеек G , играющую роль расчетной сетки. Начиная с этого подхода, основанного на ячейках, подхода к пространственной дискретизации, можно увидеть, что методика конечного интегрирования не ограничивается только трехмерными декартовыми сетками. Она позволяет рассматривать все типы координатных сеток - ортогональные (прямоугольные) и неортогональные. Отметим, что каждая грань ячеек имеет исходную ориентацию, т.е. направление, так что объединение всех этих ячеек может быть описано как направленный граф. Аналогично многоугольные грани ячеек будут связаны с направлением. Для простоты при последующем описании МКИ примем, что Ω имеет форму куба и разбиение на сетку вводится для декартовой системы координат так, что мы получаем набор ячеек [ ] [ ] [ ] }, 1 , 1 , 1 , 1 , 1 , 1 , , , , | { 1 1 1 , , 3 , , − = − = − = × × = ∈ = + + + K k J j I i z z y y x x V R V G i i i i i i k j i k j i (2.15) где узловые точки ( x i , y j , z k ) пронумерованы в соответствии с координатами i,j и k вдоль осей X, Y и Z. Это приводит к общему количеству точек Np = I·J·K для ( I− 1) · ( J− 1) · ( K− 1) ячеек сетки. 52 После определения группы ячеек G сетки, использование теории конечного интегрирования требует рассмотрения области одной ячейки V n Формулировка закона Фарадея в интегральной форме ( ) ( ) , , , 3 R A A d t r B t s d t r E A A ∈ ∀ ⋅ ∂ ∂ − = ⋅ ∫∫ ∫ ∂ r r r r r r (2.16) может быть переписана для грани A z (i,j,k) ячейки V n как обыкновенное дифференциальное уравнение ( ) ( ) ( ) ( ) ( ) k j i b dt d k j i e k j i e k j i e k j i e z y x y x , , , , , 1 , , , 1 , , )) ) ) ) ) − = − + − + + , (2.17) как показано на рис. 2.2, где скалярная величина ( ) ( ) ( ) ∫ + ⋅ = k j i k j i z y x z y x x s d E k j i e , , , , 1 , , r r ) (2.18) является электрическим напряжением вдоль одного ребра поверхности A z (i,j,k) , представляющим точное значение интеграла от электрического поля вдоль этой грани. Скалярная величина ( ) ( ) ∫ ⋅ = k j i A z z A d B k j i b , , , , r r )) (2.19) является магнитным потоком, т.е. величиной плотности магнитного потока через грань ячейки A z (i,j,k). Заметим, что направления граней ячейки влияют на знаки в (2.17). Подчеркнем, что уравнение (2.17) представляет собой точное представление выражения (2.16) для рассматриваемой поверхности ячейки. Рис. 2.2. Ячейка 1 , , − k j i V группы ячеек G с указанными электрическими напряжениями e) на ребрах A и магнитными потоками ребер b )) через эту поверхность x z y 53 Интегральная формулировка закона Фарадея (2.16) справедлива для каждой отдельной грани A(i,j,k) из G , а дискретное приближение в (2.17), естественно, распространяется на б о льшие участки граней ( ) k j i A A , , ∪ = вследствие соотношения ( ) ∑∫ ∫ = A A k j i , , . Аналогичный результат справедлив для поверхностных интегралов. Это и позволяет выбрать такой подход пространственной дискретизации конечным набором ячеек для использования в МКИ. Задавшись порядком электрических напряжений ( ) k j i e , , ) и магнитных потоков через ребра ( ) k j i b , , )) всего набора ячеек G и их представлением в виде векторов-столбцов таким образом, что сформируем сначала составляющие вдоль оси X, затем вдоль Y- и Z-направлений, получаем два вектора-столбца ( ) , R e | e | e e p p N 3 T N ,..., 1 n n , z n , y n , x ∈ = = ) ) ) ) (2.20) R b | b | b b p p N 3 T N ,..., 1 n n , z n , y n , x ∈ ⎟ ⎠ ⎞ ⎜ ⎝ ⎛ = = )) )) )) )) (2.21) Уравнения (2.17) для всех поверхностей ячеек сетки из набора G могут быть записаны в матричной форме { b b dt d e e e e e n n n n n )) )) 3 2 1) ) ) ) ) 4 4 4 4 4 3 4 4 4 4 4 2 1 ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎠ ⎞ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎝ ⎛ − = ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎠ ⎞ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎝ ⎛ ⎟ ⎟ ⎟ ⎠ ⎞ ⎜ ⎜ ⎜ ⎝ ⎛ − − C 1 1 1 1 4 3 2 1 . (2.22) Матрица С содержит только топологическую информацию об отдельных связях ребер ячейки в G и об их ориентации и состоит из коэффициентов { } 1 , 0 , 1 , − ∈ j i C . Она представляет собой дискретный оператор ротора на сетке G В качестве второго дискретного оператора рассмотрим оператор дивергенции. Его ввод следует из уравнения Максвелла, описывающего отсутствие магнитных зарядов: ( ) , 0 , 3 R V A d t r B V ∈ ∀ = ⋅ ∫ ∫ ∂ r r r (2.23) 54 которое рассматривается для ячейки k j i V , , , как показано на рис. 2.3. Рис. 2.3. Распределение шести магнитных потоков через грани, которые рассматриваются при оценке интегрирования по замкнутой поверхности при отсутствии магнитных зарядов внутри ячейки Оценка для интеграла по поверхности (2.20) для изображенной кубической ячейки дает выражение ( ) ( ) ( ) ( ) ( ) ( ) 0 1 k ,j ,i b k ,j ,i b k , 1 j ,i b k ,j ,i b k ,j , 1 i b k ,j ,i b z z y y x x = + + − + + − + + − )) )) )) )) )) )) , (2.24) которое является строгим для рассматриваемого объема. Это соотношение для одной ячейки тоже может быть расширено для всего набора ячеек G , что позволяет получить матрицу дискретной дивергенции 0 S 1 1 1 1 1 1 6 5 4 3 2 1 = ⎟⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎠ ⎞ ⎜⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎝ ⎛ ⎟ ⎟ ⎟ ⎠ ⎞ ⎜ ⎜ ⎜ ⎝ ⎛ − − − 3 2 1)) )) )) )) )) )) )) 4 4 4 4 4 4 3 4 4 4 4 4 4 2 1 b b b b b b b m m m m m m . (2.25) Матрица дискретной дивергенции p p N 3 N R S × ∈ , также как и матрица дискретного ротора |