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

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

Том 5, № 6, 2000

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

И. Б. Пллымский Новосибирский военный институт, Россия e-mail: nina@nsc.ru

A spectrum method for the calculation of stochastic, convective currents of incompressible fluid in a rectangular domain that appear when heated from below under supercritical condition equal to 1000 is suggested. The spectral characteristics of the differential problem and numerical method are compared.

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

Записанная в отклонениях от равновесного решения исходная система уравнений после обезразмеривания имеет вид [1]

и + Pr (^yШх - ) = + RaQx,

Д^ = —и,

Яг + р(^уЯх - ^хЯу) = ра я - р(1)

где (р — функция тока; ш — вихрь; Я — отклонение температуры от равновесного профиля

(полная температура равна 1 — у + Я); А/ = /хх + /уу — оператор Лапласа; И,а = —-Я —

число Рэлея; Рг = и/х — число Прандтля; д — ускорение силы тяжести; @,у,х — коэффициенты теплового расширения, кинематической вязкости и температуропроводности соответственно; Н — толщина слоя и ¿Я — разность температур на горизонтальных границах.

Система (1) решается в области П = {(х,у)|0 < х < I, 0 < у < 1} с условиями <р = ш = Я = 0 при у = 0,1, 0 < х < I на горизонтальных границах и = шх = Я = 0 при х = 0,I, 0 < у < 1 на вертикальных границах.

Подобные задачи рассматривались многими авторами [1-5]. Как правило, система (1) решалась с периодическими граничными условиями с помощью спектрального метода. Хорошо исследованы режимы конвективных течений при надкритичностях до г = Ка/Касг <

© И. Б. Палымский, 2000.

300, Racr = 657.5 [1]. Вплоть до r = 1000 исследовались симметричные решения [2]. Установлено, что при увеличении надкритичности (r > 1) нулевое решение становится неустойчивым и возникают вторичное стационарное решение, затем периодическое и двухчастот-ное. Однако по поводу существования хаотических (сложных) режимов конвективных течений имеются самые противоречивые точки зрения. С одной стороны, в итоговой работе [1] утверждается, что в двумерной задаче о конвекции нет решений, отличных от стационарных, периодических и двухчастотных. С другой стороны, хаотические режимы течений обнаружены в близкой задаче о неустойчивости плоского вращающегося слоя жидкости при подогреве снизу и малом числе Прандтля (Pr = 0.025) [3]. И, наконец, в работе [5] сообщается, что при надкритичности r > 157 и Pr = 20 найдена область хаотических режимов.

Вопрос о существовании хаотических режимов конвективных течений представляется крайне важным ввиду его прямой связи с моделированием турбулентных течений путем прямого численного решения уравнений гидродинамики без использования полуэмпирических соотношений. Для поиска ответа на него, на наш взгляд, наиболее предпочтительно исследовать область умеренных и высоких чисел Рэлея (Ra > 1000Racr) и небольших чисел Прандтля (Pr ~ 1).

Целью данной работы является создание метода численного расчета конвективных течений, работоспособного при Ra > 1000Racr и произвольных числах Прандтля.

Искомые величины w, ^ и Q запишем в виде

N M-1

w(t,x,y) = wfcm(í)pfc cos(akx) sin(nmy),

k=0 m=l N M-l ,

= V V 2 7 2 ,—Pkcos(akx)sin(nmy),

a2k2 + n2m2

k=0 m=l

N -l M -l

Q(t,x,y) = ^2 Qkm(t) sin(akx) sin(nmy),

к=1 т=1

,, Г 0.5 при к = 0, Ж, „ .

где а = п// — волновое число; рк = < 1 ^ ^ к < N 1 ' Следуя общей идеологии

метода расщепления, переход от слоя п к слою п + 1 по времени производится в два этапа. Сначала учтем линейное развитие возмущении без взаимодействия гармоник. Этап 1. 1 11

и = - А и + RaQx, А^ = -и, Цг = А Ц -— ^ (2)

2 2Рг Рг

Для эффективного решения уравнений нелинейного конвективного переноса для завихренности и и температуры Ц половина вязких членов учтена на втором этапе расчета.

После подстановки разложений искомых величин вместо (2) получим систему из двух обыкновенных дифференциальных уравнений для двух неизвестных амплитуд икт и Цкт при к = 0,1, ..., N и т = 1, 2, ..., М - 1:

икт = -2икт + RaаkQkm, Цкт = -Цкт + , S = О^2 + П2т2. (3)

Система обыкновенных дифференциальных уравнений (3) решается аналитически, без применения каких-либо аппроксимаций по времени. Приведенные ниже формулы выведе-

ны программой аналитических вычислений Maple V Release 4:

, = -F3RaPrSQ£m + (F1 + F2)_gm Q , = — F3_gm + (F1 - F2)Qm _km 2S1 ' Qkm 2S1 '

где S1 = ^/S4(1 — Pr)2 + 16SPrRaa2k2;

F1 = (S 2 + S3)S 1; F 2 = S2(Pr — 1)(S2 — S3); F3 = 4ak(S2 — S3);

—т (S2(1 + Pr) + S1) — т (S 2(1 + Pr) — S1)

S2 = exp-4Ps-; S3 = exp-4PS-•

Здесь и далее т — шаг по времени; штрих указывает на значения искомой величины после 1-го этапа расщепления.

На втором этапе учитывается нелинейный конвективный перенос, т. е. взаимодействие гармоник. В этом случае применена конечно-разностная схема переменных направлений, раннее успешно использованная для расчета турбулентных конвективных течений в прямоугольной области при подогреве сбоку [6].

Уравнения переноса для _ и Q решаются в два дробных шага, на каждом из них применяется схема А. А. Самарского для аппроксимации одномерных операторов на верхнем слое по времени и центральными разностями на нижнем. Этап 2.

pr ('y_x — 'x_y ) = 2 А _ Qt + Pr ('уQx — 'xQy) = 2PT

Для первого уравнения системы (4) запишем первый дробный шаг _ч - + А _ 1 _ 1 ,

т/2 + Л _ / У|Н 1\ _гзхх + 2_гШ'

2[ 1 +

Pr

1 (vy + 'Ц , 'у— 'Ц А 1

где

а _ -"-/ry'lryi_ ,r'y \TV\_ I , /г\

A = Pr V-2-Uijx + 2-Ш jx J — pr 'x_ijy* ' (5)

fij — fi-l,j f _ /i+l,j — fi'j

f __ J'j Ji-1,j f

fijx = H1 ' fijx

/ = fi+1,j — fi-1,j / _ fijx* n ГГ1 ' fijxx

H1 , J4X H1 fi+1,j — fi-1,j f _ fi+1,j — 2fij + fi-1,j

/ijx 2H1 ' ijxx h 12 и второй дробный шаг

j — _ij = 1 _ _ 1 n+1

т/2 + A =2_ijxx + / |'x|H2N _jy'

2 1 +

Pr

Л 1 „ -Г 1 Ух + Ы „+1 , Ух - Ы п+1 \

А _ рУу- р ^—2—_я + —2— ' ■ (6)

Здесь Н1 _ 1/У и Н2 _ 1/М — шаги разностной сетки по х и у. Если Н1 и Н2

Лух|Я2 |уу|Н 1 \ . . . .

достаточно малы ——- ^ 1, ——- ^ 1 , то схемы (5) и (6) имеют второй порядок

\ Рг Рг у

аппроксимации по пространству.

Коэффициенты (рх и (ру в разностных уравнениях (5), (6) определялись двумя способами:

'х = 'У = 'У

или

+ 'П

'П+1 + 'П . 'П+1 + '

'х = ---, 'y

2 ' гу 2

Реализация второго способа потребовала введения итерационного процесса. Конечно-разностные уравнения (5) замыкаются соотношениями

—C2j + 4u1j — 3ü0j = 0 при i = 0, —CN-2,j + 4üN-1)j- — 3CNj = 0 при i = N. (7)

Разностные уравнения (5), (6) решаются методом прогонки, а для реализации граничных условий (7) организуется внутренний итерационный процесс:

ü j1 — cj + = 0 при i = 0,

< — üN+-\j + üN-2'j4— "N'j =0 при i = N

Конечно-разностные уравнения для расчета Qn+1 полностью аналогичны уравнениям (5), (6), только на обоих дробных шагах граничные условия для Q и Qn+1 однородные.

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

Предложенный алгоритм является абсолютно устойчивым. В практических расчетах все ограничения на шаги т, H1 и H2 связаны с требуемой точностью вычислений. Численный метод имеет первый порядок аппроксимации по времени и второй — по пространственным переменным.

Чтобы показать, что предлагаемым численным методом можно рассчитывать турбулентные течения, воспользуемся методикой, разработанной для анализа численного метода, которым рассчитывались турбулентные течения в плоском канале и трубе кольцевого сечения [7]. Аналогичный подход использовался также в работе [8]. По этой методике, описанной в данных работах, сравнивались инкрименты линейного развития возмущений в дифференциальной задаче и численном методе. По близости полученных спектральных характеристик можно судить о точности численного метода.

Итак, рассматривается линейный аналог системы (1), в котором нелинейные члены отбрасываются, и решение ищется в виде

ü(t, x, y) = a exp (—At + iakx) sin (nmy), a

'(t, x, y) = — exp (—At + iakx) sin (nmy), S

Q(t, x, y) = b exp (—At + iakx) sin (nmy),

где, как и раньше, S = a2k2+n2m2, a и b — постоянные амплитуды; инкримент A находится из задачи на собственные значения.

Подобные рассмотрения проводятся и для численного метода.

Аналитические выражения для спектральных характеристик можно получить при Pr = 1:

для дифференциальной задачи

Ad = S - ,

для численного метода

. = S k /Ra 1. 1 - (т/4)а + (т2/16)а1а2 Asp =2 - VT - Т 1 + (т/4)а + (т2/16)а1а2'

где

4 . 2 akH 1 4 . 2 nmH2

а1 = --sin -; а2 =--sin -; а = а1 + а2.

H12 2 ' H22 2 '

Можно показать, что

Asp = S - km/Ra + Т2(a6k6 + n6m6) - H^a4k4 - H22n4m4 ,

V S 96 24 24

т2 H12 H22

т 6 6 6 6 H1 4 4 H2 4 4

Asp = Ad + —(a k + п m ) - a k - п m .

Видно, что Asp - Ad = 0(т2) + O(H2) и коэффициенты при т2, H12 и H22 зависят только от a, п, k и m и не зависят от чисел Рэлея и Прандтля. Для дифференциальной задачи при Pr =1:

Ad = ^ SS 2 +

для численного метода

Asp = — ln x,

1 т

где x — наибольшее собственное значение матрицы

( F1 + F2 F3RaPrS \

C C1 2S1

F3 F1 - F2

V -C22s! /

1 - (r/4)a +(t2/16)a1a2 1 - (r/4Pr)a + (т2/16Pr2)a1a2

Здесь C1 = -:——----—--; C 2 =--—--------^-, а величины S,

1 + (т/4)а + (т 2/16)a1a2 1 + (т/4Рг)а + (т2/16Pr2)a1a2'

S1 и F1, F2, F3 определены выше.

Используя эти аналитические формулы для Asp при Pr = 1 с помощью программы Maple V Release 4 и удерживая в степенных разложениях первые шесть членов, получим

Asp = Ad + 0(т 2) + O(H2).

Рассмотрим также спектральные характеристики конечно-разностного численного метода. Для простоты ограничимся изучением аппроксимации по пространству, все производные по x и y заменим соответствующими конечно-разностными выражениями, оставляя

производные по времени дифференциальными. Тогда конечно-разностный метод, использованный в [6] для расчета турбулентного течения в квадратной области при подогреве сбоку, заменяется на систему обыкновенных дифференциальных уравнений

иг = А^и + И.аЦх*, А= -и, Цг = — А^ Q - — <£х*,

Рг Рг

где А/ = /хх + /уу — разностный оператор Лапласа; /х* = ¡^'^д/ ^. Отсюда

Pr + 1 //Pr - 1\ 2 0 Rab2 Ar = -т-^— a — * ——— a2 +

2Pr V V 2Pr J Pra

sin (aкH1)

где b =-—-; a = a1 + a2.

H1

На рис. 1 представлены спектральные кривые, соответствующие первой моде m =1 как функции от ka при Ra = 1000Racr, Pr = 2, a = 1, N = 64, M = 16, т = 4 • 10-5 (сплошная линия — дифференциальная задача, штриховая — спектральный метод, значки — конечно-разностный метод). Видно, что спектральный метод более точен в области волновых чисел, отвечающих нарастающим гармоникам (А < 0).

На рис. 2 изображены среднеквадратичные нормированные значения величин Ad — Asp (сплошная линия) и Ad — Ar (штриховая), вычисленные по нарастающим гармоникам. Эти величины представлены как функции числа Прандтля при 10-2 < 102 (а) и как функции надкритичности r = Ra/Racr (б).

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

Исследуем вопрос о правильности отражения на волновой плоскости границы области неустойчивости. Покажем, что для дифференциальной задачи эта кривая в полярных координатах задается уравнением р = cos0'5 7, где 0 < 7 < п/2; р = p/Ra0'25.

Кривая, ограничивающая область неустойчивости для дифференциальной задачи, и рассчитанные численно аналогичные ей кривые для спектрального и конечно-разностного

численных методов (N = 64, M = 16, т = 4 • 10-5, Ra = 1000Racr, Pr = 1), приведены на рис. 3, где сплошная линия соответствует дифференциальной задаче, штриховая — конечно-разностному численному методу, значки — спектральному численному методу. Видно, что спектральный метод точнее передает положение границы области неустойчивости на волновой плоскости.

Для проверки правильности составления программы произведено сравнение результатов расчета с результатами расчета, полученными другим, "чисто" спектральным методом при N = 64, M = 16, a = 1, Pr = 2, Ra = 2Racr. Среднеквадратичные отклонения полей завихренности, температуры, рассчитываемых средних величин на полученном стационарном решении не превысили 1 %. Максимальный инкримент нарастания возмущений, вычисленный из соотношения

/ u)2dxdy = a exp (2Aspt),

отличался от теоретического значения на 1.4%.

а,-100

10 рг ю2 103 ю4 ю5 г

Рис. 2.

Для проверки порядка аппроксимации при Ra = 1,а = п и Pr = 2 рассчитывалось стационарное решение

ы(х, у) = cos (пх) sin (пу), Q(x, y) = sin (nx) sin (пу), <^(x, y) = ы(х, у)/2п2,

при этом в правую часть уравнений для ы и Q системы (1) вводились соответствующие массовые силы. При различных N, M и т вычислялись среднеквадратичные нормированные отклонения полей ы и Q. Эти тестовые расчеты подтвердили, что предлагаемый численный метод имеет первый порядок аппроксимации по времени и второй — по пространственным переменным.

Несколько слов о выборе числа гармоник N и M. Для создания точной картины течения необходимо (при заданном а) учесть все нарастающие длинноволновые гармоники и достаточное число затихающих коротковолновых гармоник. Как показывает простой анализ спектральных кривых Ad, все гармоники затихают, если k > 29 или m > 6, поэтому выбрано N = 64 (0 < k < 64) и M = 16 (0 < m < 16).

Описанным выше методом рассчитано стохастическое конвективное течение при Ra = 1000Racr, Pr = 2, а = 1. Вычислены следующие средние величины: основная интегральная

Tt-m/Ra0"25

Рис. 3.

характеристика конвективного теплообмена — число Нуссельта

а Г п/а N-1 М-1 т

^ =- (у (I, х, 0)<х - 1, Nu — V V (1 - (-1)к)(1т- - 1; П ^ 0 __1 ^

интеграл по области от квадрата завихренности

Г1 Г1 I N М-1

Е = ы2<1х<1у ^ Рк^

^ ^ 4 к=0 т=1

Г 0.5 при г = 0, „

где рк ^<1 . > о ,а также (8Г — среднеквадратичное отклонение температу-

ры от ее равновесного распределения; Е8Т — среднеквадратичное значение функции тока; (12 — среднее значение амплитуды во времени.

По результатам расчетов создан видеофильм о развитии во времени функции тока р и температуры Сопоставление видеофильма и графика изменения Nu во времени (рис. 4) показало, что его локальные максимумы на графике связаны с рождением и разрушением вихревых структур. Знаком □ показано рождение вихревой структуры, состоящей из четырех вихрей, а знаком х — ее разрушение. Начальные перестройки течения, связанные с выделением и преимущественным развитием наиболее быстрорастущей в линейном приближении гармоники (£ — 0.06) и последующим ее разрушением нелинейностью (£ — 0.1), показаны знаком •.

Изменение во времени амплитуды (12 представлено на рис. 5, эта гармоника вносит существенный вклад в распределение температуры.

Из тестовых расчетов видно, что все средние характеристики течения практически не изменяются при варьировании начальных данных, изменении шага по времени т, увеличении числа гармоник по х, а увеличение числа гармоник по у приводит к изменению среднего профиля температуры в приграничных точках.

Таким образом, можно сделать следующие выводы:

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

2. Предлагаемый численный метод абсолютно устойчив и имеет первый порядок аппроксимации по времени и второй — по пространству.

3. Показано, что спектральные характеристики численного метода близки к спектральным характеристикам дифференциальной задачи. Спектральные характеристики конечно-разностного метода [6] значительно хуже аппроксимируют спектральные характеристики дифференциальной задачи.

Рис. 4.

Рис. 5.

4. При Ra = 1000Racr, Pr = 2 и a = 1 рассчитано стохастическое конвективное течение. Средние характеристики течения и средний профиль температуры слабо зависят от начальных условий, числа гармоник по пространственным переменным и величины шага по времени.

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

[1] БАБЕНКО К. И., РАХМАНОВ А. И. Численное исследование двумерной конвекции. М., 1988 (Препр. / АН СССР. ИПМ; №118).

[2] Curry j. h., herring j. r., Loncaric j., Orszag s. a. Order and disorder in two-and three-dimensional Benard convection //J. Fluid Mech. 1984. V. 147. P. 1-38.

[3] ГЕРЦЕНШТЕйН С. Я., РОДИЧЕВ Е. Б. О моделях перехода к турбулентности при конвективной неустойчивости // Моделирование в механике. 1989. Т. 3(20), №4. С. 5965.

[4] ПЕТРОВСКАЯ Н. В. О применении метода Галеркина к исследованию переходов в задаче рэлеевской конвекции // Изв. АН СССР. Серия МЖГ. 1984. №2. С. 22-27.

[5] РОДИЧЕВ Е. Б., РОДИЧЕВА О. В. О двумерной турбулентности в задаче Рэлея — Бенара // Докл. АН СССР. 1998. Т. 359, №4. С. 486-489.

[6] Пасконов В.М., Полежаев В. И., ЧУДОВ Л. А. Численное моделирование процессов тепло- и массообмена. М.: Наука, 1984. 288 с.

[7] РОЖДЕСТВЕНСКИЙ Б. Л., СТОЙНОВ М. И. Алгоритмы интегрирования уравнений Навье — Стокса, имеющие аналоги законам сохранения массы, импульса и энергии. М., 1987 (Препр. / АН СССР. ИПМ; №119).

[8] НИКИТИН Н. В. Спектрально-конечно-разностный метод расчета турбулентных течений несжимаемых жидкостей в трубах и каналах // Журн. вычисл. математики и мат. физики. 1994. Т. 34, №6. С. 909-925.

Поступила в редакцию 8 февраля 2000 г.