CC BY
18
DOI: 10.5862/JCSTCS.236.3 УДК 519.612
Р.Ю. Бородулин
метод преобразования матриц для решения двумерных электродинамических задач излучения узловым методом
конечных элементов
R.Yu. Borodulin
matrix transform method for solving two-dimensional electrodynamic problems of radiation by nodal finite
element method
Представлен математический аппарат решения задач излучения узловым методом конечных элементов в двумерной области. Получены итоговые выражения для элементов локальной матрицы и методика вывода функциональной зависимости возбуждения точечного источника (для формирования вектора правой части). данным методом выполнен тестовый расчет излучения усовершенствованной модели элементарного электрического диполя в виде плоской рамки с током. достоверность результатов подтверждена высоким совпадением с аналитическим решением задачи излучения элементарного электрического диполя. Представленные результаты могут использоваться при расчетах любых осесимметричных антенн, возбуждаемых точечным источником, в том числе в задачах конструкционного синтеза осесимметричных излучателей различного назначения.
аксиально-симметричное поле излучения; осесимметричная антенна; точечный источник; элементарный электрический диполь; локальная конечно-элементная матрица; функция возбуждения; конструкционный синтез.
The article presents a mathematical apparatus for solving the radiation problems by the nodal finite elements method in a two-dimensional field. we obtain the final expression for the elements of the local matrix and the methodology for deriving the functional dependence of the excitation field of a point source (for forming the right side of the vector). the calculation test performed by this method improved the model of elementary electric dipole in the form of a flat frame with electric current. The reliability of the results confirmed a high coincidence with the analytical solution of the radiation problem of the elementary electric dipole. the results of the article can be used in the calculations of any axially symmetric antenna excited by a point source, including the problems of structural synthesis of axially symmetric emitters for different purposes.
AxiALLY symmetrical radiation FIELD; AxiALLY sYMMETRIc antenna; poINT souRcE; ELEMENTARY ELEcTRic DipoLE; LocAL Finite ELEMENT MATRix; ExciTATioN
function; structural synthesis.
Осесимметричные (аксиально-симметричные) антенны составляют большой класс излучателей. Широкое применение таких антенн в технике радиосвязи выявило ряд проблем, касающихся их массогабарит-ных показателей, решаемых при помощи теории конструкционного синтеза (кс), требующей быстроты электродинамического анализа синтезируемых излучателей [1].
Оптимизационные методы, применяемые в кс, имеют большое число итераций, и чем меньшее время будет потрачено на одну такую итерацию, тем быстрее будет найдена оптимальная конструкция антенны.
любую аксиально-симметричную антенну можно представить в двумерной области в виде половины ее сечения. это позволяет существенно экономить ресур-
4
сы ЭВМ, затрачиваемые на ее анализ, и вместо векторных уравнений использовать скалярные, имеющие меньшую размерность.
Существуют проблемы решения задач излучения узловым методом конечных элементов (МКЭ). Они заключаются в некорректности работы методик, разработанных для стационарных аксиально-симметричных полей применительно к уравнению Гельмгольца, описывающему распространяющуюся аксиально-симметричную волну, т. к. решение, получаемое МКЭ, не совпадает с решением, получаемым строгими методами для элементарного магнитного диполя (ЭМД) (рис. 1 а, б) [2-4].
Данный факт приводит к необходимости разработки метода повышения точности решения неоднородного уравнения гельмгольца в двумерной области для аксиально-симметричных волн:
дг
1 _д_
г цг дг
(Е)
д_
дг
(
1
Г д( Е ) ^
г Ц г
дг
/у
к 2р
^^ (гЕф) = -Мэ,
(1)
где г - расстояние от оси симметрии до точки наблюдения; цг — относительная магнитная проницаемость среды распространения; Е — азимутальная компонента
электрического поля; егк — относительная комплексная диэлектрическая проницаемость; Мэ — источник, порождающий волну.
Вывод «матрицы жесткости» конечного элемента для решения осесимметричных задач
Аналитическое решение задачи излучения элементарного магнитного диполя (ЭМД) для азимутальной компоненты вектора электрического поля в ближней зоне выглядит как
^ = Ф "4^ 7 + к 1 '
'ЯП 0, (2)
где тт — магнитный момент диполя; ш — круговая частота; 0 — меридиональный угол [3].
Преобразуем «матрицу жесткости» конечного элемента (термин взят из теории упругости) к виду, соответствующему уравнению (1):
(
[К ] = Ч 1
дН дИ дИ,
дг дг
дг
дН
дг
дИ
дг
дН
дг
+
а)
б)
К^Е^В/м
0.1 02 03 0.4 0.5 0.6 0.7 0.8 0.
Рис. 1. Результаты расчета действительной компоненты электрического поля круглой рамки
радиусом а = 0,05 м, имитирующей ЭМД, классическим двумерным МКЭ: а — в направлении максимума излучения; б — в виде амплитуды Е , распределенной в расчетной области
дИ1 дИ2 дИ3
д7 д7 87
дИ, Л
д7
ЭИ2
д7
дИ3
д7 /
К*. (3)
К слагаемым подынтегральной функции
добавим
к2 6„
■(гЕ ) из (1). Функция формы
треугольного элемента Н1 (для узлов с локальными номерами , = 1, 2, 3), показанная в данном уравнении, может описываться полиномами различных степеней. Так, если функция формы содержит полиномы первой степени вида р = сс с с2хс с с3х2, где в нашем случае хс = г, х2 = 7, элемент принято называть Рс-элементом. Если функция формы содержит полиномы ввторой степени вида р = сс с с2хс с с3х2 с с4хс2 с с5хсх2 с с6х22, такоИ элемент принято называть Р2-элемен=ом.
Если рассматривать Рс-элементы, то полученная нами функция формы И. примет следующий вид:
Н1 = 2юГ ^^ - Г3^ с (*2 - ^ с (г3 - г2=] , (4)
И2 = 24[(г,т.х - гс 13) с (^3 - с (гс - г,)^],(5)
И3 = 24 [(г1^2 - г2с - 12)г с (г2 - гск],(6)
то есть будет являться по сути матрицей базисных функций единичного элемента. Здесь (г- 7.-,) — координаты соответств4ю-щих вершин элемента.
Данная функция формы удовлетворяет условию
И, (г,-, 7,-) = 8,, (7)
то есть
IИ1 = 1.
(8)
Здесь 8, — дельта Кронекера. Она будет равна
8 I1' - = 1
8 =|о, -* -
(9)
После подстановки выражений (4)—(6)
в уравнение (3) получим «матрицу жесткости» элемента Ке вида
[К ] =
И ^2 ^3
^22 ^23
^31 2 3
(10)
элементы которой получены путем дифференцирования И1 :
= 4Ю
г1(72 73) / \2
-— с (г3 - г,)2
1
к12 = - ■
г2(72 - 73)(73 - 71)
4 Ю
с (г3 - г2)(г1 - г3)
^ = 41Ю
т. - 73X7 - ^^^
с (г3 - г2)(г2 - г1)
1
к21 = - ■
г1(73 - 71)(72 - ^^^
4 Ю
с (г1 - г3)(г3 - г?)
^22 = 4Ю
г2(73 - ¿1)2
с (г - г?)2
к =
23 _
_1_ 4 Ю
г3(73 - 71)(71 - ^
с (г - г3)(г2 - г1)
к л = . .
г1(71 - 72)(72 - 73)
к =
32
4 Ю [_ гс
с (г2 - г1)(г3 - г2)
йг2(7 - 72)(73 - ^
4 Ю [_ гс
с (г2 - г1)(г1 - г3)
к33 = 4 Ю
г3(71 - 72)2
с (г2 - г1)2
(11)
(12)
(13)
(14)
, (15)
(16)
(17)
(18)
(19)
с
,=1
где
1 г1 г"
А = 0,5 1 г2 г2 (20)
1 г3 г _
является площадью треугольного элемента, выраженной через узловые координаты его вершин.
Формирование матрицы H и векторов F и G производится аналогично классическому узловому МКЭ для декартовых координат.
Как ни странно, выражения (11)—(19) именно в такой интерпретации в литературе не встречались, поэтому будем считать, что они получены нами впервые путем прямой подстановки функций формы в исходное уравнение гельмгольца в цилиндрических координатах. Именно такое представление матрицы жесткости приводит к правильному снижению комплексных амплитуд поля при удалении от оси симметрии и практически полному совпадению фазы волны с фазой точного решения, что в свою очередь позволяет существенно снизить относительную погрешность действительных и мнимых составляющих Еф, что будет показано ниже.
Методика вывода функциональной зависимости возбуждения точечного источника для решения осесимметричных задач излучения
Ненулевые элементы правой части Ь решаемого СЛАУ уже рассматривались ранее, где они были равны -¿8/тт ф 0. Таким образом, за счет априорного задания данной величины, можно получить начальное приближение.
Представление функции возбуждения начнем с рассмотрения неоднородного уравнения Гельмгольца (1). Для данного уравнения правая часть (функция возбуждения) записывается как
М3 = -1шца]эст +
(21)
+-£га<Му/э ст - го/ м ст,
_ 1ШВак
• э ст • м ст
где / . , / ' — мгновенные значения век-
торов объемной плотности электрического и магнитного тока; еак, ца — абсолютные значения комплексной диэлектрической и магнитной проницаемостей среды. Так как МКЭ решается скалярное уравнение, то векторные слагаемые лишаются смысла.
Показанное выше представление функции возбуждения в виде точки нецелесообразно, т. к. в силу зависимости составляющих вектора F от величины площади элемента А при уменьшении последней в ходе учащения сетки разбиения будут уменьшаться и полученные комплексные амплитуды.
Тогда, в силу сложности вывода строгого соотношения для правой части исходного скалярного уравнения Гельмгольца, выражение (21) можно переписать с точностью до постоянного множителя у(аИ) как
Ммкэ = -шЦоЦг /эст + -Ш— Егаа^у/эст =
1 шеак (22)
= - шцоц^тт; у(а/о,
где у(ак) — скалярная величина, пропорциональная плотности стороннего тока (зависящая от положения источника относительно оси симметрии (радиуса рамки а) и ширины полоски к, образующей рамку).
Таким образом, в двумерном аксиально-симметричном представлении рамка будет представлена нитью постоянного тока, параллельной оси симметрии, длина которой равна ширине полоски к (рис. 2).
Такое представление модели возбуждения позволяет избежать зависимости полученных комплексных амплитуд поля от размера элементов.
на практике задание токовой нити в СЛАУ [К + И] • x = Р + О программируется в векторе О в виде дополнительной границы с ГУ вида (6), где к = 0, g = - шц0цг У(ак).
В выражении (22) величина функции возбуждения заменяется на Iту(а, к), так как сторонний источник в задаче представляется в виде бесконечно тонкой плоской круглой рамки. Из анализа выражения (2) видно, что величина Ет пропорциональна магнитному моменту, т. е. прямо пропорциональна стороннему току /тт, радиусу рамки и обратно пропорциональна длине волны 1. Тогда выражение (22) можно
г = 0
Рис. 2. Усовершенствованная модель ЭМД в виде плоской рамки с током
с точностью до постоянного множителя у(а, И) можно переписать как
Ммкэ = - у (а, И), (23)
где с — фазовая скорость волны в среде рас-
где А — относительная погрешность; Е^ (г, ^ = 1) — рассчитанные МКЭ комплексные амплитуды вектора электрического поля; Е®н (г, ^ = 1) — комплексные амплитуды поля, полученные аналитически согласно выражению (2).
Результаты поиска представлены в таблице.
Необходимо оговориться, что для поиска минимума функционала (24) можно использовать любой метод оптимизации. Однако, учитывая, что аппроксимация множества полученных точек некоторой непрерывной функцией все равно приведет к некоторой погрешности, в данной ситуации метод перебора вполне приемлем.
Полученные результаты лишний раз подтверждают тот факт, что с увеличением радиуса рамки при фиксированной толщи-
пространения (в нашем случае — скорость света).
Поиск значения величины j(a, h) произведем методом перебора с допуском 2 % погрешности, путем минимизации следующего функционала:
(24)
не полоски плотность тока, распределенная по полоске, определяемая выражением im y(a, h), закономерно убывает.
Интерполяцию дискретных значений функции y(a, h) для нахождения непрерывной зависимости произведем полиномом Лагранжа.
Аппроксимирующий полином Лагран-жа, пригодный для практических алгоритмов, получаемый путем интерполяции Y(a, h = const), в зависимости от a имеет вид
L(a / X) = £y. П a -(a . " Si Д1...(a / X), - (a / X),
j Ф,
На рис. 3 видно, что полученную из таблицы кривую достаточно удобно аппроксимировать полиномом вида y = C.e x + C2x_1. Итоговое выражение для
Д = mean
Re(Е?Пr, z = 1)) - Re(£;н(r, z = 1))
Re( Ефн( r, z = 1))
^ mm,
4
Результаты минимизации функционала (24) методом перебора при фиксированной к = 0,02 м
а, м 0,01 0,015 0,02 0,025 0,03 0,035 0,04 0,045 0,05 0,055 0,06
у(а, к) 940 615 480 375 330 283 252 229,5 215 198 186
Д, % 1,9 1,9 1,9 2,0 2,0 2,0 2,1 2,2 2,2 2,3 2,4
1
Рис. 3. Зависимость £1(а) и у от а/Х
трех членов, пригодное для расчета величины у(а, к = 0,02), полученное методом наименьших квадратов, приняло следующий вид:
у(а, к = 0,02) « Ь1(а) =
= 30,763 ехр(-а) + 8,977(а"1), где а — радиус витка.
(25)
На рис. 4 и 5 показаны результаты расчета поля излучения рамки радиусом а = 0,03 м с применением модифицированной матрицы жесткости (3) для реальной и мнимой частей комплексной амплитуды в сравнении с результатами, полученными согласно выражению (2). При таком радиусе рамки относительная погрешность МКЭ составляет 2 %.
Из представленных рисунков видно практически полное совпадение полученного при помощи выведенной функции у (а, к = 0,02) МКЭ-решения в сравнении со значениями комплексных амплитуд, полученными аналитически. Если увеличивать радиус рамки, погрешность может достигнуть значения 2,3 % и более, что является вполне закономерным.
Выведенные выражения для элементов «матрицы жесткости» позволяют производить электродинамический анализ осесим-метричных излучателей с высокой точностью.
Таким образом, полученная нами методика вывода функциональной зависимости возбуждения точечного источника для ре-
Рис. 4. Результаты численного расчета действительной составляющей комплексной амплитуды вектора Е1расч в сравнении с Е™, полученными строгим методом рамки радиусом 0,03 м
Рис. 5. Результаты численного расчета мнимой составляющей комплексной амплитуды вектора E(расч в сравнении с E(í¡н, полученными строгим методом рамки радиусом 0,03 м
шения задач излучения в цилиндрических координатах будет содержать следующие этапы:
фиксация ширины полоски h, желательно симметрично точке возбуждения;
получение дискретных значений функции y(a, h = const) в зависимости от радиусов витков;
нахождение C1 и C2 полинома Лагран-жа вида y = C1ex + C2x- для получения непрерывной функциональной зависимости функции возбуждения от радиуса витка путем минимизации функционала (24) с оценкой относительной погрешности любым оптимизационным методом;
оценка полученного результата путем
сравнения МКЭ-решения с решением, полученным любым строгим методом.
Выведенные выражения для элементов «матрицы жесткости» позволяют производить электродинамический анализ (конструкционный синтез) осесимметричных излучателей с высокой точностью.
Описанная методика получения функции возбуждения позволяет производить корректное моделирование источника возбуждения осесимметричных тел. Ее можно легко расширить на предмет анализа осе-симметричных антенн, существенно уменьшив при этом размерность задачи и получив выигрыш во времени.
1. Сосунов Б.В., Бородулин Р.Ю. Конструкционный синтез элементов фазированных антенных решеток // Научно-технические ведомости СПбГПУ. Информатика. Телекоммуникации. Управление. СПб., 2013. № 2(169). С. 47-54.
2. Kwon Young W. The finite element method using MATLAB (The mechanical engineering
СПИСОК ЛИТЕРАТУРЫ
series). New York: CRC Press LLC, 1997.
3. Jianming Jin. The finite element method in electromagnetics. 2nd ed. New York: John Willey and Sons Inc., 2002.
4. Иикольский В.В., Иикольскал Т.И. Электродинамика и распространение радиоволн. М.: Наука, 1989.
REFERENCES
1. Sosunov B.V., Borodulin R.Yu.
Konstruktsionnyy sintez elementov fazirovannykh antennykh reshetok [Constructural synthesis of element of a fased array of antennas]. Nauchno-tekhnicheskiye vedomosti SPbGPU.
[St. Petersburg State Polytechnical University Journal. Computer Science. Telecommunication and Control]. St. Petersburg, 2013, No. 2(169), Pp. 47-54. (rus)
2. Kwon Young W. The finite element method
Informatika. Telekommunikatsii. Upravleniye using MATLAB. The mechanical engineering series,
New York: CRC Press LLC, 1997.
3. Jianming Jin. The finite element method in electromagnetics, 2nd ed., New York: John Willey and Sons Inc., 2002.
4. Nikolskiy V.V., Nikolskaya T.I. Elektrodinamika i rasprostraneniye radiovoln [Electrodynamics and distribution of radio waves]. Moscow: Nauka Publ., 1989. (rus)
БОРОДУЛИН Роман Юрьевич — докторант Военной академии связи имени С.М. Буденного, кандидат технических наук.
194064, Россия, Санкт-Петербург, Тихорецкий пр., д. 3. E-mail: borodulroman@yandex.ru
BORODULIN Roman Yu. Military academy of communications named after S.M. Budyonny. 194064, Tikhoretsky Ave. 3, St. Petersburg, Russia. E-mail: borodulroman@yandex.ru
© Санкт-Петербургский политехнический университет Петра Великого, 2016
