Исследование линейной устойчивости адвективного течения методом сеток Текст научной статьи по специальности «Математика»

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

Том 6, № 6, 2001

ИССЛЕДОВАНИЕ ЛИНЕЙНОЙ УСТОЙЧИВОСТИ АДВЕКТИВНОГО ТЕЧЕНИЯ МЕТОДОМ СЕТОК

Е.Л. Тлрунин, К. Г. Шварц Пермский государственный университет, Россия e-mail: tarunin@psu.ru

The problem of the linear stability theory for advective flow in horizontal rotating liquid layer is solved by reduction to the boundary problem for the system of linear onedimensional partial differential equations. The distinctions of grid method realization are described for the given problem. The calculated neutral curves for normal perturbations define the dependence on the critical Grashoff number by wave number for the different Taylor number and Pr = 6.7 (water).

Исследование спектра малых нормальных возмущений конвективных течений и их линейной устойчивости сводится, как правило, к задаче на собственные значения для системы обыкновенных дифференциальных уравнений высокого порядка с переменными коэффициентами и малыми параметрами при старших производных [1]. Наиболее распространенными численными методами решения таких задач являются метод Галеркина, метод пошагового интегрирования с ортогонализацией и метод дифференциальной прогонки [2]. В данной работе задача линейной теории устойчивости решается путем сведения ее к начально-краевой задаче для системы линейных одномерных уравнений в частных производных, которая решается затем методом сеток. Аналогичная методика поиска критических чисел Рэлея для замкнутых областей использована в [3]. Возможности метода демонстрируются на примере построения нейтральных кривых для адвективного течения, возникающего во вращающемся горизонтальном слое жидкости [4]. Метод легко реализуется и, кроме того, дает полезную информацию о характере установления надкритических движений. Тому, кто освоил метод сеток для решения нелинейных задач, для построения нейтральных кривых достаточно убрать из уравнений нелинейность и зависимость от “излишних” координат.

Адвективное течение, возникающее во вращающемся плоском горизонтальном слое жидкости, в котором ось вращения совпадает с вертикальной осью координат, для случая, когда температура на границах линейно изменяется с продольной координатой, впервые описано аналитически С. Н. Аристовым [5]. В дальнейшем было построено целое семейство точных решений при различных граничных условиях. Особенностью этих течений является отсутствие вертикальной компоненты скорости. Вектор скорости в потоке ориентирован перпендикулярно к силе плавучести (имеются обе его горизонтальные компоненты). Адвективное течение во вращающемся слое имеет экмановскую спиралевидную форму [6] и наряду с изотермическим течением Экмана является основным пограничным течением в

© Е.Л. Тарунин, К.Г. Шварц, 2001.

атмосфере и океане. Работ по исследованию устойчивости адвективного течения во вращающемся горизонтальном слое жидкости до последнего времени не было. Первая работа в этой области [4], к сожалению, не лишена неточностей, которые здесь исправлены.

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

Рассмотрим плоский горизонтальный слой несжимаемой жидкости, ограниченный твердыми плоскостями Z = ±h и вращающийся с постоянной угловой скоростью Qo- Направление оси вращения совпадает с вертикальной осью координат. На обеих границах заданы температура, линейно меняющаяся с координатой х, условие прилипания и замкнутости потока

h

Z = ±h : T = Ax, A = const, v = 0, J vxdz = 0. (1)

—h

Возникающее адвективное течение описывается аналитически [3]:

u0(z) = ~^= v0(z) = ^^[z - Re/i(z)L

V Ta v Ta

z

GrPr f

To = x + GrPr To(z), po = x(z + 1) +-------vo(Z )d(,

Ta

—i

T0(z ) = vo(z), /i(z ) = |h£, A = ^ ^ <2)

где z = Z/h; uo, vo — горизонтальные компоненты вектора скорости, To — температура;

v „ „ gßAh4 / 2Qoh24 2

р0 — давление; Рг =-----число Прандтля; Ог =------------число Грасгофа; Та

X ^ V V

— число Тейлора; Ле, 1т — действительная и мнимая части комплексного значения; V —

кинематическая вязкость жидкости; X — коэффициент теплопроводности; д — ускорение

свободного падения; в — коэффициент теплового расширения [4].

Линеаризованные уравнения для возмущений скорости V, температуры 0, давления Р в безразмерном виде имеют следующее представление во вращающейся системе координат:

«V

— + Ог [^У)уо + + ^(Ц х V) = -УР + ДV + 0^, divV = 0, (3)

ЛЛ -1

— + Ог [(VV)To + (уоУ)0] = —Д0, (4)

г = ±1 : V = 0, 0 = 0; уо = (м0,г>0, 0), V = (и, г>,эд), 1z = (0,0,1). (5)

В силу большой сложности исследования устойчивости течения в трехмерной постановке будем рассматривать предельные случаи: пространственные винтовые периодические по х возмущения в виде валов с осью, перпендикулярной к оси X, и пространственные спиральные периодические по у возмущения в виде валов с осью, параллельной оси X. Заметим, что при отсутствии вращения для Та = 0 винтовые возмущения вырождаются в плоские.

2. Винтовые возмущения

Уравнения винтовых возмущений выводятся из системы (3) - (5) в предположении, что производные по у от всех функций равны нулю. Имеются все три компоненты вектора возмущения скорости, которые так же, как и возмущения температуры, являются функциями времени £ и двух пространственных переменных х, г. Учитывая дивергентность возмущений скорости, введем функцию тока возмущений

и — —

—0

dz ’

w —

—0

дж

и вихря возмущения скорости

р

du dw / d 20 д20

dz дх \ дх2 + dz2

(6)

(7)

Представим неизвестные функции в виде

р — [<£i(i, z) + ip2(t, z)] exp 0 — [01 (t, z) + i02(t, z)] exp

v — [v1(t,z) + iv2(t,z)]exp ikxx, 9 — [01(t,z) + i92(t,z)]exp ikxx, (8)

где kx — волновое число. В результате задача сведется к системе линейных уравнений в частных производных по времени t и переменной z:

dp

dt

kxGr[uo(z)p2 + Uq;02] — VTa —— — P— k^pi + kx02

dz —z2

д 20i

—z2

;— kX 01 + Pi — 0,

—p2 , i n Г ^ і // / 1 —v2 —2 p

— + kxGr[uo(z)p1 + Uo01] — VTa — —

д202 k2 , + — 0

— kx02 + P2 — 0,

22

— кжР2 — kx9b

dv

dt

dv

1 — kxGr[uo(z)v2 + v'002] — VTa—01 — — kX vb

dz dz2

dt

2 + kxGr[uo(z)v1 + vg 01] — VTa—02 — —— k^ v2,

dz

—91 — kxGr[uo(z)92 + GrPrrg 02] — Gr—01 — -1-dt dz Pr

dz2 д291

dt

с граничными условиями

01,2

д9>2 + kXGr[uo (z)91 + GrPrrg 01] — Gr—02

dz

1

Pr

dz2 d2 9

— k? 91

dz2

2 — kX92

—0

1,2

z=±1

dz

= v1,2 ,2 91 =0

z=±1 z=±1 z=±1

(9)

(10)

Одномерная система уравнений (9) с граничными условиями (10) решается по численной методике, аналогичной схемам двухполевого метода [7], который используется для

решения двумерных задач в переменных функции тока и вихря скорости. В нашем случае функция тока и вихрь скорости зависят от одной пространственной переменной поперек слоя. Вторая переменная исчезла в предположении периодичности решения от этой координаты. Уравнения для возмущения вихря и функции тока решаются с помощью классической неявной схемы:

1 »“+1 _ /^А + 1 + кЛ »“+1 + 1 »“+1 =

»1’;+1 V + т + »1>; + »1’;-1

1 «“ —

= -- к,с.г[«о(-,к, + <(-,ж,,] - ^та 1,;+12^^1’^-1 - к,^

(0 < л' < N), й,»1 = ^?,с »“+ = »“,«,

Я+« “ (/? + к,) «+1 + /12«+- = -<+1 (0 < Л < N),

«У = 0, С+1 =0,

1 »“+1 _ ГА + 1 + кЛ »“+1 + А»“+1

/2 »2’;+1 V + т + V ; ,;-1

1 «“ «“

= -^ + №[«„(г,-)»“+1 + <(2;)«+‘] - ^ **» 2,-1 + к,«;,

(0 < л < N),

»“+1 = »“,о> »“,+/ = »“.у,

/2; - (/2 + к,) «+1 +1-С-. = -»“+1 (0 < л < N),

^’Й1 = 0, </й‘ = 0, (11)

где т — шаг по времени; п = 0,1,... — номер шага по времени; = 2/N — шаг по вертикальной координате; N — число интервалов разбиения. Вихри р1, »2 на границах аппроксимируются по формуле Вудса [4]:

3 1 3 1

»’К1 = -щ«+1 - 2^¡51, »“+ = -щ«.Л- - 2¥>“+?-1 (к = 1,2). (12)

Остальные уравнения системы (9), (10) решались с помощью схемы Кранка — Николсона

[8], погрешность аппроксимации порядка 0(т2 + ^2):

1 “+1 ( 1 1 к,\ “+1 1 “+1 1 “ «“ , + 1 - 2«“ ; + V“ ;-1 #2 ;

2/1 ; ч Л + 7 + V«“+1 + 2й! , = - 7й“.; - 2/2 + т -

^“+1 _ ^“-1 -к,С,г[«о(2,К, + V»(г;)«+1] - х/тТ'Х;+12/^ ’;-1

(0 < л < N),

«’.'л1 = 0, «“# = 0,

vn+1 2h2 V2’j+1

1 1 kX

---------1-----------+ —

h2 + t + 2

1

,,n+1

V j1 + 2h2 V2j-1

1 ..n j - ^j + j + k2

t"2’ j 2h2 + 2

x vn + v2,j +

^n+1 _ rn-1

+kxGr[uo(Zj)<j + "0(Zj)rn+1] - ^2 ’j-1

2hz

(0 < j < N),

vn+1 .2 o

0, Л1 = 0,

1

0П++1 -

2Prh2 1 ’j+1

11 kX

----1---1-—

Prh2 t 2Pr

-1 -1 ЛП _______ O/^n I Z)n

fln+b 1 fln+^ J-an P1 , j+1 2P1 j + P1 j-1

#1 ’j + 2Prh2 ^ ’j-1 = - т 1 ’j

2

2Pr

9n j - kxGr[uo(Zj)^nj + PrGrTO(Zj)^+1] - Gr

2Prh2

n- 1

j-1

2h,

(0 < j < N),

ОД1 = 0,

0,

1

9n+1 -

^9 o_L1

2Prh2 2’j+1

11 kX

----1---1-—

Prh2 t 2Pr

1

0n’+1 + 2Prh2U2 ’ j-1

i fin Ofin i лп

лп+1 _ 1 rm *2 ’ j+1 2^2 ’ j + ^2 , j-1

2 j-1 t 2 ’j 2Prh2

2Pr

j + kxGr[uo(Zj)0?+1 + PrGrTO(Zj)^+1] - Gr

qijn+1 _ ?/;n-1

r 2 ’7 + 1 r 2 ’7-

2 j+1

2 j-1

2hz

(0 < j < N),

0n+1 = 0, tfnN1 = 0.

(13)

В качестве начальных возмущений для неизвестных 01? 02, v1? v2, 01, 02 берутся функции, пропорциональные sin2 nz, удовлетворяющие граничным условиям (10). Число узлов сетки, определяемое шириной экмановского пограничного слоя и точностью расчетов, изменялось от 101 до 201. Расчеты проводились при фиксированном Pr = 6.7 (вода). Число Тейлора менялось от 0 до 106.

Для построения нейтральной кривой, описывающей зависимость критического числа Грасгофа от волнового числа kx при фиксированном числе Тейлора, для каждого выбранного значения kx требуется найти такое число Грасгофа, при котором действительная часть декремента возмущений Л = А1 + iA2 равна нулю. Иными словами, решается задача о поиске корня А1 = 0 для неявной функции A1(kx, Gr, Ta). Мы строим эту функцию по точкам (дискретно) с помощью многократного решения эволюционной задачи (9), (10) методом сеток.

Для нахождения действительной части декремента возмущений Л1 прослеживалась эволюция во времени максимумов по модулю неизвестных. В качестве аппроксимации зависимости амплитуд от времени использовалась экспоненциальная формула C exp A1t. Неизвестные Си А1 определялись методом наименьших квадратов [9] по ходу вычислений уравнений (11) — (13). Для определения неизвестных необходимо производить расчеты до 20-25 единиц безразмерного времени при шаге по времени порядка h^, где hz — шаг по вертикальной координате. Использованный алгоритм экономичен и позволяет рассчитывать коэффициенты аппроксимации по ходу основных расчетов. Нулевое значение декремента возмущений уточняется методом хорд. Характер поведения возмущений во времени существенно зависит от всех параметров задачи; в области неустойчивости (Л1 > 0) все возмущения нарастают, а в области устойчивости (Л1 < 0) затухают. Судить о характере

X

f

k

устойчивости достаточно по какой-нибудь одной амплитуде. Опыт решения этой задачи показал, что наиболее информативной является амплитуда возмущения температуры.

Расчеты показали, что для винтовых возмущений результаты работы [4] были ошибочны. Причиной ошибки явился неверный знак, стоящий перед одним из слагаемых в разностной схеме. При заданном числе Прандтля наиболее опасной является колебательная неустойчивость, а не монотонная. Декремент возмущений Л в этом случае является комплексным числом. Тестовые расчеты, проведенные для случая отсутствия вращения (Та = 0), показали, что критическое число Грасгофа Сгк ^ 126.7(1 — 31.бЛ^). Из этой оценки видно, что относительная погрешность вычисления критического числа Грасгофа составляет менее 0.32% при Нг < 0.01.

Типичные зависимости амплитуды возмущений температуры в1 от времени для двух значений числа Грасгофа при фиксированных значениях числа Тейлора (Та = 1000) и волнового числа кх = 5.5 представлены на рис. 1. Первая кривая (Сг = 700) соответствует области колебательной неустойчивости с Л ~ 1.644, вторая (Сг = 630) — области устойчивости с Л1 ~ —0.379. Граница устойчивости (Л1 ~ 0) соответствует Сг = 643.5.

20

-20

. 1 А

, к I 2 '

0.0

0.5

1.0

Рис. 1. Эволюция температуры 6\ во времени при Та = 1000, к = 5.5: 1 — Сг = 700, 2 — Сг = 630.

Рис. 2. Нейтральные кривые для плоских возмущений при различных числах Тейлора (а); зависимость мнимой части декремента Л2 от Ta (б).

Мнимая часть декремента возмущения Л2 определялась после построения нейтральной кривой. Для этого выбиралось критическое число Грасгофа и соответствующее ему волновое число и прослеживалась эволюция во времени какой-либо неизвестной (напри-

мер, ui) при фиксированном значении аргумента z (например, вблизи начала пограничного слоя). В качестве аппроксимации зависимости от времени использовалась формула C sin A2t. Неизвестная C определялась по максимуму модуля прослеживаемой функции, а период — методом наименьших квадратов.

На рис. 2, а представлены нейтральные кривые зависимости критического числа Грасго-фа от волнового числа при различных значениях числа Тейлора. Как видно, вращение повышает устойчивость адвективного течения. С увеличением Та растет критическое волновое число kx, соответствующее минимуму критического числа Грасгофа. На рис. 2, б представлена зависимость мнимой части декремента возмущений от числа Тейлора при Gr = Grk. Как видно, А2 убывает при 0 < Та < 200. Для значений числа Тейлора больше 200 частота колебаний возмущений возрастает с увеличением Та.

3. Спиральные возмущения

Уравнения спиральных возмущений выводятся из системы (3) - (5) в предположении, что производные по х от всех функций равны нулю. Имеются все три компоненты вектора возмущения скорости, которые так же, как и возмущения температуры, являются функциями времени £ и двух пространственных переменных у, г. Учитывая дивергентность возмущений скорости, введем функцию тока возмущений

дФ

v = —z—, w

и вихря возмущения скорости

Ф

Sv Sw

Sz Sy

-ДФ = -

дФ

dy

д 2Ф д 2Ф

+

dy2

dz2

(14)

(15)

Представим все неизвестные функции в виде

Ф = [Фі(t, z) + гФ2(і, z)] exp ikyy, Ф = [Фі(t, z) + гФ2(і, z)] exp ikyy,

u = [u1 (t, z) + iu2(t,z)]exp ikyy, 9 = [91(t,z) + i92(t,z)]exp ikyy,

(16)

где ку — волновое число. В результате задача сведется к системе линейных уравнений в частных производных по времени £ и переменной г:

— kyGr[v0(z)^ + vo Ф2] — ^Та^ куФ1 + ky92 ■>

д 2Ф

dt -А^2 + vo Ф2^^ ^dz dz2 -У

-д~г — + Фі = 0,

дФ2 + ky Gr[vo(z^i + v0,Фl] — VTadU2 = — k^2 — ky 9і,

dt

дz дz2

&

^ — &2Ф2 + Ф2 =

иоФ2]

диі — ky Gr[vo(z)u2 + и0Ф2] + VTaдФі = — k^ui,

дz дz2

d2u2

dt

5^i

"dt~

<902

dt '

с граничными условиями

— ky Gr[vo(z)02 + GrPrrO Ф2] — Gru1 + ky Gr[v0 (z)01 + GrPrr) Ф1 ] + Gru2

dz dz2

1

Pr 1

Pr

— kyU2,

^ k2fl'

— kyfl1 ^ k2fl'

— kyfl2

(17)

Ф1,:

дФ1,;

z=±1

dz

— «1,2 — 01,2 — 0. (18)

z=±1 z=±1 z=±1

Система уравнений (17) с граничными условиями (18) решается с помощью неявных разностных схем аналогично случаю винтовых возмущений. В качестве начальных возмущений для неизвестных Ф1, Ф2, «1, u2, 01, 02 также берутся функции, пропорциональные sin2 nz, удовлетворяющие граничным условиям (18). Число Тейлора менялось в пределах от 0 до 40 000.

Расчеты подтвердили результаты работы [4]: при спиральных возмущениях возможна колебательная неустойчивость.

На рис. 3, а представлены нейтральные кривые зависимости критического числа Грасго-фа от волнового числа ky при различных значениях числа Тейлора. Заметно, что вращение существенно повышает устойчивость адвективного течения. На рис. 3, б изображена зависимость “критической” мнимой части декремента Л2 от числа Тейлора (Л2 при этом соответствует минимуму нейтральной кривой Gr(ky)). Из представленной зависимости видно, что с ростом числа Тейлора Л2 растет и, следовательно, увеличивается частота колебаний возмущений. С ростом Та критическое волновое число увеличивается. Этот факт необходимо учитывать при выборе шага по времени при расчетах методом сеток задачи устойчивости.

Gr 30 000 ■

25 000 -

20 000 -

15 000 -

10000-

5000-

0

s

1

Та = 40 000

\ \ /

\Ч Та = 10 000

Та = 5000

^2

600-

400

б

200

10

15

7 /

/

Та

0 10000 20000 30000 40000

Рис. 3. Нейтральные кривые для спиральных возмущений при различных числах Ta (а), зависимость мнимой части декремента от числа Тейлора (б).

2

Анализ нейтральных кривых для обоих случаев показал, что при малых значениях числа Тейлора наиболее опасными являются спиральные возмущения, а при Та > 60 — винтовые пространственные возмущения (рис. 4).

Рис. 4. Зависимость критического числа Грасгофа от числа Тейлора: 1 — для случая спиральных возмущений, 2 — для случая винтовых возмущений.

Выводы

1. Показана работоспособность представленного алгоритма для исследования линейной устойчивости адвективных течений при наличии вращения. Использовать метод особенно просто тем, кто освоил метод сеток.

2. Выявлена погрешность метода.

3. Исправлены результаты расчетов [4] и получены новые данные о границах устойчивости винтовых возмущений.

4. Построены нейтральные кривые для винтовых и спиральных возмущений при различных значениях числа Тейлора.

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

[1] Гершуни Г. З., Жуховицкий Е.М., Непомнящий А. А. Устойчивость конвективных течений. М.: Наука, 1989.

[2] Гольдштик М. А., Штерн В.Н. Гидродинамическая устойчивость и турбулентность. Новосибирск: Наука, 1977.

[3] Тарунин Е. Л. Определение границы конвективной устойчивости равновесия методом сеток // Алгоритмы и программы для ЭВМ. М., 1978. Деп. в ВИНИТИ 31 мая 1978, №1803-78(78).

[4] Аристов С. Н., Шварц К. Г. Об устойчивости адвективного течения во вращающемся горизонтальном слое жидкости // Изв. РАН. МЖГ. 1999. Т. 4. С. 3-11.

[5] Аристов С. Н., Зимин В. Д. Адвективные волны во вращающемся шаровом слое. Свердловск, 1986 (Препр. / АН СССР. Уральский научный центр. Ин-т механики сплошных сред).

[6] Аристов С.Н., Фрик П. Г. Динамика крупномасштабных течений в тонких слоях жидкости. Свердловск, 1987 (Препр. / УрО АН СССР. ИМСС).

[7] Тарунин Е. Л. Вычислительный эксперимент в задачах свободной конвекции. Иркутск: Изд-во Иркутского ун-та, 1990.

[8] Самарский А. А. Теория разностных схем. М.: Наука, 1989.

[9] Демидович В. П., Марон И. А., Шувалова Э.З. Численные методы анализа. Приближение функций, дифференциальные и интегральные уравнения. М.: Наука, 1967.

[10] Гершуни Г. З., ЖуховицкиЙ Е. М. Конвективная устойчивость несжимаемой жидкости. М.: Наука, 1972.

Поступила в редакцию 28 декабря 2000 г. в переработанном виде — 14 сентября 2001 г.