Numerical solution of the initial boundary value problem with vacuum boundary conditions for the magnetic field inductio.pdf Необходимость решения уравнения магнитной индукции с вакуумными граничными условиями зачастую возникает при численном решении задач магнитной гидродинамики, например при математическом моделировании явления гид-ромагнитного динамо [1-5]. При постановке задач с вакуумными граничными условиями магнитное поле считается всюду непрерывным, а в вакууме - потенциальным и соленоидальным [6]. В коммерческих программных кодах обычно используются упрощенные так называемые «псевдо-вакуумные» граничные условия, при которых на границе области полагается равной нулю тангенциальная составляющая магнитного поля [7, 8]. Решения, полученные с использованием псевдо-вакуумных условий, обычно корректно воспроизводят плотность тока, однако структура магнитного поля значительно отличается от решений, рассчитанных при вакуумных условиях [8]. Поэтому в задачах моделирования гидромаг-нитного динамо представляется целесообразным использование более реалистичных вакуумных граничных условий. При численном решении уравнения индукции, возникает проблема получения соленоидального магнитного поля. В настоящее время наиболее распространенными методами вычислительной магнитной гидродинамики являются псевдо-спектральные [9, 10] и конечно-разностные [11, 12] методы, а также метод контрольного объема [10, 13]. В рамках данных методов существуют различные подходы для обеспечения бездивергентности магнитного поля [14]: использование векторного потенциала; метод искусственного скалярного потенциала [15]; представление дискретного аналога уравнения индукции в такой форме, чтобы его решение автоматически удовлетворяло сеточному уравнению неразрывности [16-18]; метод Пауэлла, который широко применяется в астрофизических приложениях и основан на записи уравнений в форме, допускающей существование магнитных монополей [19]. Перечисленные методы обладают своими преимуще- 1 Работа выполнена при поддержке Программы ФНИ государственных академий наук на 2013-2020 гг., проект № 0065-2019-0021. 16 И.В. Бычин, А.В. Гореликов, А.В. Ряховский ствами и недостатками. Например, коррекция магнитного поля в методе скалярного потенциала нелокальна, и поэтому локальные ошибки в дивергенции магнитного поля могут мгновенно распространяться на всю область [20]. Другие методы, например метод Пауэлла, не являются, строго говоря, консервативными. Несмотря на достаточно большое число работ в этой области, разработка эффективных алгоритмов решения уравнения индукции с вакуумными граничными условиями остается актуальной задачей вычислительной магнитной гидродинамики. Целью настоящей работы является разработка и апробация алгоритма численного решения начально-краевой задачи с вакуумными граничными условиями для уравнения индукции магнитного поля в шаре на основе модификации разностных методов вычислительной гидродинамики и электродинамики для использования в ортогональных криволинейных координатах. 1. Постановка задачи и математическая модель Пусть проводящая жидкость заполняет шар G с центром в начале координат и радиусом a, который находится в неподвижной диэлектрической среде. Предполагается, что удельная электрическая проводимость жидкости σ = const. В рамках модели МГД [21-23] рассматривается задача о нахождении индукции магнитного поля B, при заданной скорости жидкости u. Задача решается в сферических координатах r, θ, φ. Полагается, что вне шара магнитное поле потенциально, и его скалярный потенциал ψ является регулярной на бесконечности гармонической функцией. На границе проводника и диэлектрика задаются условия сопряжения (так называемые вакуумные граничные условия), которые заключаются в требовании непрерывности магнитного поля и равенства нулю нормальной составляющей c плотности тока j = - rot B [21]. Таким образом, начально-краевая задача, для ко-4π торой в данной работе строится процедура численного решения, имеет вид r < a, t > 0: г > a, t ≥ 0: г = a, t ≥ 0: г ≤ a, t = 0: ∂B ∂t V2ψ= -rot(u×B)+νmrotrot B=0; (1) 0 (ψ - регулярна при r → +∞); (2) B =Vψ, (3) (rot B)r = 0; (4) B= b, где div b = 0. (5) Здесь t - время, c скорость света, νm = ---коэффициент магнитной вязкости. 4πσ Из (1) следует, что равенство нулю дивергенции начального поля b обеспечивает соленоидальность магнитного поля B и в последующие моменты времени. Численное решение начально-краевой задачи с вакуумными граничными условиями 17 2. Численный метод решения уравнения индукции В данной работе для дискретизации уравнения индукции использован метод контрольного объема [24-26] и модификация метода FDTD (Finite Difference Time Domain) [17, 27, 28] в ортогональных криволинейных координатах, учитывающая специфику рассматриваемой задачи (1) - (5). Дискретизация расчетной области Пусть {Xα} - произвольные ортогональные координаты, {eα} - соответствующий правый ортонормированный базис, {Hα} - коэффициенты Ламе (здесь и далее, нижние греческие индексы принимают значения 1, 2, 3 и определяют номер координаты или соответствующей компоненты). Предполагается, что расчетная область G - ограниченная область в пространстве, которую можно разбить системой координатных поверхностей {Xα,iα =const, α=1,2,3, iα =1,...,nα : Xα,iα-1 < Xα,iα} на непересекающиеся подобласти - контрольные объемы ^∕1∕2∕3 = {(x1, x2, X3): Xα∈ (Xa,≈α-1, xO,,,,, ), α = 1, 2, 3}, iα = 2,^,nα. В центре каждого контрольного объема содержится внутренняя расчетная точка P1ij^3 . В целях сокращения записи, в дальнейшем набор индексов i1, i2, i3 обозначается как (i). Сеточные значения компонент магнитного поля рассчитываются в точках на гранях контрольных объемов, т.е. используется дискретизация на смещённых расчетных сетках [24-28]. Для описания процедуры дискретизации на смещённых сетках удобно использовать операторы сдвига на наборах индексов hασ(i)≡hασ(i1,i2,i3)=hασ(i1),hασ(i2),hασ(i3),(6) действие которых определяется по правилу: hασ (iβ) = iβ+∣ δαβ, O = 1,2,3, σ=-1,1, (7) где δαβ - символ Кронекера, а индекс σ определяет направление и величину сдвига по индексу iβ. При этом полагается, что расчетные точки с одним полуцелым индексом Phσ /2(i) находятся в центрах граней Shσ /2(i) контрольного объема D(i); hα (i) hα (i) δShσ /2(i) - площадь соответствующей грани контрольного объема; граница кон- 3 трольного объема ∂D(i) =∪ Sh±1/2(i). Точки с двумя полуцелыми индексами hα (i) α=1 центрах рёбер lrh^'2hα/2(≈) = shχ'2(≈) ∩ shσ/2i) (Y ≠ α ≠ β) D(i) расположенных вдоль координатных линий Xγ; соответствующего ребра контрольного объема D(i), и находятся в phχ/^hl '2(≈) контрольного объема δlγ^ hχX/2^σσ^/2(i) - длина dSha.’^(i) lγ, h±^^^hα.’^(i) ∪ lβ, hγ1^hα.’^(i)' 18 И.В. Бычин, А.В. Гореликов, А.В. Ряховский Дискретный аналог уравнения индукции Уравнение (1) является следствием закона электромагнитной индукции Фарадея: ∂B (8) - = - c rot E, ∂t где E - напряженность электрического поля, которая в модели МГД имеет вид [22] 1c (9) E = - (и × B) +--rot B . c 4 πσ Для получения дискретного аналога уравнение (8) умножается на eαdSαdt (dSα -элемент площади на грани Shσ /2(.)), и затем интегрируется по поверхности hα (.) Shσ/2,.) и по промежутку времени [Ґ0, ^]. При вычислении интеграла по поверхно-hα (.) сти S σ/2 .) используется теорема Стокса, а при интегрировании по времени -hα (.) полностью неявная схема. В результате, дискретный аналог уравнения (8) для компоненты Bα можно записать в виде 13 ^δtt (Bα,hα '^(^i)- Bα,hσ '^(i)) δshσ/^а) = -cβ∑ι εαβγ∆β (δ^γ ^γ^≈α ^2,.) (10) Здесь B hσ/2(.) ≡ Ba(P σ/2(.),t) - сеточн^іе значения Ba на текущем временном слое; α,hα(.) hα(.) B0 /2 ) ≡ Ba(P /2 t00) - сеточные значения Ba на предыдущем временном слое; α,hα(.) hα(.) £■ hχ/2 σ/2..) ≡ Eγ(P /2 /2. t) - сеточные значения компоненты напряженности γ,hβha(.) hβha(.) электрического поля в точках на ребрах контрольного объема (a ≠ β ≠ γ); δt = t - t0 - шаг по времени; εaβγ - символ Леви-Чивиты; конечные разности в правой части (10) определяются по формуле δ- (δlγeY 1hασ/2(i) = ey, h1/2hσ/2(i) δlγ, hβ 2h;;/^(^i) - e,, h-12hζ,'2(i) ^lγ, h-12hσ/2(i). (11) Анализ уравнения (10) для задачи (1) - (5) показывает, что сеточные значения (rot B)r в точках на границе расчетной области не входят в уравнение (10) для приграничных расчетных точек Pn. . . Следовательно, условие (4) равенства нулю радиальной составляющей ротора магнитного поля на границе шара G не накладывает каких-либо дополнительных ограничений на вид дискретного аналога уравнения индукции для приграничных расчетных точек. Дискретный аналог уравнения неразрывности Сеточные значения div B определяются во внутренних расчетных точках P(i) контрольных объемов D(i) по формуле DIV B L ∑(β ,1/2,. ,δS./, - B ,-1/,. ,δ∖-ι/, . 1, (12) ^('■) α,h«- (’) h«- (i) α,h«- (i) h«- (i)' (i ) α=1 где δV(i) - объем D(i^ и учтено направление внешней нормали n,σ∕2,.) =σeα на hα (i) гранях контрольного объема D(.). Пусть D(.) - внутренний контрольный объем (т.е. его граница не пересекается с границей расчетной области G) и пусть в момент Численное решение начально-краевой задачи с вакуумными граничными условиями 19 времени t0 сеточные значения компонент магнитного поля Bα0,hσ /2(i) удовлетворя ют дискретному аналогу уравнения неразрывности DIV Бч = 0 . Тогда подста- P(i) новка B σ/2(.)δS /2 ,) из сеточного уравнения индукции (10) в правую часть (12) α,hα(i)hα(i) дает d'vB =-δKt- ε“β’(∆β a: V2ψ=0; (35) r = a: dΨ= B,^ ∂r r (36) С помощью преобразования обратных радиусов [32]: rr ' = a^ , ψ(r, θ, φ) = r' ω(r' , θ, φ), (37) внешнюю задачу Неймана (35), (36) можно свести к решению внутренней краевой задачи для уравнения Лапласа на ω с граничными условиями третьего рода: r’< a : V2ω = 0; (38) ∂ω (39) r = a : a--+ ω = -Br. ∂r' Краевая задача (38) - (39) численно решалась методом контрольного объема в сферических координатах на той же расчетной сетке, что и уравнение индукции (23). Сеточные значения ω вычислялись во внутренних расчетных точках P(i) контрольных объемов D(i) и в точках Ph1/2(i) ∈∂G . После этого, по (37) определялся h1 (i) потенциал ψ в точках P1/2 ∈∂G, а затем тангенциальные компоненты магнит-h1 (i) ного поля r= a: B = 1 ∂ψ , B = 1 ∂ψ (40) r = a : Bθ = , Bφ =-- θ r ∂θ φ rsin θ∂φ в точках phχ,/2hι∕2(i) ∈∂G (β ≠ 1). Поскольку тангенциальные компоненты магнитного поля на границе шара (40) находятся, как компоненты Vψ , то можно показать, что сеточные значения нормальной составляющей ротора магнитного поля (rot B)r в точках Pσ/21.х/211/2/-^ ∈∂G (γ ≠ β ≠ 1) автоматически становятся равными нулю, тем са-hγ hβ h1 (i) мым обеспечивая выполнение условия (4). Алгоритм решения начально- краевой задачи (1) - (5) Для удобства записи алгоритма сеточные значения B = Brer+Bθeθ + Bφeφ в точках на границе шара G обозначаются f (т.е. f ≡ В|r=^). Пусть в момент времени tm^1 с заданной точностью найдено численное решение B(m-1) начально-краевой задачи (1) - (5) и пусть нижний индекс k нумерует внутренние итерации, которые необходимо сделать, чтобы получить решение на следующем временном слое tm. Численное решение начально-краевой задачи с вакуумными граничными условиями 23 В качестве нулевого приближения для B на временном слое tm используется B(m - 1), т.е. при k = 0 B0 = B(m-1). В данной работе предложен и реализован следующий алгоритм получения численного решения начально-краевой задачи (1) - (5): Шаг 1: Из сеточного уравнения неразрывности (15) для приграничных контрольных объемов Dn1i2i3 , определяется нормальная (радиальная) компонента магнитного поля на границе G, т.е. вычисляется fr, k+1. Шаг 2: Решается третья внутренняя краевая задача (38), (39) на потенциал ω, для которой Br∣r,=α = fr k+1. Затем вычисляются следующие приближения для тангенциальных составляющих магнитного поля (40) fθ, k+1, fφ, k+1 на границе расчетной области. Шаг 3: Найденные fΓ, k+ι, fθ, k+ι, fφ, k+ι используются в качестве граничных условий Bk+1 Ir =a = fk+1 для нахождения следующего приближения Bk+1, которое определяется из решения уравнения индукции (23). Шаги 1 - 3 повторяются до тех пор, пока не будет достигнута сходимость на данном временном слое. В качестве критериев сходимости решения задачи (1) -(5) использовались следующие условия: maxIDiv Bk+J ≤ε1; (41) G (42) Ilf k+1- f k∣∣2 ≤ε2. При достижении сходимости B(m) = Bk+1. Для решения систем линейных алгебраических уравнений, полученных в результате дискретизации уравнений начально-краевой задачи (1) - (5), использована открытая программная библиотека итерационных решателей LIS (http://www.ssisc.org/lis/). Дискретные аналоги уравнений индукции магнитного поля и его потенциала решались стабилизированным метод бисопряжённых градиентов с предобуславливанием по методу Якоби. 4. Тестирование алгоритма В качестве теста рассматривалась задача о размагничивании проводящего шара единичного радиуса, которая является частным случаем задачи (1) - (5) при u ≡ 0, ν m = 1, r = 1: r < 1, t > 0: ∂B --+ rot rot B = 0 ; ∂t (43) r > 1, t ≥ 0: V2ψ=0; (44) r = 1, t ≥ 0: B =Vψ, (45) (rotB)r = 0; (46) r ≤ 1, t = 0: B= b, где div b = 0. (47) 24 И.В. Бычин, А.В. Гореликов, А.В. Ряховский Если начальное магнитное поле b = br(r, θ)er+ bθ(r, θ)eθ (r ≤ 1) имеет вид , cos θ (sin πr '∖ , sin θ Ґ sin πr . br =------1--cos πr I, bθ=-------1--cos πr-πr sin πr I, 22 (48) то поле (49) и потенциал (r ≤ 1) e-π t cos θ ψ(r, θ, t) = - 2 2 2 2π I (50) являются аксиально-симметричным решением задачи (43) - (47) [33]. Расчеты проводились до момента времени T = 1 на сетках (nI × nθ), где nI, nθ -количество контрольных объемов вдоль соответствующей координатной линии. На каждом временном слое критерии (41), (42) выполнялись с точностью ε1 = ε2 = 10-9. Точность численного решения для компонент индукции магнитного поля оценивалась по максимальным абсолютным ∆F и средним относительным δF погрешностям: ∆F = max F -FaI, δF = - ∫ F _Fa^ d3∙100% , g^(0,TV G Fan (51) где F - сеточные значения компонент BI, Bθ; Fan - точные значения BI, Bθ, опре деляемые по (49); Fan = - ∫ ∣Fan ∖d З, V - объем шара G. VG Точность численного решения для потенциала магнитного поля оценивалась на границе G по формулам ∆ψ= max ∣ψ-ψaJ, δψ= - ∫ ∣ψ,ψa^ dS∙100%, ∂G^(fi,T] S ∙^g (52) ψan где ψ - сеточные значения потенциала; ψan - точные значения потенциала (50); ψan = -1 ∫ ∣ψan I dS; S - площадь сферы δG. S ∂ G Результаты тестов на численную сходимость на последовательности вложенных расчетных сеток, приведенные в табл. 1, демонстрируют, что предложенный алгоритм решения начально-краевых задач (1) - (5) с вакуумными граничными условиями занимает промежуточное место по порядку аппроксимации между схемами первого и второго порядков по пространственным переменным. На рис. 1 - 3 представлены максимальные абсолютные погрешности для компонент индукции магнитного поля и потенциала, полученные на вложенных расчетных сетках. Максимальные абсолютные и средние относительные погрешности Сетка Пг × nθ ∆Br ∆Bθ ∆ψ δBrlt=0.25 , % δBθlt ==0^25, % δψlt=0.25 , % 12 × 32 2.12∙10-3 3.14∙10-3 7.00∙10-5 0.92 0.79 1.02 22 × 62 7.78∙10-4 1.04∙10-3 2.43∙10-5 О.32 0.29 О.35 42 × 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 Заключение Получен дискретный аналог уравнения индукции магнитного поля в произвольных ортогональных криволинейных координатах, решение которого удовлетворяет условию соленоидальности. Разработан новый алгоритм численного решения начально-краевой задачи с вакуумными граничными условиями для уравнения индукции магнитного поля в шаре. Сравнение результатов тестовых расчетов с аналитическим решением демонстрирует сходимость построенной разностной схемы.

Скачать электронную версию публикации