CC BY
38
Литература
1. The Black Sea Experiment, science and global security. - 1990. - V. 1. - № 3-4. - Р. 323-333.
2. Gordon M.S., Goldhagen P., Rodbell K.P., Zabel T.H., Tang H.H.K., Clem J.M., Bailey P. Measurement of the flux and energy spectrum of cosmic-ray induced neutrons on the ground // IEEE Transactions on Nuclear Science. - 2004. - V. 51. - Is. 6. - P. 3427-3434.
3. Gozani T. The role of neutron based inspection techniques in the post 9/11/01 era // 5th Topical Meeting on Industrial Radiation and Radioisotope Measurement Applications. NIM B. - 2004. - V. 213. - P. 460-463.
4. Wald M.L. Shortage Slows a Program to Detect Nuclear Bombs // The New York Times. - 2009, November 22.
Постоварова Дарья Владимировна - Радиевый институт им. В.Г. Хлопина, младший научный сотрудник,
УДК 517.946:518.12:538.3:538.5
ОБ УЧЕТЕ ИСТОЧНИКОВ ЭЛЕКТРОМАГНИТНОГО ПОЛЯ В СОВМЕСТНОМ МЕТОДЕ КОНЕЧНЫХ И ГРАНИЧНЫХ ЭЛЕМЕНТОВ
И.М. Ступаков, М.М. Корсун, М.Э. Рояк
Рассматриваются два способа учета источников поля в математической модели, используемой при моделировании трехмерных электромагнитных процессов совместным методом конечных и граничных элементов. Один из способов основан на введении в расчетную область поверхности разрыва полного скалярного потенциала, во втором используется представление напряженности магнитного поля в виде суммы градиента неполного потенциала и напряженности поля токовых обмоток.
Ключевые слова: метод конечных и граничных элементов, электромагнитные процессы.
Введение
Одним из эффективных методов моделирования трехмерных нестационарных электромагнитных процессов в технических конструкциях с высоким контрастом проводимостей и магнитных проницаемо-стей является векторный метод конечных элементов. Основной проблемой применения этого метода является плохая сходимость итерационных методов решения получаемых систем линейных алгебраических уравнений (СЛАУ). Поэтому для решения нестационарных задач электромагнетизма, в которых электромагнитное поле формируется за счет вихревых токов, циркулирующих в отдельных объектах с высокой проводимостью, окруженных непроводящей средой, более эффективными являются вычислительные схемы, основанные на совместном использовании векторного и скалярного магнитных потенциалов [1]. В тех случаях, когда непроводящая среда занимает существенную часть расчетной области, можно заметно снизить вычислительные затраты, используя для аппроксимации векторного потенциала векторный метод конечных элементов, а для скалярного потенциала метод граничных элементов [2]. В данной статье будут рассмотрены некоторые особенности учета источников электромагнитного поля при использовании такого подхода.
Математическая модель
Система уравнений Максвелла, используемая для описания нестационарных электромагнитных процессов, при несущественности токов смещения принимает вид:
rot H = JCT +ctE , (1)
5B
rot E =--, (2)
at
div B = 0, (3)
где H и E - напряженность магнитного и электрического полей, JCT - вектор плотностей сторонних токов, ст - удельная проводимость среды, B - индукция магнитного поля, связанная с напряженностью соотношением B = цН, ц - магнитная проницаемость среды.
Среду, в которой изучается электромагнитное поле, разделим на две (возможно, неодносвязных) подобласти QCT и Qo, причем в подобласти Qo удельная проводимость ст = 0. В подобласти QCT индукцию магнитного поля BCT представим через вектор-потенциал A с помощью соотношения
BCT = rot A , (4)
ОБ УЧЕТЕ ИСТОЧНИКОВ ЭЛЕКТРОМАГНИТНОГО ПОЛЯ...
которое автоматически удовлетворяет уравнениям (2) и (3). Уравнение (1) с учетом соотношения (4) преобразуется к виду:
rotf —rot A 1 + ст— = Jст . (5)
U J dt W
Представим в подобласти Q0 напряженность магнитного поля H0 в виде
H0 =-gradU , (6)
где U - полный скалярный потенциал, тогда для выполнения уравнения (3) необходимо, чтобы скалярный магнитный потенциал U в области Qo удовлетворял уравнению
- div (ц grad U ) = 0. (7)
Чтобы решение уравнений (5) и (7) стало решением исходной системы (1)-(3), на границе S между подобластями Qo и QCT нужно выполнить условия непрерывности нормальных составляющих B и тангенциальных составляющих H
B0 n = B°n , H0 х n = H°x n , (8)
где n - любая нормаль к рассматриваемой границе S.
Обратим внимание, что для корректности рассматриваемой модели требуется, чтобы внутри области скалярного потенциала нельзя было бы провести такой замкнутый контур, что интеграл от всех токов, пересекающих натянутую на этот контур поверхность, был бы не равен нулю [1]. Это ограничение создает некоторые дополнительные сложности, связанные как с учетом сторонних, так и с учетом вихревых токов.
Для вихревых токов эта проблема возникает в том случае, если в конструкции имеется замкнутый контур, в котором будет наводиться ток, интегральное значение которого неизвестно. Для разрыва получающегося в таком случае контура обхода ненулевого тока в скалярном потенциале приходится исключать часть расчетной области из Q0 и вводить ее подобласть , задавая в ней нулевую проводимость. Проблема учета сторонних токов связана с тем, что во многих практических задачах известен полный ток в некотором контуре - источнике поля. При этом распределение тока в самом контуре часто либо вообще не представляет интереса, либо известно. Рассмотрим один из способов учета такого источника поля. В подобласти Q0 зададим поверхность Sp разрыва скалярного магнитного потенциала так, чтобы обеспечить отсутствие контура, по которому можно было бы обойти ненулевой ток. Зададим во всех точках поверхности Sp скачок скалярного потенциала, равный известному суммарному току I в контуре:
u\Sp = US+P +1, (9)
где U|sp и U|sp - значения функции U с разных сторон поверхности Sp. Для потенциала U потребуем выполнение дополнительного условия:
- dU
ц —
dn
+ dU
= ц —
SP dn
dU
где
dn
dU и -
SP dn
S+P
(10)
- значения нормальных составляющих H с разных сторон поверхности SP, ц+ и
S+P
ц - значения коэффициента ц с разных сторон поверхности 8Р. Выполнение условия (10) обеспечивает непрерывность нормальных составляющих индукции В на 8Р, непрерывность касательных
составляющих Н на 8Р следует из (9). Некоторым недостатком этого способа учета токового контура является то, что он позволяет простым способом учитывать только поверхностно распределенные токи - в этом случае саму область токового контура можно не включать в расчетную область, задав на поверхности контура краевое условие отсутствия нормальной составляющей магнитной индукции. Если же токи не являются поверхностными и их распределение существенно, то сам контур приходится включить в подобласть .
Рассмотрим другой способ учета источников электромагнитного поля в совместной векторно-скалярной постановке, основанный на использовании неполного скалярного потенциала. Представим в
подобласти напряженность магнитного поля Н0 не в виде (6), а в виде
Н0 = Нс -ягаа^ , (11)
где Hc - напряженность магнитного поля от токов JCT в вакууме (которая может быть вычислена, например, по формуле Био-Савара [1]), а U - неполный скалярный магнитный потенциал. Тогда для определения U необходимо вместо однородного уравнения (7) решать неоднородное уравнение
- div (| grad U ) = - div (|Hc ). (12)
Однако чаще всего в подобласти Q0 магнитная проницаемость совпадает с магнитной проницаемостью вакуума, поэтому в (12) правая часть становится нулевой. В этом случае Hc необходимо учитывать только в уравнениях (8). Заметим, что при необходимости введения в подобласть Q0 непостоянного значения магнитной проницаемости (например, при моделировании конструкций с шихтованным железом, в котором не создаются вихревые токи), можно использовать и более сложную модель в Q0 на основе двух скалярных потенциалов [1].
Совместное использование векторного метода конечных элементов и метода граничных элементов
Аппроксимация по времени уравнения (5) приводит к векторному уравнению
rotp-rot Aj + yA = F , (13)
где коэффициент y и вектор-функция F определяются разностной схемой аппроксимации по времени.
Пусть n - внешняя по отношению к Qa нормаль к границе S между подобластями Q0 и Qa . Тогда, учитывая представление индукции и напряженности через векторный полный и неполный скалярный потенциалы, условия непрерывности (8) принимают вид
— rot A х n
vI
S
f Hc 5U
И H -n- —
v a», s
fc
(( - gradU) х n)
(14)
= (rotA-n)j , (15)
где H = 0 при использовании полного скалярного потенциала.
Будем считать, что граница расчетной области Q либо может быть бесконечно удаленной, либо на ней в качестве краевых условий могут быть заданы краевые условия равенства нулю касательных,
либо нормальных составляющих магнитной индукции. Обозначим через Sх участок границы, на котором задано условие равенства нулю касательных составляющих H , а через S" - участок границы расчетной области Q , на котором задано условие равенства нулю нормальных составляющих B .
Введем следующие обозначения для скалярных произведений, соответствующие интегралам по
объему и по границе: (v, w) =Jv - wdV , (v, w) =Jv - wdT . С учетом введенных обозначений эквива-
у г Г
лентная вариационная постановка для векторного уравнения (13) и условия сопряжения (14) принимает вид [1]:
f —rotA,rotY I + (gradUхn,Y)S + (yA,Y)Q = (F,Y)q +(hc хn,Y) . (16)
v I 'Qa
Уравнение (16) должно выполняться для VY е И0°(Qa), где под H0ot (Qa) понимается пространство вектор-функций Y, определенных на Qa, для которых функция rot Y является суммируемой с
квадратом, при этом касательные всех функций Y должны быть равны нулю на границе Sn n Sa (Sa -граница области Qa).
В том случае, когда коэффициент магнитной проницаемости в области Qo является постоянным, уравнение (7) эквивалентно уравнению Лапласа
- div (grad U) = 0, (17)
тогда, используя представление решения этого уравнения в виде суммы потенциалов, можно получить граничные интегральные уравнения и построить аппроксимацию скалярного потенциала в подобласти Qo с помощью метода граничных элементов [2].
Используя оператор Стеклова-Пуанкаре S, который связывает значения потенциала на границе и потоки через границу, запишем вариационное уравнение соответствующее условию (16)
ОБ УЧЕТЕ ИСТОЧНИКОВ ЭЛЕКТРОМАГНИТНОГО ПОЛЯ.
(и\80,V)0 =(А• п-цИс-п,V) Уу е и\ (£0), (18)
где £0 - граница области ОД, Н (£0) - гильбертово пространство скалярных функций V,
определенных на границе £0, имеющих интегрируемые с квадратом первые производные и равные нулю х 0
на границе £ п £ (более подробно вывод уравнения (18) приводится в [3]).
Построив дискретные аппроксимации уравнений (16) и (18), получим систему линейных уравнений, решение которой позволяет найти искомые потенциалы и и А .
Результаты численного моделирования
На примере решения модельной задачи проиллюстрируем возможности применения рассмотренных методов учета источников тока в совместном методе конечных и граничных элементов. Пусть в расчетной области содержится объект в виде усеченного металлического цилиндра толщиной 5 и радиусом
50 мм с удельной электропроводностью ст = 10 См/м и магнитной проницаемостью ц = 10 . Источником
возбуждения электромагнитного поля является петля, имеющая радиус 25 мм и расположенная соосно с объектом на высоте 10 мм над ним. В начальный момент времени токовый импульс отсутствует, в следующие моменты равен 1 А. Поскольку задача является двумерной в цилиндрических координатах, то результаты трехмерного моделирования будем сравнивать с двумерным решением, которое получено с погрешностью, не превышающей 0,1%.
Будем измерять г-компоненту магнитной индукции В в момент времени / = 50 мкс. Линия измерений расположена на высоте 5 мм над цилиндром вдоль оси у (рис. 1, 2).
Рис. 1. Погрешность расчетов Вх на грубых сетках: - результаты, полученные при непосредственном задании источника в векторном потенциале; - результаты, полученные при учете источника тока
на основе задания поверхности разрыва Бр; -О - результаты, полученные при учете источника тока
на основе метода с неполным потенциалом
На представленных графиках показана погрешность расчетов г-компоненты магнитной индукции В относительно двумерного решения. Анализ результатов показывает, что реализованные методы дают погрешность не более 7% на достаточно грубых сетках, дальнейшее дробление сеток уменьшает погрешность до 1,5%. При этом метод, основанный на использовании неполного потенциала, дает наименьшую погрешность при меньших вычислительных затратах за счет отсутствия неизвестных, требуемых для дискретизации источника.
Bz, %
0
-0,5
-1
-1,5
0 0,01 0,02 0,03 0,04 0,05 У, м
Рис. 2. Погрешность расчетов Вг на подробных сетках: - результаты, полученные при непосредственном задании источника в векторном потенциале; - результаты, полученные при учете источника
тока на основе задания поверхности разрыва Б? ; - результаты, полученные при учете источника
тока на основе метода с неполным потенциалом
Заключение
Анализируя рассмотренные модели и вариационные постановки, заметим, что учет источников тока с использованием неполного скалярного потенциала технологически удобнее при использовании граничных элементов, поскольку в отличие от полного скалярного потенциала не требует введения дополнительных неизвестных в подобласти с граничными элементами на поверхности разрыва скалярного
потенциала. Дополнительные же вычисления поля Нс требуются только на границе подобластей ОД и Ост , что не приводит к существенному увеличению общих вычислительных затрат на решение задачи.
Работа выполнена при поддержке ФЦП «Научные и научно-педагогические кадры инновационной России» на 2009 - 2013 годы.
Литература
1. Соловейчик Ю.Г., Рояк М.Э., Персова М.Г. Метод конечных элементов для решения скалярных и векторных задач: Учеб. пособие. - Новосибирск: Изд-во НГТУ, 2007. - 869 с.
2. Steinbach O. Numerical Approximation Methods for Elliptic Boundary Value Problems. - New York: Springer Science, 2008. - 388 p.
3. Корсун М.М., Ступаков И.М., Рояк М.Э. Об использовании граничных элементов при моделировании электромагнитных процессов с существенным влиянием вихревых токов // Научный вестник Новосибирского ГТУ. - 2010. - № 2 (39). - С. 101-109.
Ступаков Илья Михайлович Корсун Мария Михайловна Рояк Михаил Эммануилович
Новосибирский государственный технический университет, аспирант, [email protected]
Новосибирский государственный технический университет, аспирант, [email protected]
Новосибирский государственный технический университет, доктор технических наук, профессор, [email protected]
