научная статья по теме АЛГОРИТМ ТИПА ПРЕДИКТОР–КОРРЕКТОР ДЛЯ ЧИСЛЕННОГО РЕШЕНИЯ УРАВНЕНИЯ ИНДУКЦИИ В ЗАДАЧАХ МАГНИТНОЙ ГИДРОДИНАМИКИ ВЯЗКОЙ НЕСЖИМАЕМОЙ ЖИДКОСТИ Математика

Текст научной статьи на тему «АЛГОРИТМ ТИПА ПРЕДИКТОР–КОРРЕКТОР ДЛЯ ЧИСЛЕННОГО РЕШЕНИЯ УРАВНЕНИЯ ИНДУКЦИИ В ЗАДАЧАХ МАГНИТНОЙ ГИДРОДИНАМИКИ ВЯЗКОЙ НЕСЖИМАЕМОЙ ЖИДКОСТИ»

ДОКЛАДЫ АКАДЕМИИ НАУК, 2015, том 464, № 5, с. 525-528

^ МАТЕМАТИЧЕСКАЯ

ФИЗИКА

УДК 517.9+519.87+532+536

АЛГОРИТМ ТИПА ПРЕДИКТОР-КОРРЕКТОР ДЛЯ ЧИСЛЕННОГО РЕШЕНИЯ УРАВНЕНИЯ ИНДУКЦИИ В ЗАДАЧАХ МАГНИТНОЙ ГИДРОДИНАМИКИ ВЯЗКОЙ НЕСЖИМАЕМОЙ ЖИДКОСТИ © 2015 г. Академик РАН В. Б. Бетелин, В. А. Галкин, А. В. Гореликов

Поступило 19.05.2015 г.

Получено точное 3Б-решение системы уравнений МГД, удовлетворяющее условиям скольжения по границе параллелепипеда. Проведен анализ эффективности итерационной процедуры, основанной на последовательном расщеплении вычислений для магнитного поля и гидродинамического течения на основе детального сравнения приближений с точным решением задачи. Результаты тестовых расчетов на последовательности вложенных сеток соответствуют второму порядку аппроксимации вычислительного метода по пространственным координатам.

DOI: 10.7868/S0869565215290034

Исследуется численный алгоритм решения уравнения индукции в задачах магнитной гидродинамики, позволяющий получить магнитное поле, с заданной точностью удовлетворяющее условию соленоидальности. Система уравнений магнитной гидродинамики (МГД) вязкой несжимаемой жидкости в изотермическом случае в декартовых координатах {ха, а = 1,2,3} имеет вид

ди + (иУ)и = --Vр + v Ди + —(го1В х В), (1) дг р 4пр

Шу и = 0, (2)

— = го1;(и х В) + V тАВ, (3)

дг

Шу В = 0, (4)

где г — время; и — скорость жидкости; р — давление; В — индукция магнитного поля; р — плотность; V — кинематическая вязкость; vm — магнитная вязкость.

Дискретные аналоги уравнений (1)—(4) получены методом контрольного объема [1] с использованием неявной схемы и схемы со степенным законом для аппроксимации конвективных и диффузионных потоков на гранях контрольных объемов. При дискретизации использовались разнесенные расчетные сетки, т.е. компоненты скорости и индукции магнитного поля рассчиты-

Научно-исследовательский институт системных исследований Российской Академии наук, Москва Политехнический институт Сургутского государственного университета E-mail: val-gal@yandex.ru

вались на гранях основных контрольных объемов. Вычисления основаны на итерационной процедуре расщепления с последовательным определением полей скорости течения и давления, а затем магнитного поля. В частности, в настоящей работе использован алгоритм PISO (Pressure-Implicit with Splitting of Operators) [1, 2]. Используя тождество

rot(u X B) = u divB + (BV)u - (uV)B - B div u, получаем уравнение индукции в виде d в

dBa + div (uBa -vmgradBJ = div(B«J, а = 1,2,3. (5)

dt

При построении дискретного аналога уравнение (5) интегрируется по контрольному объему D с центром в точке P и по времени с использованием полностью неявной схемы. Граница контрольного объема д D + ориентирована по внешней нормали n и состоит из шести гладких поверхностей Sb (b = 1,2,...,6), ортогональных координатным осям. Индексы n и n + 1 обозначают сеточные значения функций в последовательные моменты времени tn и tn+1 (At = tn+1 - tn — шаг по времени). Дискретный аналог уравнения индукции (5) записывается в виде

f (С1 - впа) = H(Bn:1) + Qa, а = 1,2,3, (6) At

6

где H(ВГ1) = -1X ab(Bf(Pb) - Bf(P)) - конечно-

KDb=1

разностный оператор, аппроксимирующий конвективные и диффузионные потоки на гранях контрольного объема, Pb — соседние с P точки; коэффициенты ab рассчитываются по схеме со

526 БЕТЕЛИН и др.

(а) (б)

4

2 0 -2 -4

-6

Рис. 1. Структура течения и изолинии компоненты В2 (Ву) индукции магнитного поля в сечении х2 =п при t = 0.5 (а — численное решение, б — аналитическое решение).

степенным законом [1] (Ь = 1,2,..., 6); Qa — линеаризованные члены источника; УБ — объем области Б.

Дискретный аналог уравнения неразрывности получается путем интегрирования (4) по контрольному объему Б:

6

|ё1уВйУ = |В • пйБ * £В"ь+1АБь = Аа В"а+1,

Б

где Вь = В • п|

дБ

Ь=1

Бь

А а ва+1 = 0,

Дt

(В* - ВП) = Н (В*) + Оа, а = 1,2,3.

(8)

А а В** = 0.

А а Ф =А а В*,

(10)

где А а — разностный аналог оператора Лапласа.

Предполагается, что на границе расчетной области О задана нормальная составляющая вектора магнитной индукции В , интегрально удовлетворяющая условию соленоидальности:

Вп йБ = 0,

о**| _ о*

Вп \ЯО — Вп

— Вп, тогда

дО

дВ" 'дО

= 0,

(11)

дО

(7)

где А а — разностный аналог частной производной д

дха

На первом шаге алгоритма (предиктор) находится первое приближение магнитного поля В* из уравнения (6):

1

На втором шаге производится корректировка поля В* так, чтобы новое поле В** было соленои-дальным:

т.е. потенциал ф определяется как решение краевой задачи на уравнение Пуассона (10) с однородными граничными условиями второго рода (11). Шаги алгоритма повторяются до тех пор, пока не будет достигнута сходимость в уравнениях (6), (7) на данном временном слое. В частности, максимальное по модулю сеточное значение дивергенции В должно быть меньше некоторого наперед заданного е , т.е.

тах |ШуВ(Р)| <е,

V РеО

где

6

ё1у В(Р) «V- |В • пйБ «V-£ВьАБь.

дБ

(9)

Поле В** ищется в виде В** = В* + 5 В, причем поправочное поле 5 В полагается потенциальным: 5 В = ^гаё ф. Уравнение на потенциал ф получается из условия соленоидальности В**:

Ь=1

При достижении сходимости В,П+1 = В**.

Авторами получена серия точных решений системы уравнений МГД (1)—(4), которые использовались для тестирования разработанного численного алгоритма. Одно из них представлено ниже. В частности, при р = 1 г/см3, V = Vт = = 0.1см2/с точное решение системы уравнений МГД (1)—(4) имеет вид

ДОКЛАДЫ АКАДЕМИИ НАУК том 464 № 5

2015

АЛГОРИТМ ТИПА ПРЕДИКТОР-КОРРЕКТОР 527

Таблица 1. Максимальные абсолютные погрешности

«1 X n2 X n3 AUj аи2 Au3 abj aB2 aB3

30 x 15 x 30 0.11 0.097 0.13 0.25 0.40 1.17

60 x 30 x 60 0.033 0.027 0.037 0.076 0.14 0.35

120 x 60 x 120 0.0087 0.0069 0.0098 0.021 0.043 0.096

Таблица 2. Максимальные относительные погрешности (%)

n1 x n2 x n3 s«! s«2 s«3 sB! sb2 sB3

30 x 15 x 30 3.026 2.35 3.31 2.58 1.64 3.09

60 x 30 x 60 0.81 0.61 0.88 0.70 0.46 0.83

120 x 60 x 120 0.20 0.15 0.22 0.18 0.12 0.21

= {х1 е (0, п), х2 е (0, п), х3 е (0, п)} со следующими начальными и граничными условиями:

и t=0 = и=0 ' В|t=0 = Вап| 1 =0' и|дв = Иап|ВО' В1дв = Вап|дв '

Отметим, что нормальные составляющие и и В на границах расчетной области в равны нулю, т.е. выполняется условие непротекания. Расчеты проводились до момента времени Т = 1.5 на равномерных сетках (п1 х п2 х п3), где па — количество контрольных объемов вдоль соответствующей координатной линии. На каждом шаге по времени условие соленоидальности для и и В выполнялось с точностью в = 10-4. Шаг по времени выбирался так, чтобы погрешность аппроксимации по времени была существенно меньше погрешности аппроксимации по пространственным переменным (во всех представленных ниже расчетах

А t = 10-5). Точность численного решения оцени-(б)

2

u an = exp(-0.9tf + exp(-0.6f pan =-^,

Ban = 2>/nU an,

где

f1 =

'2 sin(2x1)cos(2x2)sin(x3) - sin(2x1)sin(2x2)cos(x3) ^ sin(2x1)sin(2x2)cos(x3) - 2cos(2x1)sin(2x2)sin(x3) ч2 cos(2x1)sin(2x2)sin(x3) - 2sin(2x1)cos(2x2)sin(x3)

f2 =

' 2 sin(x1) cos(2x2) sin(x3) - sin(x1) sin(2x2) cos(x3) ^ sin(x1)sin(2x2)cos(x3) - cos(x1)sin(2x2)sin(x3) vcos(x1)sin(2x2)sin(x3) - 2sin(x1)cos(2x2)sin(x3)

В качестве одного из тестов рассматривалась начально-краевая задача для системы уравнений (1)-(4) в параллелепипеде G =

(a)

-8 -4 0 4 8

л,,™™,,™ тз. t Bz)

Рис. 2. Структура течения и изолинии компоненты В3 (Вг) индукции магнитного поля в сечении =п при t = 0.5

(а — численное решение, б — аналитическое решение). ДОКЛАДЫ АКАДЕМИИ НАУК том 464 № 5 2015

528

БЕТЕЛИН и др.

валась по максимальным абсолютной A F и относительной 5F погрешностям:

A F = max |F - Fan|,

G 8(0,T ]

5F = ma^-1 Г F -Fa^ d3-100%,

(0,T] Vg { Fan

G

где

Fan =— [| Fan\d VG = I'd S — объем области G. Vg 1 1

G G G

Результаты тестов на численную сходимость на последовательности вложенных расчетных сеток, приведенные в табл. 1, 2, демонстрируют второй порядок аппроксимации на точном решении. На рис. 1, 2 представлены результаты сравнения

численных расчетов на сетке 120 х 60 х 120 с аналитическим решением (12).

Работа выполнена при поддержке РФФИ (гранты 13—01—12051-офи_м, 15-41-00013-р_урал_а).

СПИСОК ЛИТЕРАТУРЫ

1. Патанкар С. Численные методы решения задач теплообмена и динамики жидкости. М.: Энерго-атомиздат, 1984. 152 с.

2. Issa R.I. Solution on the Implicitly Discretised Fluid Flow Equations by Operator-Splitting // J. Comput. Phys. 1985. V. 61. P. 40-65.

3. Issa R.I., Gosman A.D, Watkins A.P. The Computation of Compressible and Incompressible Recirculating Flows by a Non-Iterative Implicit Scheme // J. Comput. Phys. 1986. V. 62. P. 66-82.

ДОКЛАДЫ АКАДЕМИИ НАУК том 464 № 5 2015

Для дальнейшего прочтения статьи необходимо приобрести полный текст. Статьи высылаются в формате PDF на указанную при оплате почту. Время доставки составляет менее 10 минут. Стоимость одной статьи — 150 рублей.

Показать целиком