CC BY
8
2020 Математика и механика № 64
УДК 519.63 М8С 76W05
Б01 10.17223/19988621/64/2
И.В. Бычин, А.В. Гореликов, А.В. Ряховский
ЧИСЛЕННОЕ РЕШЕНИЕ НАЧАЛЬНО-КРАЕВОЙ ЗАДАЧИ С ВАКУУМНЫМИ ГРАНИЧНЫМИ УСЛОВИЯМИ ДЛЯ УРАВНЕНИЯ ИНДУКЦИИ МАГНИТНОГО ПОЛЯ В ШАРЕ1
Разработан алгоритм численного решения начально-краевой задачи с вакуумными граничными условиями для уравнения индукции магнитного поля в шаре. При дискретизации использован метод контрольного объема и модификация ЕБТО-метода, учитывающая специфику рассматриваемой задачи. Внешняя задача Неймана для уравнения Лапласа на потенциал магнитного поля в вакууме решалась с использованием преобразования обратных радиусов. Представлены результаты тестирования.
Ключевые слова: уравнение индукции, вакуумные граничные условия, численное решение.
Необходимость решения уравнения магнитной индукции с вакуумными граничными условиями зачастую возникает при численном решении задач магнитной гидродинамики, например при математическом моделировании явления гидромагнитного динамо [1-5]. При постановке задач с вакуумными граничными условиями магнитное поле считается всюду непрерывным, а в вакууме - потенциальным и соленоидальным [6]. В коммерческих программных кодах обычно используются упрощенные так называемые «псевдо-вакуумные» граничные условия, при которых на границе области полагается равной нулю тангенциальная составляющая магнитного поля [7, 8]. Решения, полученные с использованием псевдо-вакуумных условий, обычно корректно воспроизводят плотность тока, однако структура магнитного поля значительно отличается от решений, рассчитанных при вакуумных условиях [8]. Поэтому в задачах моделирования гидромагнитного динамо представляется целесообразным использование более реалистичных вакуумных граничных условий.
При численном решении уравнения индукции, возникает проблема получения соленоидального магнитного поля. В настоящее время наиболее распространенными методами вычислительной магнитной гидродинамики являются псевдоспектральные [9, 10] и конечно-разностные [11, 12] методы, а также метод контрольного объема [10, 13]. В рамках данных методов существуют различные подходы для обеспечения бездивергентности магнитного поля [14]: использование векторного потенциала; метод искусственного скалярного потенциала [15]; представление дискретного аналога уравнения индукции в такой форме, чтобы его решение автоматически удовлетворяло сеточному уравнению неразрывности [16-18]; метод Пауэлла, который широко применяется в астрофизических приложениях и основан на записи уравнений в форме, допускающей существование магнитных монополей [19]. Перечисленные методы обладают своими преимуще-
1 Работа выполнена при поддержке Программы ФНИ государственных академий наук на 2013-2020 гг., проект № 0065-2019-0021.
ствами и недостатками. Например, коррекция магнитного поля в методе скалярного потенциала нелокальна, и поэтому локальные ошибки в дивергенции магнитного поля могут мгновенно распространяться на всю область [20]. Другие методы, например метод Пауэлла, не являются, строго говоря, консервативными.
Несмотря на достаточно большое число работ в этой области, разработка эффективных алгоритмов решения уравнения индукции с вакуумными граничными условиями остается актуальной задачей вычислительной магнитной гидродинамики. Целью настоящей работы является разработка и апробация алгоритма численного решения начально-краевой задачи с вакуумными граничными условиями для уравнения индукции магнитного поля в шаре на основе модификации разностных методов вычислительной гидродинамики и электродинамики для использования в ортогональных криволинейных координатах.
1. Постановка задачи и математическая модель
Пусть проводящая жидкость заполняет шар G с центром в начале координат и радиусом a, который находится в неподвижной диэлектрической среде. Предполагается, что удельная электрическая проводимость жидкости с = const. В рамках модели МГД [21-23] рассматривается задача о нахождении индукции магнитного поля B, при заданной скорости жидкости u. Задача решается в сферических координатах r, 9, ф. Полагается, что вне шара магнитное поле потенциально, и его скалярный потенциал у является регулярной на бесконечности гармонической функцией. На границе проводника и диэлектрика задаются условия сопряжения (так называемые вакуумные граничные условия), которые заключаются в требовании непрерывности магнитного поля и равенства нулю нормальной составляющей
c
плотности тока j = — rot B [21]. Таким образом, начально-краевая задача, для которой в данной работе строится процедура численного решения, имеет вид r < a, t > 0:
dB
--rot (u X B) + v mrotrot B = 0; (1)
dt
r > a, t > 0:
r = a, t > 0:
V2y = 0 (у - регулярна при r ^ +<»); (2)
r < a, t = 0:
B = Vy, (3)
(rot B)r = 0; (4)
B = b, где div b = 0 . (5)
c2
Здесь t - время, c - скорость света, vm =--коэффициент магнитной вязкости.
4яст
Из (1) следует, что равенство нулю дивергенции начального поля b обеспечивает соленоидальность магнитного поля B и в последующие моменты времени.
2. Численный метод решения уравнения индукции
В данной работе для дискретизации уравнения индукции использован метод контрольного объема [24-26] и модификация метода FDTD (Finite Difference Time Domain) [17, 27, 28] в ортогональных криволинейных координатах, учитывающая специфику рассматриваемой задачи (1) - (5).
Дискретизация расчетной области
Пусть {ха} - произвольные ортогональные координаты, {еа} - соответствующий правый ортонормированный базис, {На} - коэффициенты Ламе (здесь и далее, нижние греческие индексы принимают значения 1, 2, 3 и определяют номер координаты или соответствующей компоненты). Предполагается, что расчетная область G - ограниченная область в пространстве, которую можно разбить системой координатных поверхностей
{ха,га = С0П^ а = 1,2А 'а = Ъ — Па : Ха,га-1 < }
на непересекающиеся подобласти - контрольные объемы
Айг3 = {(x1, Х2 , x3): Ха 6 (ха,га-1, Ха,га ), а = 1 2, 3}, га = 2,^,па.
В центре каждого контрольного объема содержится внутренняя расчетная точка P ■ 3. В целях сокращения записи, в дальнейшем набор индексов /ь i2, i3 обозначается как (i). Сеточные значения компонент магнитного поля рассчитываются в точках на гранях контрольных объемов, т.е. используется дискретизация на смещённых расчетных сетках [24-28]. Для описания процедуры дискретизации на смещённых сетках удобно использовать операторы сдвига на наборах индексов
К (i) - hi ('!, i2, /3) = hi (О, hi (i2), hi ('з) , (6)
действие которых определяется по правилу:
hl Op ) = 'в +i 8ар , а = 1,2,3, l=-l1, (7)
где 5ар - символ Кронекера, а индекс с определяет направление и величину сдвига по индексу /р. При этом полагается, что расчетные точки с одним полуцелым индексом Phl/2 находятся в центрах граней Shl/2 ) контрольного объема D^;
ha (/) ha (/)
5S,a/2) - площадь соответствующей грани контрольного объема; граница кон-
ha (/)
3
трольного объема dD^ = М S,±1/2(. Точки с двумя полуцелыми индексами
а=1
P*^/2(г) находятся в центрах рёбер /у, ^/2(/) = /2(г) П Shl/2(/) (Y Ф а Ф e)
контрольного объема D(f) расположенных вдоль координатных линий xY; 5/ ,x/2,a/2() - длина соответствующего ребра контрольного объема D^, и
У, he Ца (i) 5ShS/2(г) = lY, h±1/2hl/2(г) М lp, h±1/2hi/2(0 '
Дискретный аналог уравнения индукции
Уравнение (1) является следствием закона электромагнитной индукции Фарадея:
dB
— = - с rot E, (8)
dt
где E - напряженность электрического поля, которая в модели МГД имеет вид [22]
E = -^(u х B) + — rot B . (9)
с 4яст
Для получения дискретного аналога уравнение (8) умножается на eadSadt (dSa -элемент площади на грани S,a/2( )), и затем интегрируется по поверхности
ha (i)
S,a /2 и по промежутку времени [t0, t]. При вычислении интеграла по поверхно-
ha (i)
сти S,a/2(.) используется теорема Стокса, а при интегрировании по времени -
ha (i)
полностью неявная схема. В результате, дискретный аналог уравнения (8) для компоненты Ba можно записать в виде
1 / \ 3
St (() - Сго-)) 5V(i) = -СХ ^ Ав ( £^,2(() . (10)
Здесь B ,a/2(.) = Ba(P,G/2( t) - сеточные значения Ba на текущем временном слое;
a,ha (i) ha (i)
B0 /2 = Ba(P /2 t0) - сеточные значения Ba на предыдущем временном слое;
a,^a (i) ha (i)
E hX/2,a/2 ) = EY(P /2 /2 t) - сеточные значения компоненты напряженности
У, "р ha (l) ' "р ha (i)
электрического поля в точках на ребрах контрольного объема (a Ф в Ф у); 5t = t - t0 - шаг по времени; eapy - символ Леви-Чивиты; конечные разности в правой части (10) определяются по формуле
АР (5/уEy /2(i) = EY, hl'h/2(i-) 5/y, /2(i) - Ey, tf'h/2(0 Sy, V/2hg/2(0. (11)
Анализ уравнения (10) для задачи (1) - (5) показывает, что сеточные значения (rot B)r в точках на границе расчетной области не входят в уравнение (10) для приграничных расчетных точек Pni ^ . Следовательно, условие (4) равенства нулю
радиальной составляющей ротора магнитного поля на границе шара G не накладывает каких-либо дополнительных ограничений на вид дискретного аналога уравнения индукции для приграничных расчетных точек.
Дискретный аналог уравнения неразрывности
Сеточные значения div B определяются во внутренних расчетных точках P.) контрольных объемов D^) по формуле
DIV B\P = —Y (( .SS.^,, - B . V (12)
P) SVn Л a,>£2(0 Ц[2(1) a,ha1/2(i) ha (i) J v '
(i) a=1
где 8V(i) - объем D({) и учтено направление внешней нормали n/2 ) =aea на
ha (i)
гранях контрольного объема D(i). Пусть D(i) - внутренний контрольный объем (т.е. его граница не пересекается с границей расчетной области G) и пусть в момент
времени /0 сеточные значения компонент магнитного поля Аст/2(,) удовлетворяют дискретному аналогу уравнения неразрывности Б1У В0 = 0 . Тогда подста-
|ре)
новка В /2 )/2) из сеточного уравнения индукции (10) в правую часть (12)
а, На (1) Ни. (1)
дает
сЫ
Б!У В
^(г) а,р,у=1
X 8ару(др(5/уЕу)к1,2(1)-Др(8/уЕу)й-1/2(г.)), (13)
где с учетом (11) правая часть имеет вид
Дв (8/У ^(Г, "ДР К еу)^-1/2(1) = (8/у ЕУ ^ - (8/у ЕУ ^ —20) - (14)
- ^ ЕУ ^ ^^ + ^ ЕУ ^ ^-"2й-2(1).
Из определения (7), следует, что при а Ф в операторы сдвига Н? 2 и Щ 2 (х, сте{-1,1}) коммутируют, следовательно, выражение (14) симметрично по индексам а и в, и его свертка с символом Леви-Чивиты в правой части (13) равна нулю, т.е.
Б1У В|р = 0. (15)
Таким образом, решение дискретного аналога уравнения индукции (10) удовлетворяет уравнению неразрывности (15) для внутренних контрольных объемов на текущем временном слое, при условии, что оно соленоидально на предыдущем временном слое. Для контрольных объемов, примыкающих к границе О, уравнение неразрывности (15) можно использовать для нахождения нормальной компоненты магнитного поля на границе расчетной области.
Аппроксимация напряженности электрического поля на ребрах контрольного объема
Из (9) следует, что компоненты напряженности электрического поля в ортогональных криволинейных координатах вычисляются по формуле
Еу=-1 X М«тНрВр-^т^Щр) ). аб)
У с нД т Р Р Нт дХт ) У '
Для удобства записи и аппроксимации конечно-разностных выражений в (10) вводятся новые зависимые переменные:
Фр = Нр Вр, р = 1, 2, 3. (17)
По аналогии с вычислительной гидродинамикой, можно определить суммарные (конвективные и диффузионные) потоки [24] переменной Фр по направлению хт в
точках Р х/2?ч на ребрах I „/2 ?ч контрольного объема:
НХ'Н/2(о ^ ^ у, НХ Н? (о ^
I ЧГ V дфр
JYJ.хп.„„... = -^|итФ„- т р
утрН/2(1) нД т р Нт дхт Тогда, конечные разности (11) можно записать в виде
,х/2н?/2
Н?'2(1)
Р(1)
-СДР ('2(0 = X £ухр АР•тр |*2/2(0 ; (19)
т,р=1
ГДе АР^утр\ь1 /2(0 = ^утр |й1/2й^/2(0 - •/™|,-1/2,а/2,.л • (20)
Уравнение индукции (10) в терминах переменных Фр принимает вид
11 \ 3 3
§(/2(0 - Ф^2/2(0 /2(0 = Рх ^ х 81-р АР^ 1*2/2(г) ' (21)
в,7=1 т,р=1 /2(;)
где А,2/2 = а ( ) • (22)
(г) Н 2/2
а, (0
Вычисление сумм в правой части (21) с учетом (20) позволяет записать дискретный аналог уравнения индукции в форме
§7(/20) -Фа,^2/2(0 5)/2(г) =ла[Фа,*2/2(0 , Фа,*в±1*2/2(0 , ^й^1*2/2(г) ] +
+ [Фр,*±1/2(0, Фр,*2*±1/2(0, ФуЛ±1/2,) , Фу,йа2й±1/2,) ] ( а ^ р ^ у ), (23)
где
Ла [Фа,й2/2(0 , Фа,й±1*2/2(0 , Фа,й±1й2/2(0 ] = ^ 1*-1/2й2'2(0 '2(0
+ ^РУа |*-1/2й2/2(г) ^вУа |*1/2*2 '2(гу (24)
^а[Фр,й±1/2(г) , ФР,^1/2(0 , Фу,й±1/2(г) , Фу,й>±1/2(0 ] = УУаРр/2й2'2(0 ~
-^УаР 1йр1/2й2/2(0 +'/Рау 1лф/2^2/2(0 - ^Рау |й-1/2*2'2(0 • (25)
Различные схемы аппроксимации задач конвекции - диффузии подробно рассмотрены в [24-26]. В данной работе на этапе отладки и тестирования нового кода для моделирования МГД-течений использована схема со степенным законом, которая проста в реализации, устойчива и позволяет получать результаты, точность которых сопоставима с численными схемами второго порядка [29-31] при относительно небольших значениях числа Пекле.
Схему со степенным законом можно записать с использованием функции
2(х,у) = (х/у)©(g(х/у)) + х®(х), (26)
где g (х ) = (1 - 0,11X )5, ©( х) = {0
0, х < 0, х > 0.
Аргументами 2(х, у) являются величины, определяющие интенсивность конвекции и диффузионную проводимость в расчетных точках на ребрах контрольных объемов.
Суммарные потоки в конечно-разностном операторе Ла (24) аппроксимируются по формулам
*/2Й2 /2(О = Е'Ра */2*2 /2(0 Фа,й2/2(0 + +Х ЯУРа /2й2/2(г) (Фа,*2/2(0 - Фа,*Х*2/2(0 ), (27)
где
К,
ир5/у
уР<* Н/2н?/2(1) н
А„
ут 8/т
нГн-ю уРа 'Н^2«" на8/р
(28)
,х/2на/2
(1)
(8/ х/2,а/2( ) - длина дуги координатной линии хв, между точками РН?/2(,) и
\ Р, Нр На (1) На (1)
Нх/2Н?/2(1)'
(29)
Ян /2(1));
°УРа 1нрХ /2Н?/2(1) = 2 (-Х ЕуРа, ^УРа )1 Н
Формулы для аппроксимации потоков в (25):
•«РН/2Н?/2(1) = ^Р Н/2Н?/2(1) фв,нх/2(1) + Н/2Н?/2(1) (ФР,Н?НХ/2(1) - ^«Гю^ (30) где
К1
ма8/у
уар |
нХ/2н?/2(1) нв
д„
V т8/у
УаР ы/2н?/2(1) н 8г
НрХ/2Н?/2(1) 15 " НР8/а
(31)
%х/2Н?/2(1)
(8/а,НХ/2Н?
2(1)
- длина дуги координатной линии ха, между точками Р
Х /2
Р
?нХ/2,
(1)
);
= 2 (^ар, ^уар )| Н
УаР \НХ/2Н? 12(1) К уар'^уар; Ц/2Н?/2(1) ' С учетом (27) - (32) выражения для Ла и приобретают вид
Ла[Фа,Н?/2(1) , Фа,Н±1Н?/2(1) , Фа,Н±'Н?/2(1) ] =
(1)
(32)
= (°УРа 1н1/2н?/2(1) + аУРа 1нр1/2н?/2(1) + аРУа 1н}/2Н?/2(1) + ЯРУа IН-1/2Н?/2(1))Фа,Н?/2(г) +
+(-Я
УРа 1н1/2н?/2(1) + ЕУРа 1н-1/2Н?/2(1) 1н}/2Н?/2(1) + 1н71/2н?/2(1 ))Фа,Н?/2(г)
+ ЯУР«Нр/2Н?/2(1) Фа,н1н?/2(0 + °УР« 1нр1/2Н?/2(1) Фа,Нр'Н?/2(1) +
+аРУа 1н}/2Н?/2(1) Фа,н}н?/2(1) + I;-1
ф
2(1) а,
2(1)
(33)
2(1 ) :
, Ф„
2(l), Фу,Н±1/2(1-) , Ф,Н? Нл±1/2(1)-
= Е)/аР 1н1/2н?/2(0 фР,н1/2(1) + ""у-РЦ/2н?/2(1) (ф|
■"в/2(0 ФР,н1/2(1) )
Е'аР 1нр1/2н?/2(1) Фр,Нр1/2(1) аУ«Р 1нр1/2н?/2(1) (Фр,Н?Н-1/2(1) Фр,Нр1/2(1) } + + ^1н;/2н?/2(0 Фу,Н{/2(0 + ^Раг 1н;/2Н?/2(0 (Фу,Н?Н}/2( 1) "фу,н;/2(0 ) -
РРаУ 1н-1/2Н?/2(0 Фу,Н-1/2(1) ЯР<*У 1л-
1/2Н?/2(1) (Фу,Н?л-1/2(1) Фу,л-1/2(1) ^
где а Ф в Ф У.
(34)
и
3. Численное решение начально-краевой задачи (1) - (5)
Внешняя краевая задача на потенциал
Если на границе шара G найдена нормальная компонента магнитного поля Br, то потенциал у в каждый момент времени 7 > 0 в области г > а является регулярным на бесконечности решением внешней задачи Неймана для уравнения Лапласа:
г > а:
V2y = 0; (35)
r = a:
дш=Br (36)
dr
С помощью преобразования обратных радиусов [32]:
rr' = a2, y(r, 6, ф) = r 'ю(г' , 6, ф), (37)
внешнюю задачу Неймана (35), (36) можно свести к решению внутренней краевой задачи для уравнения Лапласа на ю с граничными условиями третьего рода:
Г < a : V2ro = 0; (38)
r' = a : a —+ o=-Br. (39)
dr' r
Краевая задача (38) - (39) численно решалась методом контрольного объема в сферических координатах на той же расчетной сетке, что и уравнение индукции (23). Сеточные значения ю вычислялись во внутренних расчетных точках P( ,) контрольных объемов D(j) и в точках Р,ц2,Л е dG. После этого, по (37) определялся
( О
потенциал у в точках P1/2 е dG, а затем тангенциальные компоненты магнит-
«1 (0
ного поля
1 дш 1 дш r = a : B6= —-, Вф =--- (40)
r 56 т r sin 6 дф
в точках P¥/2й;/2(0 edG (Р ф 1).
Поскольку тангенциальные компоненты магнитного поля на границе шара (40) находятся, как компоненты Vy, то можно показать, что сеточные значения нормальной составляющей ротора магнитного поля (rot B)r в точках Pa/2 х/2,i/2,., edG (у ф в ф 1) автоматически становятся равными нулю, тем са-
«Y «в «1 (i)
мым обеспечивая выполнение условия (4).
Алгоритм решения начально-краевой задачи (1) - (5)
Для удобства записи алгоритма сеточные значения B = Brer+B0ee + Вфеф в точках на границе шара G обозначаются f (т.е. f = B| ). Пусть в момент времени tm-1
с заданной точностью найдено численное решение B(m-1) начально-краевой задачи (1) - (5) и пусть нижний индекс к нумерует внутренние итерации, которые необходимо сделать, чтобы получить решение на следующем временном слое tm.
В качестве нулевого приближения для B на временном слое tm используется B(m-1), т.е. при к = 0 B0 = B(m-1). В данной работе предложен и реализован следующий алгоритм получения численного решения начально-краевой задачи (1) - (5):
Шаг 1: Из сеточного уравнения неразрывности (15) для приграничных контрольных объемов D - з, определяется нормальная (радиальная) компонента
магнитного поля на границе G, т.е. вычисляется f, к+1.
Шаг 2: Решается третья внутренняя краевая задача (38), (39) на потенциал ю, для которой Br\r, = fr к+1. Затем вычисляются следующие приближения для
тангенциальных составляющих магнитного поля (40) f к+1, уф, к+1 на границе расчетной области.
Шаг 3: Найденные fr, к+1, f к+1, Уф, к+1 используются в качестве граничных условий Bk+1 |r = fk+1 для нахождения следующего приближения B^, которое определяется из решения уравнения индукции (23).
Шаги 1 - 3 повторяются до тех пор, пока не будет достигнута сходимость на данном временном слое. В качестве критериев сходимости решения задачи (1) -(5) использовались следующие условия:
max|DIV Bk+j| <Sj; (41)
G
|f к+1 -f к||2 <^2. (42)
При достижении сходимости B(m) = Bk+1.
Для решения систем линейных алгебраических уравнений, полученных в результате дискретизации уравнений начально-краевой задачи (1) - (5), использована открытая программная библиотека итерационных решателей LIS (http://www.ssisc.org/lis/). Дискретные аналоги уравнений индукции магнитного поля и его потенциала решались стабилизированным метод бисопряжённых градиентов с предобуславливанием по методу Якоби.
4. Тестирование алгоритма
В качестве теста рассматривалась задача о размагничивании проводящего шара единичного радиуса, которая является частным случаем задачи (1) - (5) при
u = 0, vm = 1, r = 1: r < 1, t> 0:
dB
— + rotrot B = 0; (43)
dt
r > 1, t > 0:
V2y = 0; (44)
r = 1, t > 0:
B = Vy, (45)
r < 1, t = 0:
(rotB)r = 0 ; B = b, где div b = 0
(46)
Если начальное магнитное поле b = br(r, 0)er+ b0(r, 0)e0 (r < 1) имеет вид
cos 8 (sin nr Л , sin 8 (sin nr . Л
b =--I--cos nr I, b8=--I--cos nr-nr sin nr I, (48)
(nr)2 \ nr ) 2(nr)2 \ nr )
то поле
B(r, 8, t) = b(r, 8) (49)
и потенциал (r < 1)
2
„ ^ ч t cos 8
у (r, 8, t) =--(50)
2n r
являются аксиально-симметричным решением задачи (43) - (47) [33].
Расчеты проводились до момента времени T = 1 на сетках (nr х n0), где nr, n0 -количество контрольных объемов вдоль соответствующей координатной линии. На каждом временном слое критерии (41), (42) выполнялись с точностью е1 = е2 = 10-9. Точность численного решения для компонент индукции магнитного поля оценивалась по максимальным абсолютным AF и средним относительным 5F погрешностям:
ДF = max |F-F\, 5F = — f F-Fa^ dS-100%, (51)
°®(0,T ]' ^ V]a Fan
где F - сеточные значения компонент Br, B0; Fan - точные значения Br, B0, определяемые по (49); Fan = — f |Fan | d S, V - объем шара G.
VG
Точность численного решения для потенциала магнитного поля оценивалась на границе G по формулам
Ду = max Iw-wJ, 8у = — f l^-lsd. dS-100%, (52)
9G®(0,T] S dG Van
где у - сеточные значения потенциала; - точные значения потенциала (50); У an =1 f |van I dS ; S - площадь сферы 5G.
S SG
Результаты тестов на численную сходимость на последовательности вложенных расчетных сеток, приведенные в табл. 1, демонстрируют, что предложенный алгоритм решения начально-краевых задач (1) - (5) с вакуумными граничными условиями занимает промежуточное место по порядку аппроксимации между схемами первого и второго порядков по пространственным переменным. На рис. 1 - 3 представлены максимальные абсолютные погрешности для компонент индукции магнитного поля и потенциала, полученные на вложенных расчетных сетках.
Максимальные абсолютные и средние относительные погрешности
Сетка nr X n0 ДВг ДВ0 Ду SBrlt=0.25 ' % SB811 =0.25 ' % 1 =0.25 ' %
12 x 32 2.12-10-3 3.14-10-3 7.00-10-5 0.92 0.79 1.02
22 x 62 7.78-10-4 1.04-10-3 2.43-10-5 0.32 0.29 0.35
42 x 122 2.51-10-4 3.19-10-4 1.29-10-5 0.17 0.16 0.18
Рис. 1. Максимальная абсолютная погрешность для радиальной компоненты индукции магнитного поля: кр. 1 - сетка 12 х 32; кр. 2 - сетка 22 х 62; кр. 3 - сетка 42 х 122 Fig. 1. Maximum absolute error for the radial component of the magnetic field induction: 1, grid 12 х 32; 2, grid 22 х 62; 3, grid 42 х 122
Рис. 2. Максимальная абсолютная погрешность для меридиональной компоненты индукции
магнитного поля: кр. 1 - сетка 12 х 32; кр. 2 - сетка 22 х 62; кр. 3 - сетка 42 х 122 Fig. 2. Maximum absolute error for the meridional component of the magnetic field induction: 1, grid 12 х 32; 2, grid 22 х 62; 3, grid 42 х 122
Рис. 3. Максимальная абсолютная погрешность для потенциала на сфере (r = 1): кр. 1 - сетка 12 х 32; кр. 2 - сетка 22 х 62; кр. 3 - сетка 42 х 122 Fig. 3. Maximum absolute error for the potential on the sphere (r = 1): 1, grid 12 х 32; 2, grid 22 х 62; 3, grid 42 х 122
Заключение
Получен дискретный аналог уравнения индукции магнитного поля в произвольных ортогональных криволинейных координатах, решение которого удовлетворяет условию соленоидальности. Разработан новый алгоритм численного решения начально-краевой задачи с вакуумными граничными условиями для уравнения индукции магнитного поля в шаре. Сравнение результатов тестовых расчетов с аналитическим решением демонстрирует сходимость построенной разностной схемы.
ЛИТЕРАТУРА
1. Olson P.L., Glatzmaier G.A., Coe R.S. Complex polarity reversals in a geodynamo model // Earth Planet. Sci. Lett. 2011. V. 304. P. 168-179. DOI: 10.1016/j.epsl.2011.01.031.
2. Olson P., Driscoll P., Amit H. Dipole collapse and reversal precursors in a numerical dynamo // Phys. Earth. Planet. Inter. 2009. V. 173. P. 121-140. DOI: 10.1016/j.pepi.2008.11.010.
3. Jouve L., Brun A.S., Arlt R., Brandenburg A., Dikpati M., Bonanno A., Kapyla P.J., Moss D., RempelM., Gilman P., KorpiM.J., Kosovichev A.G. A solar mean field dynamo benchmark // Astron. Astrophys. 2008. V. 483. P. 949-960. DOI: 10.1051/0004-6361:20078351.
4. Glatzmaier G.A. Computer simulations of Jupiter's deep internal dynamics help interpret what Juno sees // Proc. Nat. Acad. Sci. 2018. V. 115 Iss. 27. P. 6896-6904. DOI: 10.1073/pnas. 1709125115.
5. Jones C. A dynamo model of Jupiter's magnetic field // Icarus. 2014. V. 241. P. 148-159. DOI: 10.1016/j.icarus.2014.06.020.
6. Jones C.A. Convection-driven geodynamo models // Phil. Trans. R. Soc. Lond. A. 2000. V. 358. P. 873-897. DOI: 10.1098/rsta.2000.0565.
7. Jackson A.A. et al. A spherical shell numerical dynamo benchmark with pseudo-vacuum magnetic boundary conditions // Geophysical Journal International. 2014. V. 196 Iss. 2. P. 712-723. DOI: 10.1093/gji/ggt425.
8. Bandaru V., Boeck T., Krasnov D., Schumacher J. A hybrid finite difference/boundary element procedure for the simulation of turbulent MHD duct flow at finite magnetic Reynolds number // Journal of Computational Physics. 2016. V. 304. P. 320-339. DOI: 10.1016/j.jcp. 2015.10.007.
9. Glatzmaier G.A. Numerical simulations of stellar convective dynamos I. The model and method // Journal of Computational Physics. 1984. V. 55. P. 461-484. DOI: 10.1016/0021-9991(84)90033-0.
10. ReshetnyakM.Yu. Simulation in Geodynamo. Saarbrucken: Lambert, 2013. 180 p.
11. Kageyama A., Miyagoshi T., Sato T. Formation of current coils in geodynamo simulations // Nature. 2008. V. 454. P. 1106-1109. DOI: 10.1038/nature07227.
12. Галанин М.П., Попов Ю.П. Квазистационарные электромагнитные поля в неоднородных средах: Математическое моделирование. М.: Наука; Физматлит, 1995. 320 с.
13. Hejda P., ReshetnyakM. Control volume method for the dynamo problem in the sphere with the free rotating inner core // Studia geoph. et geod. 2003. V. 47. P. 147-159. DOI: 10.1023/ A:1022207823737.
14. Куликовский А.Г., Погорелов Н.В, Семенов А.Ю. Математические вопросы численного решения гиперболических систем уравнений. М: Физматлит, 2001. 608 с.
15. Бетелин В.Б., Галкин В.А., Гореликов А.В. Алгоритм типа предиктор - корректор для численного решения уравнения индукции в задачах магнитной гидродинамики вязкой несжимаемой жидкости // Доклады Академии наук. 2015. T. 464. № 5. C. 525-528.
16. Колмычков В.В., Мажорова О.С., Федосеев Е.Э. Численный метод решения уравнений МГД // Препринты ИПМ им. М.В.Келдыша. 2009. № 30. 28 с.
17. Toth G. The V- B = 0 constraint in shock-capturing magnetohydrodynamics codes // Journal of Computational Physics. 2000. V. 161. P. 605-652. DOI: 10.1006/jcph.2000.6519.
18. Попов М.В., Елизарова Т.Г. Моделирование трехмерных МГД-течений в рамках магнитной квазигазодинамики // Препринты ИПМ им. М.В.Келдыша. 2013. № 23. 32 с.
19. Powell K.G., Roe P.L., Linde T.J., Gombosi T.I., De Zeeuw D.L. A Solution-Adaptive Upwind Scheme for Ideal Magnetohydrodynamics // Journal of Computational Physics. 1999. V. 154. p. 284-309. DOI: 10.1006/jcph.1999.6299.
20. Teyssier R., Commercon B. Numerical Methods for Simulating Star Formation // Frontiers in Astronomy and Space Sciences. 2019. V. 6. Article 51. DOI: 10.3389/fspas.2019.00051.
21. Ландау Л.Д., Лифшиц Е.М. Теоретическая физика: учеб. пособ. для вузов: в 10 т. Т. VIII. Электродинамика сплошных сред. М.: Физматлит, 2005. 656 с.
22. Куликовский А.Г., Любимов Г.А. Магнитная гидродинамика: учеб. пособие. 3-е изд. М.: Логос, 2011. 324 с.
23. Cowling T.G. Magnetohydrodynamics. Bristol: Adam Hilger, 1976. 136 p.
24. Patankar S. Numerical Heat Transfer and Fluid Flow. Washington DC: Hemisphere Publishing, 1980. 197 p.
25. Versteeg H.K., Malalasekera W. An Introduction to Computational Fluid Dynamics. Harlow: Pearson Education Limited, 2007. 503 p.
26. Ferziger J.H., Peric M. Computational Methods for Fluid Dynamics. Berlin: Springer, 2002. 423 p.
27. Taflove A., Hagness S. Computational Electrodynamics: The Finite-Difference Time-Domain Method, 2nd edition. Boston: Artech House, 2000. 866 p.
28. ГригорьевА.Д. Методы вычислительной электродинамики М.: Физматлит, 2012. 432 c.
29. Ates A., Altun O., Kilicman A. On A Comparison of Numerical Solution Methods for General Transport Equation on Cylindrical Coordinates // Applied Mathematics & Information Sciences. 2017. V. 11. No. 2. P. 433-439. DOI: 10.18576/amis/110211.
30. Dewangan S.K., Sinha S.L. Comparison of various numerical differencing schemes in predicting non-newtonian transition flow through an eccentric annulus with inner cylinder in
rotation // International Journal of Mechanical and Industrial Engineering. 2013. V. 3. Iss. 1. P. 119-129.
31. Arshad K., Rafique M., Majid A., Jabeen S. Analyses of MHD Pressure Drop in a Curved Bend for Different Liquid Metals // Journal of Applied Sciences. 2007. V. 7. Iss. 1. P. 72-78. DOI: 10.3923/jas.2007.72.78.
32. Тихонов А.Н., Самарский А.А. Уравнения математической физики: учебник для ун-тов. 4-е изд., испр. М.: Наука, 1972. 746 с.
33. Гореликов А.В. Аксиально-симметричное решение задачи о размагничивании проводящего шара // Вестник кибернетики. 2018. № 4(32). C. 9-15.
Статья поступила 07.09.2019 г.
Bychin I.V., Gorelikov A.V., Ryakhovskii A.V. (2020) NUMERICAL SOLUTION OF THE INITIAL BOUNDARY VALUE PROBLEM WITH VACUUM BOUNDARY CONDITIONS FOR THE MAGNETIC FIELD INDUCTION EQUATION IN A BALL. Vestnik Tomskogo gosudarstvennogo universiteta. Matematika i mekhanika [Tomsk State University Journal of Mathematics and Mechanics]. 64. pp. 15-30
DOI 10.17223/19988621/64/2
Keywords: induction equation, vacuum boundary conditions, numerical solution.
The aim of this work is to develop and test an algorithm for numerical solution of the initial-boundary value problem with vacuum boundary conditions for the equation of magnetic field induction in a ball based on a modification of finite-difference methods of computational fluid dynamics and electrodynamics for use in orthogonal curvilinear coordinates.
The problem of the numerical solution of the magnetic field induction equation in the ball G is considered in this paper using the magnetohydrodynamics (MHD) model. It is assumed that outside the ball the magnetic field induction B = grad y, and the potential у is a harmonic function regular at infinity. At the boundary dG , the vacuum boundary conditions are set. They include the requirement that the magnetic field is continuous and the normal component of the current density is equal to zero.
The discretization of the induction equation is carried out using the control-volume method and the FDTD (Finite Difference Time Domain) method modified for orthogonal curvilinear coordinates. When integrating over time, a completely implicit scheme is used. The total (convective and diffusion) flows are approximated according to the power-law scheme.
The proposed algorithm consists in sequentially executing the following steps at each time layer:
Step 1: The approximation of the radial component of the magnetic field at dG is determined from the discretized continuity equation for the boundary control volumes.
Step 2: The external Neumann problem for the Laplace equation is solved using the Kelvin transformation. The approximation for the tangential components of the magnetic field at the dG is calculated using the found potential y.
Step 3: The found components of B are used as boundary conditions for finding a solution to the induction equation.
The steps 1 to 3 are repeated until convergence is achieved at a given moment in time. The algorithm was tested on the problem of the magnetic field diffusion in a conducting ball, which has an analytical solution.
In this paper, the discretization of the magnetic field induction equation the solution of which satisfies the solenoidality condition is obtained in arbitrary orthogonal curvilinear coordinates. A new algorithm for numerical solution of the initial-boundary value problem with vacuum boundary conditions for the equation of magnetic field induction in a ball is developed. The comparison between the numerical solution and the analytical solution demonstrates the convergence of the constructed finite-difference scheme.
Financial support. The work was supported by the FSR Program of the state academies of Sciences for 2013-2020, project № 0065-2019-0021.
AMS Mathematical Subject Classification: 76W05
Igor V. BYCHIN(Surgut State University, Branch of Federal State Institution "Scientific Research Institute for System Analysis of the Russian Academy of Sciences" in the Surgut, Surgut, Russian Federation). E-mail: [email protected]
Andrey V. GORELIKOV(Candidate of Physics and Mathematics, Surgut State University, Branch of Federal State Institution "Scientific Research Institute for System Analysis of the Russian Academy of Sciences" in the Surgut, Surgut, Russian Federation) E-mail: [email protected]
Alexei V.RYAKHOVSKII (Candidate of Physics and Mathematics, Surgut State University, Surgut, Russian Federation). E-mail: [email protected]
REFERENCES
1. Olson P.L., Glatzmaier G.A., Coe R.S.. (2011) Complex polarity reversals in a geodynamo model. Earth Planet. Sci. Lett. 304. pp. 168-179. DOI: 10.1016/j.epsl.2011.01.031.
2. Olson P., Driscoll P., Amit H. (2009) Dipole collapse and reversal precursors in a numerical dynamo. Phys. Earth. Planet. Inter. 173. pp. 121-140. DOI: 10.1016/j.pepi.2008.11.010.
3. Jouve L., Brun A.S., Arlt R., Brandenburg A., Dikpati M., Bonanno A., Käpylä P.J., Moss D., Rempel M., Gilman P., Korpi M. J., Kosovichev A.G. (2008) A solar mean field dynamo benchmark. Astron. Astrophys. 483. pp. 949-960. DOI: 10.1051/0004-6361:20078351.
4. Glatzmaier G.A. (2018) Computer simulations of Jupiter's deep internal dynamics help interpret what Juno sees. Proc. Nat. Acad. Sci. 115 (27). pp. 6896-6904. DOI: 10.1073/pnas. 1709125115.
5. Jones C. (2014). A dynamo model of Jupiter's magnetic field. Icarus. 241. pp. 148-159. DOI: 10.1016/j.icarus.2014.06.020.
6. Jones C.A. (2000) Convection-driven geodynamo models. Phil. Trans. R. Soc. Lond. A. 358. pp. 873-897. DOI: 10.1098/rsta.2000.0565.
7. Jackson A.A. et al. (2014) A spherical shell numerical dynamo benchmark with pseudo-vacuum magnetic boundary conditions. Geophysical Journal International. 196 (2). pp. 712723. DOI: 10.1093/gji/ggt425.
8. Bandaru V., Boeck T., Krasnov D., Schumacher J. (2016) A hybrid finite difference/boundary element procedure for the simulation of turbulent MHD duct flow at finite magnetic Reynolds number. Journal of Computational Physics. 304. pp. 320-339. DOI: 10.1016/j.jcp.2015. 10.007.
9. Glatzmaier G.A. (1984) Numerical simulations of stellar convective dynamos I. The model and method. Journal of Computational Physics. 55. pp. 461-484. DOI: 10.1016/0021-9991(84)90033-0.
10. Reshetnyak M.Yu. (2013) Simulation in Geodynamo. Saarbrucken: Lambert.
11. Kageyama A., Miyagoshi T., Sato T. (2008) Formation of current coils in geodynamo simulations. Nature. 454. pp. 1106-1109. DOI: 10.1038/nature07227.
12. Galanin M.P., Popov Yu.P. (1995) Kvazistacionarnye elektromagnitnye polya v neodnorodnykh sredakh: Matematicheskoe modelirovanie. [Quasistatonary electromagnetic fields in inhomogeneous media: Mathematical simulation]. Moscow: Nauka; Fizmatlit.
13. Hejda P., Reshetnyak M. (2003) Control volume method for the dynamo problem in the sphere with the free rotating inner core. Studia Geoph. et Geod. 47. pp. 147-159. DOI: 10.1023/A:1022207823737.
14. Kulikovskij A.G., Pogorelov N.V, Semenov A.Yu. (2001). Matematicheskie voprosy chislennogo resheniya giperbolicheskikh sistem uravneniy. [Mathematical problems of numerical solution of hyperbolic systems of equations]. Moscow: Fizmatlit.
15. Betelin V.B., Galkin V.A., Gorelikov A.V. (2015). Algoritm tipa prediktor-korrektor dlya chislennogo resheniya uravneniya indukcii v zadachah magnitnoj gidrodinamiki vyazkoj
neszhimaemoj zhidkosti [Predictor-corrector type algorithm for the numerical solution of the induction equation in problem of magnetic hydrodynamics of viscous incompressible fluid]. Doklady Akademii nauk. 464(5). pp. 525-528. DOI: 10.7868/S0869565215290034.
16. Kolmychkov V.V., Mazhorova O.S., Fedoseev E.E. (2009) Chislennyy metod resheniya uravneniy MGD [Numerical Method for MHD Equations]. Preprinty IPMim. M.V. Keldysha. 30. 28 pages.
17. Toth G. (2000) The V- B = 0 constraint in shock-capturing magnetohydrodynamics codes. Journal of Computational Physics. 161. pp. 605-652. DOI: 10.1006/jcph.2000.6519.
18. Popov M.V., Elizarova T.G. (2013) Modelirovanie trekhmernykh MGD-techeniy v ramkakh magnitnoy kvazigazodinamiki [Simulation of 3D MHD flows within the framework of magneto-quasi-gasdynamics]. Preprinty IPM im. M.V. Keldysha. 23. 32 pages.
19. Powell K.G., Roe P.L., Linde T.J., Gombosi T.I., De Zeeuw D.L. (1999) A solution-adaptive upwind scheme for ideal magnetohydrodynamics. Journal of Computational Physics, 154. pp. 284-309. DOI: 10.1006/jcph.1999.6299.
20. Teyssier R., Commercon B. (2019) Numerical methods for simulating star formation. Frontiers in Astronomy and Space Sciences. 6(51). DOI: 10.3389/fspas.2019.00051.
21. Landau L.D., Lifshitz E.M. (1984) Electrodynamics of Continuous Media (Volume 8 of A Course of Theoretical Physics). Oxford: Pergamon Press.
22. Kulikovskiy A.G., Lyubimov G.A. (2011) Magnitnaya gidrodinamika [Magnetohydrodyna-mics]. Moscow: Logos.
23. Cowling T.G. (1976)Magnetohydrodynamics. Bristol: Adam Hilger.
24. Patankar S. (1980) Numerical Heat Transfer and Fluid Flow. Washington DC: Hemisphere Publishing.
25. Versteeg H.K., Malalasekera W. (2007) An Introduction to Computational Fluid Dynamics. Harlow: Pearson Education Limited.
26. Ferziger J.H., Peric M. (2002) Computational Methods for Fluid Dynamics. Berlin: Springer.
27. Taflove A., Hagness S. (2000) Computational Electrodynamics: The Finite-Difference TimeDomain Method. 2nd edition. Boston: Artech House.
28. Grigorev A.D. (2012) Metody vychislitelnoy elektrodinamiki [Methods of computational electrodynamics]. Moscow: Fizmatlit.
29. Ates A., Altun O., Kilicman A. (2017) On a comparison of numerical solution methods for general transport equation on cylindrical coordinates. Applied Mathematics & Information Sciences. 11 (2). pp. 433-439. DOI: 10.18576/amis/110211.
30. Dewangan S.K., Sinha S.L. (2013) Comparison of various numerical differencing schemes in predicting non-Newtonian transition flow through an eccentric annulus with inner cylinder in rotation. International Journal of Mechanical and Industrial Engineering. 3 (1). pp. 119-129.
31. Arshad K., Rafique M., Majid A., Jabeen S. (2007) Analyses of MHD pressure drop in a curved bend for different liquid metals. Journal of Applied Sciences. 7 (1). pp. 72-78. DOI: 10.3923/jas.2007.72.78.
32. Tihonov A.N., Samarskij A.A. (1972) Uravneniya matematicheskoy fiziki [Equations of mathematical physics]. Moscow: Nauka.
33. Gorelikov A.V. (2018) Aksialno-simmetrichnoe reshenie zadachi o razmagnichivanii provodyashchego shara [An axisymmetric solution of the problem about demagnetization of a conducting ball]. Vestnikkibernetiki. 4(32). pp. 9-15.
Received: September 7, 2019
