Стабилизированное течение электропроводящей жидкости в круглой трубе Текст научной статьи по специальности «Физика»

Вычислительные технологии

Том 10, № 3, 2005

СТАБИЛИЗИРОВАННОЕ ТЕЧЕНИЕ ЭЛЕКТРОПРОВОДЯЩЕЙ ЖИДКОСТИ В

КРУГЛОЙ ТРУБЕ

А. М. Бубенчиков, А. В. Колесникова Томский государственный университет, Россия e-mail: [email protected], [email protected]

A simple numerical algorithm for calculation of a stabilized flow of conductive liquids in a round pipe is presented. A five-point finite difference scheme for two unknown variables and Gauss — Seidel method were used to solve the complete problem. The influence of Ha number on the velocity and induced magnetic field patterns is investigated.

Введение

Стабилизированным по длине канала является такое течение, когда поверхность скоростей перестает меняться от сечения к сечению, линии тока становятся параллельными оси трубы, а осевой градиент давления превращается в постоянную величину. Такие течения еще называют автомодельными, поскольку в этом случае можно исключить из рассмотрения аксиальную координату и тем самым уменьшить число пространственных координат хотя бы на единицу (а при наличии еще и осевой симметрии — на две единицы). При ламинарном стабилизированном течении непроводящей жидкости поверхность скоростей представляет собой параболоид Пуазейля. Если рассматривать проводящую жидкость, а трубу поместить в поперечное магнитное поле и постепенно увеличивать напряженность поля, то поверхность скоростей будет подвергаться такому изменению, что параболоид с вершиной на оси канала будет сжиматься в направлении осевой линии. Если продолжать увеличивать напряженность магнитного поля, то наряду с отмеченной деформацией наблюдается сжатие поверхности скоростей и по оси, параллельной вектору напряженности магнитного поля. Несмотря на простоту ситуации, обусловленную автомодельностью, к настоящему времени не было вычислительной технологии, разработанной в рамках метода конечных разностей, с использованием которой был бы рассчитан параболоид скоростей, сжатый под действием магнитного поля в двух взаимно-перпендикулярных направлениях.

В настоящей работе, опираясь на имеющийся опыт решения задач такого рода [1-3], мы представили простейшую вычислительную процедуру, используемую для расчета неосесимметричного автомодельного течения.

© Институт вычислительных технологий Сибирского отделения Российской академии наук, 2005.

1. Постановка задачи

Рассмотрим стабилизированное течение проводящей жидкости в трубе круглого поперечного сечения. Неизвестными величинами в рассматриваемой задаче являются: поле скоростей V, электрическое напряжение E и магнитная индукция B. Вводим цилиндрическую систему координат (x, r, в) такую, чтобы ось x совпадала с осью трубы. Даже в более общем случае бесконечной цилиндрической трубы произвольной формы поперечного сечения удается показать [4], что решение задачи можно свести к определению двух функций — Vx = U(r, в) и Bx = B(r, в), а электрическое поле из уравнений исключить. Как и в задаче о стабилизированном течении вязкой непроводящей жидкости, в данном случае распределения всех механических и физических величин будут зависеть только от полярных координат r и в. Более того, осевой градиент давления, входящий в соответствующую проекцию уравнения импульсов, сохраняет свою величину вдоль оси трубы, т. е. будет константой. В этих условиях электрический ток вдоль трубы невозможен, так что из закона Ома при V = (U, 0, 0) имеем jx = 0, т. е. индуцированные токи будут располагаться в плоскостях нормальных сечений трубы, а индуцированная составляющая магнитного поля будет направлена по оси канала. Принимая во внимание наличие однородного внешнего поперечного магнитного поля Be = (0, Bo cos в, -Bo sin в) и условие соленоидальности магнитного поля divB = 0, заключаем, что суммарные компоненты магнитного поля будут равны

Bx = B(r, в), Br = Bo cos в, Be = -Bo sin в. (1)

При этом, как следует из уравнения Ампера rotB = j0j, величину B(т',в)/^0 можно рассматривать как "функцию тока" плоского поля векторов j = (jr ,je):

• 1 dB . dB

Jo jr = r^, ^oje = - -qt' (2)

Рассмотрим векторное уравнение Навье — Стокса:

p(V -V)V = -Vp + jV2V + — rotB x B. (3)

Jo

Запишем проекции уравнения (3) на оси цилиндрической системы координат (x,r,в). В рассматриваемых нами условиях, когда V = (U, 0, 0), получим

dp _~ Bo cos в dB Bo sin в dB

Tx = jVu + W; (4)

(5)

др _ В дВ др _ В дВ дг ^о дг' дб ^о дб ' Два последних соотношения можно переписать в виде

д В2 д В2

¥г[р _0, Ге{" +2^ _0-

Полученные равенства говорят о том, что сумма р + В2/2может быть функцией лишь координаты х. Замечая, что В зависит только от г и в, заключаем, что др/дх может зависеть лишь от х. Однако правая часть (4) может быть функцией лишь г и в, но так как левая часть есть осевой градиент др/дх, зависящий лишь от х, каждая из этих частей

в отдельности может быть только постоянной величиной. Представим эту постоянную в виде

dp Ap

дХ = "17' (6)

где l0 — характерный масштаб длины, а Ap — величина продольного перепада давления на этой длине. Таким образом, уравнение (4), в левой части которого стоит константа вида (6), является первым из двух искомых уравнений. Второе уравнение для В(т,9) можно получить из уравнения индукции следующего вида [5]:

— —У2В = arot(V х В). (7)

Уо

Учитывая, что в рассматриваемом случае В = (В(т,в),В0 cos 9, —В0 sin 9) и V = (U(т,9), 0,0) при проектировании (7) на ось x, найдем

1 ^л ( ndU sin9 0U\

--V В = аВо( cos 9----— . (8)

у0 \ дт т д9 J

Две другие проекции (7) на оси т и 9 удовлетворяются автоматически.

Если в качестве характерных масштабов длины, скорости и магнитной индукции выбрать величины

- Ap / Q \ 1/2

т * = /о, U* = -0-, В* = loAp^o\ л , (а)

у W

то безразмерный вид уравнений (4) и (8) будет следующим:

1 d2U 1 д ( dU\ тт ( пдВ sin9дВ\

+ т— + Яа cos9----— = —1; (9)

т2 дв2 т дт \ дт J \ дт т дв

1 д2В 1 д ( дВ\ ( BU sin9дU\ , ч

+ т— + Яа cos9----— =0. (10)

r2 д92 r dr \ dr) \ dr r 39

В записи уравнений (9), (10) и ниже все величины уже безразмерные. Однако для них оставлены прежние обозначения.

Эти уравнения будем решать при граничных условиях [4]:

r = 1, U = 0 (11)

— условие прилипания;

dB B

r =1, тг + - = 0 (12)

dr p

— условие сопряжения магнитного поля на границе двух областей. Здесь p = ——, а\, а2 —

a\R

проводимость жидкости и материала стенки соответственно, 8 — толщина стенки, R — радиус канала.

Линейными преобразованиями из уравнений (9) и (10) нетрудно получить следующую систему:

1 d2(U + B) + 1 в (rS(U + B) • + Ha / 9d(U + B) _ sin*d(U + B) •. = _1,

т2 д92 т дт \ дт ) \ дт т д9

1 8Чи - В) + 1 а (Г8(и-В^] - со8 /хи-В) - втв Щ- В)

дв2

г дг

дг

дг

в

Вводя обозначения Ш — и + В и X — и — В, эту систему можно переписать в виде

= -1, (13)

1 д2Ш 1 д ( дШ\ тт ( Ш втвдШ

+ --гг- г+ ИМ сое в

г2 в2 г г

г

г

г дв

г2 в2 г г г г г в

-1.

(14)

2. Решение задачи

Аппроксимируем дифференциальные члены, входящие в (13) и (14), следующими разно-

стями:

сое в

дШ г

1 д2Ш _ - + + о(Дв2),

в2

г2Дв2

1 д ( дШ

г дг \ дг

_ гг+1/2Шг+Ц - (гг+1/2 + гг-1/2)Шг,з + гг-1/2Шг-1,з + 0(Дг2)

ггДг2 '

вШ в дШ (вШ вз + | вШ вз |) Шг,з - Шг3-1

г дв

2гг

(яп вз - | яп вз |) - Ш^з

2гг

(С08 вз + | сов вз |) Шг+1'' - Ш

Дв

Дв + о(Дв),

2Дг

- + (сов вз - | сов вз |) Ш -Шг-1* + о(Дг).

2Дг

Подставляя данные представления в (13), получим

Шг

,з _ (аг,з-1Шг,з-1 + аг,з+1Шг,з+1 + аг-1,зШг-1,з + аг+1,зШг+1,з + 1)/а

(15)

где

аг,з-1

Чз'+1

1 „ | э1п вз | + э1п вз

(ггДв)2

+ Иа1

2ггДв

1 „ | э1п вз | - э1п вз

(ггДв)2

+ Иа-

2ггДв

аг-1,з _

гг-1

г 2

аг+1:

з ггДг2

_

з ггДг2

+ Иа

+ Иа

| сов вз | - сов вз

2Дг

| сов вз | + сов вз 2Дг

аг,з — аг,з-1 + аг,з+1 + аг-1,з + аг+1,з • Формулу для вычисления получаем аналогично, и вид ее будет схожим с (15):

,з (сг,3-1 Хг ,з-1 + сг,з+1 Хг ,з+1 + сг-1,з Хг -1,з + сг+1,з Хг+1,з + 1)/сг,з •

(16)

(17)

При этом Сг,з-1 — аг,з+1, Сг,з+1 — аг,з-1, сг-1,з — аг+1,з, сг+1,з — аг-1,з, сг,з — аг,з, г — 1...М - 1,

3 — 1 ••• М - 1; ^о — ^м, Шг,о — Шг,ы, г — 1 ...М - 1•

г

г

Уравнения (15), (17) используются для нахождения значений Ш и Z во внутренних узлах. Разностная реализация условий (12) на границе с порядком о(Дг4) дает

WN.

p(9(WN-!j - Zn-i,j) - 9(Wn-2,j - Zn-2,j) + Wn-3- Zn

j 11p + 6дг

ZNjj = -WNjj, j = 0 ...M. (18)

По построению коэффициенты разностной схемы удовлетворяют критерию Скарбор-роу [8], поэтому итерационный процесс Гаусса — Зейделя, состоящий в простом пересчете значений W и Z по формулам (15), (16) и (18), является сходящимся.

Как видно из выражений (16), разностные уравнения (15) и (17) не могут быть использованы при i = 0, так как r = 0. Поэтому следует получить отличные от (15), (17) формулы, которые можно было бы использовать для определения W0j, Z0j. Более того, не только (15), (17), но и их дифференциальные прообразы — уравнения (13), (14) — по тем же соображениям не могут быть применены непосредственно в точке r = 0. Если записать определяющие уравнения в других координатах, то на их основе можно получить упомянутые формулы. Запишем (13), (14) в декартовых координатах, направив ось у по вектору напряженности внешнего однородного магнитного поля:

V2W = -(Ha^ + Л , V2Z = Ha^ - 1, (19)

V dy J dy

д2 д2

где V2 = + ; Ha — число Гартмана. дх2 ду2

Рассмотрим поведение решений этих уравнений в малой окрестности точки r = 0. Пусть r < h, где h — шаг разностной сетки по радиальной координате. Для примера рассмотрим первое из уравнений (10). Его правую часть обозначим через f. Положим, что в этой окрестности значение f = f0 — величине правой части на осевой линии канала (при r = 0). Тогда в указанной окрестности рассматриваемое уравнение приближенно можно аппроксимировать следующим образом:

VW = f0, (20)

где f0 = const. Если теперь положить

W = Wi + W2, (21)

где

f0 f0

Wi = ^(х2 + у2) = ^4 r2 (22)

есть частное решение (20), то для определения W2 будем иметь уравнение Лапласа

AW2 = 0. (23)

Проинтегрируем (21) по длине lh элементарной окружности радиуса h с центром в начале координат. С учетом (22) получим

1 Г 1 Г f0h2

1 Wdlh = — ф W2dlh + J—. (24)

2nh J 2nhJ 4

lh lh

С другой стороны, записывая (21) в точке г — 0, получим

Шо — Ш2о. (25)

Далее, так как Ш2 — гармоническая функция, по теореме о среднем для потенциальных функций первый член правой части (24) равен значению функции Ш2 в нуле

1 ' Ш2(ИН — Ш2о. (26)

2пЬ

Тогда из (24)-(26) находим

1С №

Шо - }-Т- (27)

Первое слагаемое в (27) имеет о(1), а второе — о(Ь2). В рассмотренных ниже примерах влияние второго слагаемого ничтожно, поэтому в дальнейшем мы будем работать без него, простейшая же аппроксимация первого слагаемого в (27) (интеграла) дает

1 М

Фо — мТ, Ф 1,к, Ф — (28)

к=о

Следовательно, при стабилизированном движении значение всех магнитогидродинами-ческих параметров на оси трубы есть среднее арифметическое соответствующих значений на ближайшей к оси сеточной окружности.

В работе [1] обоснование формулы (28) получено способом, отличным от представленного выше. Кроме того, в отмеченной работе найдены соотношения для расчетных величин на оси канала, обеспечивающие более высокий порядок аппроксимации относительно шага по радиусу.

Таким образом, решение исходной дифференциальной задачи (9)-(12) эквивалентно решению двух систем разностных уравнений для Ш и X, определяемых соотношениями (15)-(18), (28). Как видим, при р — 0 система для определения Ш никак не связана с аналогичной системой разностных уравнений для X, а с ростом р увеличивается взаимная обусловленность соответствующих систем уравнений. Поэтому безразмерный комплекс ( , характеризующий проводимость стенки канала, можно рассматривать как параметр связности получающихся систем разностных уравнений. Как оказалось, от этого параметра существенным образом зависит сходимость итерационного процесса. Как показывает опыт расчетов, сходимость к некоторому стационарному распределению наблюдается при любом соотношении шагов по г и в, но не при любом значении р. Так, при р — 0 наблюдается быстрая сходимость к точному решению, а с ростом р она замедляется или не достигается

В

вовсе. При аппроксимации производной ——, входящей в (12) с первым порядком точности

г

относительно шага Дг, получен сходящийся вычислительный процесс лишь при р — 0.01. Если указанную производную в (12) аппроксимировать с порядком о(Дг4), как и записано в (17), то в вопросе сходимости можно продвинуться до значения р — 0.1.

В процессе вычислений использованы формулы (15)—(18), (28), которые в своей совокупности обеспечивают нахождение всех значений Шг,з и , г — 0...М, з — 0...М, если заданы их исходные распределения в виде какого-либо нулевого приближения.

3. Результаты расчетов

Предварительно были выполнены тестовые проверки, и при На = 0 получены параболоид Пуазейля и тривиальное распределение для второй искомой функции: В = 0.

На рис. 1, 2 показаны поверхности скоростей и индуцированного магнитного поля, рассчитанные при На = 10. Как видим из рис. 2, по диаметру, параллельному силовым линиям поля, профиль скорости перестраивается качественно подобно тому, как это происходит в условиях простейшей магнитогидродинамической задачи для плоского канала. Поверхность В = В (у, г) является антисимметричной относительно оси г, перпендикулярной к Ве (см. рис. 1). Далее рассчитывались собственно магнитогидродинамические

Рис. 1. Поверхность, определяющая инду- Рис. 2. Поверхность скорости, отвечающая цированное магнитное поле, полученная при случаю На =10, у = 0.001. На = 10, у = 0.001.

а б

Рис. 3. Профили скоростей, полученные при На = 4 в двух сечениях: а — параллельном вектору напряженности внешнего магнитного поля (зависимость от у), б — перпендикулярном к этому вектору (зависимость от г); сплошные кривые — расчет, пунктир — решение [4].

течения, т. е. различные варианты, отвечающие условию На = 0. На рис. 3 представлены результаты вычислений, выполненных по описанному выше алгоритму (сплошная линия), и решение [6], найденное с использованием метода Фурье (пунктирная линия). Как видим, имеется хорошее согласие результатов аналитического и численного решений.

Как видим из рис. 4, а, характерное значение напряженности магнитного поля содержит величину перепада давления. В то же время известно, что при движении проводящей жидкости в магнитном поле коэффициент трения Л как безразмерный перепад давления прямо пропорционален числу Гартмана, т. е. безразмерной величине напряженности поля. Поэтому рассчитываемое безразмерное значение В можно считать уже отнесенным к напряженности поля. В связи с этим не вполне очевидным представляется характер изменений в распределении В с ростом или уменьшением напряженности поля, т.е. чисел На.

Как показали расчеты, представленные на рис. 4, с ростом значений На экстремум в распределении В смещается к периферии, т.е. к стенке канала. Однако если На < 5,

а б

Рис. 4. Профили потенциала, полученные при р = 0.01 в двух сечениях: а — параллельном вектору напряженности внешнего магнитного поля (зависимость от у), б — перпендикулярном к этому вектору (зависимость от г); На = 2, 2.5, 5,., 6.5, 8,10 (кривые 1 — 7 соответственно).

В 0.06

0.04

0.02 12 3

|| 1 1 1 1 1 1 1 1 1 1 -0.5 / 0 0.5 1 2

[1 -0ГО2

\\\ / 0 04

- 0.06

Рис. 5. Влияние р на распределение индуцированной составляющей магнитного поля при На = 10: р = 0, 0.05, 0.1 (кривые 1 — 3 соответственно).

то величина экстремума возрастает с увеличением числа Гартмана, а при Ha > 5 картина меняется на противоположную. На рис. 5 показано влияние параметра p, характеризующего проводимость стенки канала, на распределение величины индуцированного магнитного поля B. Как видно, с увеличением p экстремум в распределении B растет и его положение смещается к периферии, но распределения относительной малой величины скорости U/U, полученные при различных p и одном и том же числе Гартмана, остаются неизменными.

Заключение

Как показал проведенный анализ, математическое описание внутренних стационарных автомодельных течений электропроводящей жидкости допускает разделение на квазисимметричные части на уровне определяющих уравнений. Более того, для случая непроводящих стенок p = 0 удается разделить задачу и на уровне граничных условий. С увеличением проводящих свойств стенки (с ростом неотрицательных значений p) увеличивается степень связанности частей задачи через граничные условия. Применение метода конечных разностей с симметрической аппроксимацией дифференциальных членов и монотонной аппроксимацией всех других дифференциальных слагаемых уравнений позволяет получить системы алгебраических уравнений, которые эффективно разрешаются с помощью простейшей итерационной процедуры Гаусса — Зейделя, по крайней мере для малых значений p. В результате численного решения задачи для случаев высоких величин напряженности магнитного поля подтвержден эффект сжатия поверхности скоростей по направлению силовых линий поля. Проведен также параметрический анализ влияния чисел Гартмана на характер распределения индуцированной составляющей вектора магнитной индукции.

Список литературы

[1] Бубенчиков А.М., Ливаев Р.З. Некоторые автомодельные задачи магнитной гидродинамики // Вест. ТГУ. Бюл. оперативной науч. информации. 2001. № 4. С. 32-52.

[2] Бубенчиков А.М., Клевцова А.В., Фирсов Д.К. Течение электропроводящих жидкостей в тонких трубках в поперечном магнитном поле // Мат. моделирование. 2003. Т. 15, № 9. С. 76-87.

[3] Бубенчиков А.М., Клевцова А.В., Харламов С.Н. Закрученный поток проводящей жидкости в узких трубах при наличии магнитного поля // Мат. моделирование. 2004. Т. 16. С. 109-122.

[4] Лойцянский Л.Г. Механика жидкости и газа. М.: Наука, 1973. 848 с.

[5] Тананаев А.В. Течения в каналах МГД-устройств. М.: Энергоатомиздат, 1981. 420 с.

[6] Владимиров В.С., Жаринов В.В. Уравнения математической физики. М.: Физмат, 2003. 400 с.

[7] Magnetohydrodynamic pipe flow. Pt 1 // J. Fluid Mech. 1962. Vol. 13, N 4. P. 505-512.

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

Поступила в редакцию 29 ноября 2004 г., в переработанном виде — 27 декабря 2004 г.