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

Том ХЫ

УЧЕНЫЕ ЗАПИСКИ ЦАГИ 2010

№ 2

УДК 534.83 + 532.511 536.25

РАСЧЕТЫ ТЕРМОАКУСТИЧЕСКОЙ КОНВЕКЦИИ НА МНОГОПРОЦЕССОРНОЙ ЭВМ

А. А. ГОРБУНОВ, С. А. НИКИТИН, В. И. ПОЛЕЖАЕВ

Численно исследуется одномерная и двумерная термоакустическая конвекция в замкнутом объеме, заполненном углекислым газом при нормальных условиях (совершенный газ) или вблизи критической точки (газ Ван-дер-Ваальса), в диапазоне чисел Рейнольдса 103 — 105 на параллельной ЭВМ Межведомственного Суперкомпьютерного Центра РАН. Достигнутые ускорения расчетов по сравнению с однопроцессорной ЭВМ равны 60 для одномерной задачи и 300 для двумерной.

Исследованы поля течения и температуры в зависимости от параметров задачи, показана необходимость использования подробной разностной сетки в обоих направлениях для разрешения тонких пограничных слоев и фронтов (104 узлов сетки для числа Рейнольдса, равного 105).

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

Термоакустические (ТА) течения относятся к классу естественных течений негравитационного типа, возникающих в сжимаемой среде при подводе тепла. Иногда к этому классу течений относят и течения, возникающие при взаимодействии тепловой гравитационной конвекции с заданными акустическими возмущениями. Несмотря на то, что ТА течения известны давно (по-видимому, первые упоминания об этом механизме течения содержатся в [1]), они до настоящего времени мало изучены. Теплообмен, обусловленный ТА течениями, представляет особый интерес в невесомости и вблизи термодинамической критической точки, что привлекло внимание к ним в эпоху развития космических исследований [2, 3].

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

В 70 — 80-х годах были существенно продвинуты приближенные теоретические модели и начато применение ТА явлений в криогенной технике (см., например, [6, 7]). Модель осреднен-ных течений, вызываемых взаимодействием заданных акустических возмущений и тепловой конвекции, рассмотрена в [8].

Однако работы, относящиеся к численному моделированию ТА, выполненные в последние годы, немногочисленны [9, 10]. Интерес к численному моделированию ТА течений возник в последнее время в связи с расчетами трехмерных конвективных течений сжимаемых небуссинесков-ских сред, в том числе вблизи критического состояния (влияние на формирование небуссинесков-ских трехмерных структур, вклад в теплообмен и др.) [11]. В первых одномерных расчетах ТА течений совершенного газа [4, 5], а также околокритических жидкостей [12] использованы

недостаточно подробные для таких течений сетки, что было обусловлено ограниченным быстродействием ЭВМ в те годы. Высокие требования к сеточному разрешению предполагают использование высокопроизводительных многопроцессорных ЭВМ, и в этом отношении в данной работе продолжен подход, предпринятый в [13]. Ввиду значительных требований к сеточному разрешению этих течений использована параллельная ЭВМ МВС-100К Межведомственного супер-компьютерного центра РАН.

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

Эр + Э(ри) + Э(ру) =о

дt Эх Эу

du _ 1 dp + 1

dt ру dx p Re

4 d2u d2u 1 d2v

3 dx dy 3 dxdy

dv _ 1 dp + 1

dt ру Эу p Re

4 d2v d 2v 1 d 2u

3 dy2 dx2 3 ЭхЭу

d-L__(T_ 1)T (

dt

Y (Y-1)

du dv

- + -

\

p ^ dT )p ^ dx dy ) p Re Pr

Y ( d 2T d 2T"

dx2 + dy2

+

+

p Re

(du + dv dy dx

+

(du dv dx dy

+

dy

+

p T 9

p =pT — совершенный газ; p _—-—-----------p2 — газ Ван-дер-Ваальса (VDW).

1 -p/3 8

Здесь u, v — компоненты вектора скорости; p — плотность; p — давление; Т — температура; d/dt = d/dt + ud/dx + vd/dy.

При приведении системы к безразмерному виду приняты масштабы: длины — L, температуры — Тс (критическая температура), давления — RpcTc (pc — критическая плотность),

скорости — (yRTc)112, времени — L/(yRT;)1/2. В уравнениях появляются безразмерные параметры:

1/2

Y = cp/cv — для совершенного газа; Pr = |icpA, — число Прандтля, Re = pc(yRTc) L/ц — число Рейнольдса. Так как в дальнейшем будут рассматриваться и сравниваться газ VDW и совершенный газ, в приведенном выше безразмерном виде для совершенного газа Tc, pc заменяются на То, po — начальные значения температуры и плотности.

Начальные условия: неподвижный, однородный по температуре, плотности и давлению газ находится при температуре То (То = 1.01 Тс для газа VDW). В начальный момент времени температура левой стенки (x = 0) повышается до Т1 (Т1 = 1.1 Т0 для совершенного газа и Т1 = 1.1 Тс для газа VDW), что приводит к возникновению в области ТА конвекции, затухающей со временем. Граничные условия для компонент скорости на всех границах — условия прилипания, температура левой стенки равна Т1, температура правой стенки (x = 1) равна Т2 = То, остальные границы — теплоизолированы.

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

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

По разработанной методике и программам на параллельной ЭВМ МВС-100К Межведомственного суперкомпьютерного центра (www.jscc.ru) были проведены расчеты задач 1Б и 2Б ТА конвекции. Исследовались зависимости времени расчета задач от числа процессоров. На рис. 1, а представлена зависимость ускорения вычислений от количества процессоров при расчетах одномерной задачи на сетке 55 440 узлов, а на рис. 1, б — аналогичная зависимость для двумерной задачи на сетке 1680 • 1680 узлов. Видно, что существует число процессоров, при котором ускорение максимально (для Ш задачи — максимальное ускорение 65 на 800 процессорах, а для 2Б задачи — максимальное ускорение 300 на 256 процессорах). Аналогичные результаты были ранее получены на ЭВМ МВС-1000М и теоретически объяснены в [13]. На рис. 1, б обращает на себя внимание то, что в диапазоне 100 300 процессоров ускорение выше теоретической

прямой, это, по-видимому, связано с влиянием кэшей процессоров.

Ускорение Ускорение

а) б)

Рис. 1. Ускорение вычислений для 1D (а) и 2D (б) задач термоакустики в зависимости от числа процессоров (ядер)

для ЭВМ МВС-100К:

1 — теоретическая кривая; 2 — реальное ускорение

Результаты расчетов. Выполнено моделирование затухания одномерной и двумерной

ТА конвекции в совершенном газе и газе VDW (двуокись углерода СО2, у = 1.3, Pr = 0.726) для

3 5

значений числа Рейнольдса: Re = 10 — 10 и HIL = 0.25 для 2D задач.

Рассмотрим результаты для одномерной задачи (Re = 10 ). Профили температуры на начальном этапе (t < 1, рис. 2) показывают существенный вклад ТА течения «на хвосте» температурной волны, что обусловлено эффектами вязкой диссипации и адиабатического сжатия. При этом отличия указанных эффектов (проявляющихся в характерных поднятиях, а затем выравниваниях «на хвосте» температурной волны) в совершенном газе и газе VDW очень существенны. На промежуточном этапе (t < 10, рис. 3) эффект адиабатического сжатия («piston» эффект [12]) приводит к полному выравниванию упомянутых температурных неоднородностей в центре слоя, наиболее явно проявляющихся в газе VDW. Отметим также для более отчетливого понимания

T - T

результатов, что в газе VDW удаление от критической точки составляло є =--------------— = 0.01, причем

Tc

на боковой стенке є = 0.1 (рис. 2, 3).

На рис. 4, а, б показаны зависимости теплового потока (т. е. числа Нуссельта) на правой границе (х = 1) от времени на промежуточном интервале прогрева совершенного газа и газа VDW в одномерном и двумерном случаях при Re = 103. Видно резкое увеличение и резкий спад тепло-

Рис. 2. Профили температуры на начальном этапе (Г = 0.5) развития ТА течения в совершенном (1) и УБШ (2) газах

Рис. 3. Профили температуры на промежуточном этапе (Г = 10) развития ТА течения в совершенном (1) и УБШ (2) газах

Ыи

Ыи

Рис. 4. Зависимость № (х = 1) от времени в совершенном газе (а) (1Б задача — 1, 2Б задача — 2) и газе УБШ (б)

(1Б задача)

вого потока, когда ТА волна периодически достигает правой границы. В двумерном случае (кривая 2 на рис. 4, а) заметно существенное ослабление конвекции и уменьшение теплообмена в связи с влиянием трения на боковых стенках. Важно отметить, что перенос тепла к границе х = 1 на порядок больше в газе VDW (рис. 4, б), чем в совершенном газе, а скорость волны в газе VDW меньше, чем в совершенном газе. Аналогично и влияние уравнения состояния и размерности задач на затухание ТА течения (рис. 5). Здесь следует отметить, что максимальные значения скорости ТА течений, имеющие место в начальный период времени, более чем на порядок меньше скорости звука. Видно также, что из-за трения на боковых стенках затухание намного быстрее в двумерном случае (Г = 40), чем в одномерном (Г = 100), причем в газе VDW затухание медленнее (Г = 180).

Значения безразмерной (через скорость звука) скорости, равное 10 4, принято как условие, при котором ТА конвекция предполагается затухшей. Это связано с тем, что обычное значение числа Рейнольдса (определенное через вязкость и характерный размер) равно 0.1, т. е. течение оказывается в области так называемых вязких медленных течений, где роль конвективных членов пренебрежимо мала.

Рис. 5. Затухание ТА течений в совершенном газе (1Б задача — 1, 2Б задача — 2) и газе УБШ (1Б задача — 3)

Рис. 6. Изолинии модуля вектора скорости для 2D ТА течения совершенного газа при Яе = 103 и Г = 0.8 (а, ишах = 8.5 • 10-3), Г = 1.5 (б, ишах = 4.7 • 10-3), Г = 10 (в, мшах = 2.9 • 10~4)

Величина первого импульса этого теплового потока была выбрана для оценки требуемой разностной сетки (изменения этого значения менее чем на 1% при увеличении сетки вдвое, означали, что разрешение достаточно для последующих расчетов). Ниже в таблице приведены результаты этих методических расчетов, на основе которых были выбраны сетки по оси х (средний столбец). В 2Б случаях использовалась квадратная сетка (Их = ку = к) с числом узлов по х, как в 1Б задачах. Шаг по времени был равен к^[к.

Re (размер сетки по x)/Nu

103 (100)/5.12 (200)/5.43 (400)/5.45

104 (500)/16.1 (1000)/17.65 (2000)/17.75

105 3000/53.4 (6000)/57.81 (12 000)/58.02

Интересные результаты относятся к полю скорости ТА течений на разных режимах во времени. Здесь (рис. 6) приведены изолинии поля модуля скорости для трех характерных моментов времени: а — до столкновения фронта волны с противоположной стенкой (t = 0.8); б — возвратное течение после первого столкновения с противоположной стенкой (t = 1.5); в — структура ТА течения после многократного отражения от противоположной стенки (t = 10). Последнее имеет очень сложный характер, в котором разделяются пристеночные пограничные слои, течения в ядре и противоположное ему возвратное течение между ядром и пограничным слоем. Отметим также, что для разрешения структуры течения в пристеночном пограничном слое применяемые в работах [8, 9] разностные сетки в поперечном сечении (50 узлов) недостаточны. Вместе с тем следует заметить, что диапазон определяющего критерия подобия — числа Рейнольдса, отнесенного к скорости звука, на практике может на несколько порядков превышать максимальное значение, достигнутое в данной работе в двумерном случае, поэтому потребности в продолжении методических исследований для численного решения этого класса задач остаются.

Выводы. Применение высокопроизводительных, в том числе многопроцессорных ЭВМ, в которых достигнуто ускорение вычислений от 60 (1D) до 300 (2D), позволяют получить качественно новые результаты по структуре термоакустических течений и переноса ими тепла, выявить характерные режимы течений во времени и при изменении основного критерия подобия числа Re совершенного газа и газа VDW вблизи критической точки CO2 (є = 0.01), а также характерные эффекты теплообмена и затухания при наличии вязкой диссипации, адиабатического сжатия в одномерных и двумерных течениях. Отметим, что большая часть этих эффектов не была известна и не могла быть реализована при применявшемся ранее пространственном разрешении ТА течений.

Работа выполнена при финансовой поддержке Ведущей научной школы № 2496.2008.8, Российского фонда фундаментальных исследований (гранты 06-01-00281 и 09-01-00230), а также программы Президиума РАН № 17.

ЛИТЕРАТУРА

1. Релей. Теория звука. — М. — Л.: ОГИЗ, 1944, 476 с.

2. Повицкий А. С., Любин Л. Я. Основы динамики и тепломассообмена в жидкостях и газах при невесомости. — М.: Машиностроение, 1972, 252 с.

3. Полежаев В. И. Конвекция и процессы, тепло- и массообмена в условиях космического полета // Изв. РАН. МЖГ. 2006. № 5, с. 67 — 88.

4. Полежаев В. И. Численное решение системы одномерных нестационарных уравнений Навье — Стокса для сжимаемого газа // Изв. АН СССР. МЖГ. 1966. № 6, с. 34 — 43.

5. Васин В. Г., Полежаев В. И. Разностные схемы для уравнений Навье — Стокса сжимаемого газа и расчет термоакустических волн. — М.: ИПМ АН СССР. Препринт № 124,

1977, 47 с.

6. R o tt N. Thermoacoustics / In: Advances in Applied Mechanics 1980. V. 20, p. 135 — 175.

7. Wheatly J., Cox A. Natural engines // Physics Today. 1985, N 5, p. 50 — 57.

8. Любимов Д. В. О тепловой конвекции в акустическом поле // Изв. РАН. МЖГ.

2000. № 2, c. 28 — 36.

9. F a r o u k B., O r a n E. S., F u s e g i T. Numerical study of thermoacoustic waves in enclosure // Physics of Fluids. 2000. V. 12, N 5, p. 1052 — 1061.

10. Aktas M. K., Farouk B., Narayan P., Weatley A. A numerical study of the generation and propagation of thermoacoustic waves in water // Physics of Fluids. 2004. V. 16,

N 10, p. 3786 — 3794.

11. Горбунов А. А., Полежаев В. И. Метод возмущений и численное моделирование конвекции для задачи Рэлея в жидкости с произвольным уравнением состояния. —

М.: ИПМех РАН. Препринт № 897, 2008, 46 c.

12. Zappoli B.Durand-Daubin A. Heat and mass transport in a near supercritical fluid // Physics of Fluids. 1994. V. 6, N 5, p. 1929 — 1936.

13. Горбунов А. А., Никитин С. А., Полежаев В. И. Численное моделирование конвекции Релея — Бенара в сжимаемом газе на параллельных вычислительных системах //

Вестник Нижегородского ун-та. Сер. Матем. моделир. и опт. управл. 2005, вып. 1 (28), с. 75 — 82.

Рукопись поступила 13/III2009 г.