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

ГАЗОДИНАМИКА И ГОРЕНИЕ

УДК 532.510.5

М.М. Горохов, И.Г. Русяк, А.А. Бас, А.В. Корепанов, А.В. Милюков

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

ПРОСТРАНСТВЕННЫХ ТУРБУЛЕНТНЫХ ТЕЧЕНИЙ

Ижевский государственный технический университет

Введение

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

1. Основные возможности и структура программного комплекса

Основными возможностями программного комплекса являются:

- определение аэродинамических характеристик поверхностей в турбулентных потоках;

- моделирование эрозионного горения твердого гранулярного топлива;

- аэрация городских территорий;

- расчет загрязнения от автотранспорта в условиях городской застройки;

- расчет параметров конвективного теплообмена зданий;

- моделирование распределения загрязнения в воздушной среде.

В состав взаимосвязанных составных частей програм м но-вычис л отельного комплекса входят: блок численного моделирования, модуль распараллеливания вычислений, база данных (БД), блок структурных запросов к базе данных, топологическая основа региона, геоинформационная система, блок визуализации данных (рис. І).

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

БД

О

Блок структурных запросов

Г еоинформационная система О Блок численного моделирования О Блок визуализации

Топологическая основа региона

Модуль

распараллеливания

вычислений

Рис. 1. Структурная схема программного комплекса

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

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

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

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

Ниже представлены формы поверхностей, около которых могут быть определены параметры пространственных течений (рис. 2),

2, Математические модели

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

■^ + ¥р¥¥ = -¥| р + УУ| +

ы к ъ )

< + 2У(цёеГУ), (1)

— + Ур¥ = 0, ог

где р - плотность газа; V - вектор скорости; I - время; р—давление; р = ри + ц, - эффективная вязкость, определяемая как сумма молекулярной и турбулентной составляющих.

Для учета распределения температуры в потоке к уравнениям (1) следует добавить уравнение энергии:

+ ¥рУ Т = V (л¥Г). (2)

Здесь Т - температура; X - Хт + %, - коэффициент теплопроводности представлен в виде молекулярной и турбулентной составляющих.

'-т.

*

!::Щ|

ЯД/'. '..--А \ \ 1 Г I /

-Л_Л „.Л-Г !

. *, I ■ *. I Г./

* Iгч

..Л I

\

Рис. 2. Типы поверхностей, около которых исследуются параметры течений

Система уравнений (1) — (2) замыкается уравнением состояния. Для связи коэффициента молекулярной вязкости \іт и абсолютной температуры Т используется зависимость Саттерлзнда [1].

Область численного интегрирования ограничена твердой поверхностью, входной, выходной и внешней границами, расположенными во внешнем течении. На твердой поверхности задаются температура, условия прилипания и непротекания, на входной и внешней границах выставляются параметры набегающего потока, на выходной - давление и «мягкие» граничные условия для прочих переменных.

Для описания турбулентных режимов течения при решении задач обтекания инертных тел и горящих гранул используется двухпараметрическая модель турбулентности [2].

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

£Р±. + урУТ = V(%УТ) + ~~/г (С,, Т) +

+|1./;(с2,г), о)

-^Ь-УрУс, = ії-ус, 1- /;(с.,г),

а \ 8с ') 1 * ;

‘ УрУС2 =у\^Сгу /{С,., Т)- (4)

-/2(С,г).

Здесь & - число Шмидта; с -теплоемкость продуктов горения при постоянном давлении; Qj - тепловой эффект реакции первой (і = 1) и второй (і = 2) стадий химического превращения; С, /'.(С., Г.) - концентрация продуктов и скорость химической реакции соответствующей (г = 1, 2) стадии [3].

На поверхности горения ставятся условия теплового баланса конденсированной и газовой фаз. Скорость горения твердого топлива вычисляется по формуле Мержанова-Дубовицкого [4].

Математическая модель распространения загрязняющих веществ в приземном слое атмосферы, а также модель расчета аэрационного режима зон с городской застройкой опираются на уравнения (1) - (2) и уравнение для переноса I концентрации примеси [5]:

^ + Щ=Ч(к,\\), (5)

где qt - концентрация i примеси: kd - коэффициент обмена, определяемый согласно [6]. Проекция вектора скорости на вертикальную ось координат определяется как Voc = D{ur\x +Щу +(w- ww)ti„), vM = 2ptg rf j9\im - скорость осаждения;pt,rf -плотность, радиус частиц примеси; g - ускорение свободного падения.

Для получения логарифмического распределения скорости ветра в приземном слое атмосферы используется однопараметрическая модель турбулентности [7].

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

3. Методика численного расчета

Для построения конечно-разносты»-;у сеток область интегрирования в направлен ора скорости делилась секущими плоское 1а каждой

плоскости строилась адаптированная к геометрии границ конечно-разностная сетка на основании комплексного метода граничных элементов [8]. Впервые такой подход к построению конечно-разностных сеток был применен в работе [9]. Построенные на основании данного метода линии сетки по сути являются линиями тока и эквипотенциальными линиями безвихревого течения. Это свойство является достаточно значимым, поскольку до появления в потоке отрывных зон сетка адаптирована не только к геометрии твердой поверхности, но и к течению потока газа, что позволяет свести к минимуму влияние «схемной» вязкости. Линии конечно-разностных сеток сгущались к твердой поверхности таким образом, чтобы значение сеточного числа Рейнольдса, вычисленного по параметрам набегающего потока в ближайшем от поверхности узле, соответствовало единице (рис. 3),

Для численного решения системы уравнений использовался метод SIMPLE [101, реализованный на параллельных вычислительных системах. Для применения параллельных алгоритмов и синхронизации параллельных вычислений область интегрирования делилась на равные части, локальные разделы, содержащие одинаковое количество узлов. Количество разделов соответствовало количеству используемых процессоров. Разделы состоят из совокупностей секущих плоскостей, располагающихся параллельно

Рис, 3. Конечно-разностные сетки

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

Для организации кластера использована архитектура «клиент-сервер». Задачей «клиента» является решение уравнений в отведенном для «клиента» разделе области интегрирования на текущем итерационном слое. Задачами «сервера» являются синхронизация процессов «клиентов», сбор и визуализация данных, принятие решения о переходе к следующему итерационному слою или о завершении вычислительного процесса, представленного в виде блок-схемы (рис. 4).

Для организации обмена данными между клиентом и сервером были задействованы компоненты Delphi 5.0 и Builder 5.0 [11]: TServerSocket, TClientSocket для пересылки текстовых сообщений; TNMStrm, TNMStrmServ - для пересылки массивов.

При проведении вычислений сервером являлась вычислительная система Pentium IV1400 МГц

НАМ 512 МВ, а в качестве клиентов использовалось 10 вычислительных систем Сеіегоп 1000 МГц ВАМ 256 МВ объединенные в ЛВС 100 МЬ.

Для решения систем линейных алгебраических уравнений, являющихся дискретными аналогами уравнений количества движения и переноса концентраций, применялся метод продольно-поперечной прогонки, а при решении дискретного аналога уравнения для поправок к давлению [10] - модифицированный сильно-неявный алгоритм [12] в совокупности с методом сопряженных градиентов [13].

Сходимость метода численного решения была установлена путем измельчения разностной сетки и варьирования размеров области численного интегрирования. В качестве тестовой рассматривалась задача обтекания сферы в диапазоне изменения значения числа Рейнольдса от 10 до 106. Из расчетных значений коэффициента полного сопротивления сферы от числа Рейнольдса и стандартной кривой сопротивления [14, 15] видно (рис . 5), что в диапазоне изменения значения числа Рейнольдса от 10 до 2><105 расчетные и экспериментальные значения удовлетворительно согласуются с аппроксимационными зависимостями. При числе Рейнольдса ~3 х 105, соглас-

СЕРВЕР

Прием потоков от клиентов

Принятие решения о продолжении процесса

£

КЛИЕНТ

Прием потока и разрежения на начало процесса

Решение уравнений гидромеханики д^"

Отправка потока на сервер

СЕРВЕР

Расформирование кластера

Рис. 4. Блок-схема распараллеливания вычислительного алгоритма

Таблица 1

0.01

10 102 103 104 105 106 Р!е

Рис. 5. Зависимость значения коэффициента сопротивления шара от числа Рейнольдса:

1 - экспериментальные данные [14-16]; 2 - расчетные значения

Же Сх

расчетное экспериментальное (1)

251 300 0.297 0.313

298 500 424 500 0,172 0.151

0,209 0.143

но экспериментальным данным, ламинарный пограничный слой на сфере переходит в турбулентный. В этом случае расчеты показали, что происходит уменьшение значения коэффициента полного сопротивления. Этот результат соответствует экспериментальным данным [1, 16]. Представлено сравнение расчетных и экспериментальных значений Сх в диапазоне чисел Рейнольдса от 2.51><105до 4.24хЮ5 (табл. 1).

Показано также (рис. 6) изменение поля течения около сферы при различных значениях числа Рейнольдса.

Рис. 6. Влияние числа Рейнольдса на картину течения около сферы: а - Яе = 1.71 х105; б~ Яе = 4.35x10*

Как видно из рисунков, размер циркуляционной зоны за сферой увеличивается с ростом числа Рейнольдса. При достижении критического значения числа Рейнольдса -3*105 происходит смещение вниз по течению циркуляционной зоны и уменьшение ее размера, что подтверждают экспериментальные данные [1, 16]. Наибольшие значения турбулентной вязкости соответствуют циркуляционной зоне (рис. 7).

Рис, 7. Распределение турбулентной вязкости, отнесенной к молекулярной вязкости при ,% = 4.35-105

Для тестирования параллельного алгоритма была решена задача вязкого пространственного обтекания сферы при числе Рейнольдса 1О4, Сходимость вычислительного алгоритма контролировалась установлением значения величины полного сопротивления-сферы Сх и поправки к давлению ар. При распараллеливании алгоритма возрастает количество итераций М, требуемых для установления численного решения. Приводим данные зависимости количества N и значений Сх от числа используемых процессоров кластера q, а также астрономическое время проведения расчетов (табл. 2). Во всех случаях величина Сх удовлетворительно согласуется с экспериментальными данными [1].

Таблица 2

е N Сх 1<р с

1 304 0.473 13 020

2 340 0.476 8 137

4 443 0.484 3 921

6 530 0.491 2 726

8 591 0.496 2 246

10 612 0.499 2 081

Ускорение параллельного алгоритма оценивалось согласно [17]: 5^ = где ^l - астрономическое время выполнения алгоритма на одном процессоре, / -астрономическое время выполнения алгоритма на д процессорах.

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

Показаны картины течения около различных поверхностей (рис. 8-14). В частности, представлены (рис. 8) мгновенные линии тока для случая течения около сферы при Не = 104. Отметим, что в расчетах в качестве граничного условия задавалась «закрутка» потока. Угловая скорость «закрутки» составляла -40 % скорости набегающего потока. Показано обтекание инертной сферы (рис. 8, а), обтекание сферы с вдувом (рис. 8, б-г). Вдув с поверхности приводит к увеличению области влияния обтекаемого тела (рис. 8, а, б). Структура течения внутри отрывной зоны (рис. 8, в), а также линии тока в

а

В

б

в

г

Рис, 8. Мгновенные линии тока при обтекании сферы

отрывной зоне и линии, «выходящие» с поверхности тела, показывают (рис. 8, г), что вдув с поверхности приводит к смещению отрывной зоны вниз по направлению течения.

а

Мгновенные линии тока для случая обтекания эллипсоида вращения представлены на следующих рисунках: обтекание инертного тела (рис. 9, а), обтекание тела с вдувом (рис. 9, б-г),

б

Рис, 9. Мгновенные линии тока при обтекании эллипсоида вращения

Расчеты эрозионного горения гранулы топлива I/,. =2 м/с (рис, 10, а, о), и,. = 400 м/с (рис. 10, в, г).

проведены для температуры внешнего течения Тж = Скорость £/х = 400 м/с соответствует турбулентно-

= 2 300 К, давления рх =10 МРа и скоростей потока му режиму течения. Начальная угловая скорость в

а б

Рис, 10. Мгновенные линии тока при обтекании криволинейной поверхности

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

а

ния обтекаемого тела, а неравномерность вдува и поверхности гранулы приводит к «закрутке» и потере симметричности течения.

Структура течения около последовательно расположенных горящих гранул зависит от скорости внешнего течения; 1]т = 2 м/с (рис. 11, а, б), 11х- 200 м/с (рис. 11, в, г). Продукты горения оттесняют линии тока

б

Г

Рис. 11. I едение около последовательно расположенных поверхностей

от поверхности обтекаемых тел, а «закрутка» течения вызвана влиянием расположенной ниже по потоку7 гранулы, Изотермы около горящих сфер (рис. 12) указывают на то, что даже при наличии угловой скорости потока картина течения остается симметричной.

Как было указано выше, программный комплекс способен решать задачи аэрации жилой застройки и экологического мониторинга города. Процесс деформации воздушного потока в зоне жилой застройки показывает (рис. 13, а), что режим течения между

зданиями при выбранном направления ветра существенно разнороден. Здесь имеет место как сквозная продувка района застройки с возрастанием скорости потока, так и процесс вихреобразования, сопряженный с вертикальным подъемом воздушных масс (рис. 13, б). Расчет распространения концентраций от городских магистралей произведен с учетом заданного направления и силы ветра (рис. 14, а), а расчет месячного рассеяния примеси - с учетом розы ветров (рис. 14, б).

Рис. 12. Изотермы около горящих сфер: а-их = 2 м/с; б - = 200 м/с

Рис. 13. Течение в зоне застройки

щ

ШШШГ,

ш

вШШвШшЖ

Рис. 14. Распространение загрязнений от городских автомагистралей

Заключение

Предлагаемый в работе подход к моделированию турбулентных течений, положенный в основу программно-вычислительного комплекса, может быть

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

Литература

Лойцанский Л.Г. Механика жидкости и газа. М., 1987,

Курбацкий А.Ф. Моделирование турбулентных течений (обзор) // Изе. СО АН СССР. 1989, Вып. 6,

Горохов М.М., Русяк И.Г, Моделирование эрозионного горения гранулированного топлива // ФГВ. 2001. Т. 87. № 3. Мержанов А.Г., Дубовицкий Ф.И. К теории стационарного горения пороха // ДАН СССР. 1959. Т, 129. N8 1,

Берлянд М.Е, Современные проблемы атмосферной диффузии и загрязнения атмосферы. Л., 1975,

Атмосферная диффузия и загрязнение воздуха / Под ред, А.С. Монина. М,, 1972,

Секундов А.Н, Применение дифференциального уравнения для турбулентной вязкости и анализ плоских неавтомодельных течений // Изв. АН СССР. МЖГ. 1971. № 5.

Громадка Т., Лей Ч. Комплексный метод граничных элементов в инженерных задачах. М., 1990.

Тененев В.А. и др. Численное исследование горения частиц алюминия в двухфазном потоке // РАН. Математическое моделирование, 1997. Т. 9. № 5.

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

11, Архангельский А. Программирование в C++ Builder 6. 2003.

12, Андерсон Д. и др. Вычислительная гидромеханика и теплообмен, Т. 2. М., 1990.

13, Ильин В.П., Лившиц И.З, Симметризованный метод Стоуна //' ЖВМиМФ. 1994, Т. 34, № 11.

14, Двухфазные моно- и лолидисперсные течения газа с частицами / Под ред. Л.Е. Стернина. М., 1980.

15, Шрайбер А,А, и др. Гидромеханика двухкомпонентных потоков с твердым лолидисперсным веществом. Киев, 1980.

16, Ван-Дайк М. Альбом течений жидкости и газа. М., 1986,

17, Ортега Д. Введение в параллельные и векторные методы решения линейных систем, М,, 1991,