авторефераты диссертаций БЕСПЛАТНАЯ  БИБЛИОТЕКА

АВТОРЕФЕРАТЫ КАНДИДАТСКИХ, ДОКТОРСКИХ ДИССЕРТАЦИЙ

<< ГЛАВНАЯ
АГРОИНЖЕНЕРИЯ
АСТРОНОМИЯ
БЕЗОПАСНОСТЬ
БИОЛОГИЯ
ЗЕМЛЯ
ИНФОРМАТИКА
ИСКУССТВОВЕДЕНИЕ
ИСТОРИЯ
КУЛЬТУРОЛОГИЯ
МАШИНОСТРОЕНИЕ
МЕДИЦИНА
МЕТАЛЛУРГИЯ
МЕХАНИКА
ПЕДАГОГИКА
ПОЛИТИКА
ПРИБОРОСТРОЕНИЕ
ПРОДОВОЛЬСТВИЕ
ПСИХОЛОГИЯ
РАДИОТЕХНИКА
СЕЛЬСКОЕ ХОЗЯЙСТВО
СОЦИОЛОГИЯ
СТРОИТЕЛЬСТВО
ТЕХНИЧЕСКИЕ НАУКИ
ТРАНСПОРТ
ФАРМАЦЕВТИКА
ФИЗИКА
ФИЗИОЛОГИЯ
ФИЛОЛОГИЯ
ФИЛОСОФИЯ
ХИМИЯ
ЭКОНОМИКА
ЭЛЕКТРОТЕХНИКА
ЭНЕРГЕТИКА
ЮРИСПРУДЕНЦИЯ
ЯЗЫКОЗНАНИЕ
РАЗНОЕ
КОНТАКТЫ


Pages:   || 2 |
-- [ Страница 1 ] --

ФЕДЕРАЛЬНОЕ АГЕНТСТВО ПО ОБРАЗОВАНИЮ

_

МОСКОВСКИЙ АВИАЦИОННЫЙ ИНСТИТУТ

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

А. А. ЮН, Б. А.

КРЫЛОВ

РАСЧЕТ И МОДЕЛИРОВАНИЕ

ТУРБУЛЕНТНЫХ ТЕЧЕНИЙ

С ТЕПЛООБМЕНОМ, СМЕШЕНИЕМ,

ХИМИЧЕСКИМИ РЕАКЦИЯМИ И

ДВУХФАЗНЫХ ТЕЧЕНИЙ В ПРОГРАММНОМ

КОМПЛЕКСЕ FASTEST-3D

Допущено Министерством образования и науки Российской Федерации в качестве учебного пособия для студентов высших учебных заведений, обучающихся по специальности «Авиационные двигатели и энергетические установки» направления подготовки «Двигатели летательных аппаратов»

Москва Издательство МАИ 2007 ББК 27.5.4 Ю49 Ю49 А.А. Юн, Б.А. Крылов. Расчет и моделирование турбулентных течений с теплообменом, смешением, химическими реакциями и двухфазных течений в программном комплексе Fastest-3D: Учебное пособие. - М.:

Изд-во МАИ, 2007. – 116 с.: ил.

ISBN 978-5-7035-1854- Приведены основные направления моделирования процессов в турбулентных течениях на примере трехмерного программного комплекса “FASTEST-3D”.

Для студентов инженерных факультетов, а также для читателей, интересующихся моделированием турбулентных течений.

Рецензенты: Dr.-Ing. А. Мальцев;

Dr.-Ing. Д. Горынцев ISBN 978-5-7035-1854- © Московский авиационный институт (государственный технический университет) © Юн А. А., Крылов Б. А., Тем. план 2007, доп.

Юн Александр Александрович Крылов Борис Анатольевич Расчет и моделирование турбулентных течений с теплообменом, смешением, химическими реакциями и двухфазных течений в программном комплексе Fastest-3D Редактор М. С. Винниченко Подписано в печать 29.08. Бум. офсетная. Формат 60х84 1/16 Печать офсетная.

Усл. печ. л. 6,74. Уч.-изд. л. 7,25. Тираж 250 экз.

Зак. 3736/2239. С.629.

Издательство МАИ «МАИ», Волоколамское ш., д. 4, Москва, А-80, ГСП-3 125993.

Типография Издательства МАИ «МАИ», Волоколамское ш., д. 4, Москва, А-80, ГСП-3 125993.

Оглавление Предисловие Условные обозначения 1. Теоретические основы 1.1. Основные уравнения 1.2. Двухфазные течения 1.3. Горение 1.4. Турбулентность 2. Моделирование турбулентных течений 2.1. Прямое численное моделирование (DNS) 2.2. Моделирование крупных вихрей (LES) 2.3. Моделирование на базе осредненных уравнений Рейнольдса (RANS) 2.3.1. Предположение Буссинеска 2.3.2. Алгебраические модели 2.3.3. Модели с одним дифференциальным уравнением 2.3.4. Модели с двумя дифференциальными уравнениями 2.3.5. Нелинейные модели и алгебраические модели рейнольдсовых напряжений 2.3.6. Дифференциальные модели рейнольдсовых напряжений 2.3.7. Моделирование процессов теплообмена и смешения 2.3.8. Моделирование двухфазных течений 2.3.9. Моделирование горения 2.3.10. Моделирование пристеночных течений 3. Численные методы 3.1. Метод конечных объемов 3.2. Преобразование координат 3.3. Дискретизация диффузионного члена 3.4. Дискретизация конвективного члена 3.5. Дискретизация временного члена 3.6. Дискретизация источникового члена 3.7. Методы интерполяции 3.8. Общая система уравнений 3.9. Связь уравнений для скорости и давления 3.10. Решение системы линейных уравнений 4. Программный комплекс FASTEST-3D 4.1. Установка и структура FASTEST-3D 4.2. Связь сеточных файлов с FASTEST-3D 4.3. Выбор моделей, схем дискретизации и других параметров 4.4. Введение граничных условий, запуск программы 4.5. Вывод результатов 4.6. Пример расчета Библиографический список Приложение ПРЕДИСЛОВИЕ Появление быстродействующих ЭВМ в последние 30 лет резко изменило характер применения основных принципов теоретической гидромеханики и теплопередачи при решении инженерных задач. Наряду с развитием традиционных методов, таких как аналитические и экспериментальные методы, бурно развивается третий метод исследований – вычислительная гидрогазодинамика (Computational Fluid Dynamics – CFD).

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

Данное методическое пособие призвано ввести читателей в практический курс вычислительной гидрогазодинамики на базе трехмерного расчетного комплекса FASTEST-3D (Flow Analysis Solving Transport Equation Simulating Turbulence 3 Dimensional), любезно предоставленного техническим университетом Дармштадта (TU Darmstadt). FASTEST-3D позволяет охватить полный спектр численных расчетов для описания ламинарных и турбулентных течений с теплообменом, смешением, химическими реакциями, двухфазных течений в сложных инженерных конфигурациях, кроме этого, открытый код предоставляет неограниченные возможности как добавления новых моделей для описания физики процессов, так и улучшения математической составляющей. В настоящее время FASTEST-3D поддерживается и расширяется рядом технических вузов Германии в рамках проекта SFB, финансируемой правительством и ведущими транснациональными компаниями, такими как Siemens, Rolls-Roys, Alstom. Участие ведущих технических вузов Германии позволяет программному коду иметь последние наработки в области турбулентных течений и горения. Знание данного кода, а также методов разработки позволяет просто перейти к использованию таких коммерческих кодов, как Fluent, CFX, ANSYS и других не только в рамках пользователя, но и разработчика.

FASTEST-3D требует установки операционной системы, базирующейся на Linux, также требуется установка программ создания геометрии и сетки и программы для визуализации полученных результатов. Требования к ресурсам компьютера зависят от сложности конфигурации и выбранных моделей и могут варьироваться в пределах от Pentium 100 (32Mb) до суперкомпьютера.

Авторы благодарят за помощь в рецензировании методического пособия Dr.-Ing. Мальцева А., Dr.-Ing. Горынцева Д. Отдельная благодарность техническому университету Дармштадта (TU Darmstadt) в лице Prof. Dr.-Ing. M. Schafer, Prof. Dr.-Ing. J. Janicka, Prof. Dr.-Ing. A. Sadiki, Dr.-Ing. D. Sternel.

УСЛОВНЫЕ ОБОЗНАЧЕНИЯ Латинские Символ Размерность Определение I, II, III, IV,V инварианты первого, второго и т.д.

порядка компоненты тензора анизотропии м aij с2 рейнольдсовых напряжений модельные коэффициенты Ci, ci кориолисовы силы Н C ji удельная теплоемкость при Дж cp ( кгК ) постоянном давлении м диаметр частиц dp диффузионный перенос пассивного D скаляра весовая функция F1, F силы Н Fi демпфирующие функции f, fi ускорение свободного падения м gi с турбулентная кинетическая энергия м k с + безразмерная турбулентная k кинетическая энергия, k + = k U м характерный масштаб l турбулентности статистическое давление Па p член генерации турбулентности Pk число Прандтля Pr турбулентное число Прандтля Prt газовая постоянная Дж R ( кгК ) число Рейнольдса Re число Рейнольдса для частиц Re p турбулентное число Рейнольдса Ret молекулярное число Шмидта Sc компоненты тензора скоростей де Sij c формации температура К T с время t вектор скорости, компоненты u, ui м с скорости динамическая скорость, w / w м u с + безразмерная скорость u компоненты тензора рейнольдсовых uiu j м с2 напряжений компоненты потока скаляра u j объем м V м декартовы координаты xi м расстояние от стенки y y+ обезразмеренное расстояние от стенки, u y / Греческие Символ Размерность Определение i, i, i, i модельные коэффициенты молекулярный коэффициент диффузии пассивного скаляра, ij символ Кронекера скорость диссипации м с2 турбулентной кинетической энергии 0 скорость диссипации м с2 турбулентной кинетической энергии на стенке + безразмерная скорость диссипации турбулентной кинетической энергии, + = U i м локальные координаты k м колмогоровский масштаб длины коэффициент теплопроводности Дж ( сКм ) ij тензор перераспределения величина пассивного скаляра константа Кармана, t коэффициенты динамической кг ( мс ) молекулярной и турбулентной вязкости, t коэффициенты кинематической м с2 молекулярной и турбулентной вязкости плотность кг м площадь м турбулентное число Шмидта или Прандтля с турбулентный масштаб времени k с колмогоровский масштаб времени ij компоненты тензора м с2 рейнольдсовых напряжений ij компоненты тензора вектора c завихренности k, i скорость вращения c удельная скорость диссипации переменная Операторы Оператор Определение ( ) осреднение Рейнольдса ( ) осреднение по Фавру div ( ) дивергенция grad ( ) градиент Нижний индекс Оператор Определение ( )b сгоревшие компоненты ( ) F топливо ()O окислитель ( ) P продукт горения ( )u несгоревшие компоненты Аббревиатуры Аббревиатуры Определение AS Абе, Шуга BML Брэй, Мосс, Либби CDS центрально-разностная схема (Central Difference Scheme) CFD вычислительная гидрогазодинамика (Computational Fluid Dynamics) CLS Крафт, Лаундер, Шуга CV контрольный объем (Control Volume) DES изолированное моделирование крупных вихрей (Detached Eddy Simulation) DH Дейли, Харлоу DNS прямое численное моделирование (Direct Numerical Simulation) EARSM явная алгебраическая модель рейнольдсовых напряжений (Explicit Algebraic Reynolds Stress Model) EASFM явная алгебраическая модель переноса скаляра (Explicit Algebraic Scalar Flux Model) EDC концепция вихревой диссипации (Eddy Dissipation Concept) EVM модель вихревой вязкости (Eddy Viscosity Model) GGDH основная градиенто-диффузионная гипотеза (General Gradient-Diffusion Hypothesis) GL Гибсон, Лаундер GS Гатски, Специале IARSM неявная алгебраическая модель рейнольдсовых напряжений (Implicit Algebraic Reynolds Stress Model) KM Ким, Моин LES моделирование крупных вихрей (Large Eddy Simulation) LRR Лаундер, Рис, Роди LU нижняя- верхняя декомпозиция (Lower Upper decomposition) MPI Message Passing Interface PDF функция плотности вероятности RANS моделирование на базе осредненных уравнений Рейнольдса (Reynolds Averaging based Numerical Simulations) RSM модель рейнольдсовых напряжений (Reynolds Stress Model) RSFM модель рейнольсового потока скаляра (Reynolds Scalar Flux Model) RST тензор рейнольдсовых напряжений (Reynolds Stress Tensor) SJ Съерген, Йоханссон SGS Шуга, Гатски, Специале SIMPLE полуявный метод для связи уравнений через давление (Semi Implicit Method for Pressure Linked Equations) SIP строгая явная процедура (Strongly Implicit Procedure) SMC замыкание второго уровня (Second Moment Closure) SSG Специале, Саркар, Гатски SST перенос сдвиговых напряжений (Shear Stress Transport) UDS разностная схема «вперед» (Upwind Difference Scheme) URANS моделирование на базе нестационарных осредненных уравнений Рейнольдса (Unsteady RANS) WJ Валлин, Йоханссон WWJ Виркстром, Валлин, Йоханссон ZFK Зельдович-Франк-Каменецкий Глава 1. ТЕОРЕТИЧЕСКИЕ ОСНОВЫ Уравнения динамики сплошных сред основаны на универсальных законах сохранения: массы, импульса и энергии.

1.1. Основные уравнения Уравнение неразрывности В дифференциальной форме уравнение неразрывности или сохранения массы [4] записывается в виде ui + = 0, (1.1) t xi I II где первый член описывает нестационарность потока, второй член представляет собой конвективный перенос.

Для несжимаемых жидкостей уравнение неразрывности записывается в виде ui = 0. (1.2) xi Уравнение изменения количества движения Уравнение количества движения записывается следующим образом ui ui u j p ij + gi, + = + (1.3) t x j xi x j V I III II IV где первый член описывает нестационарность потока, второй конвективный перенос, третий и четвертый члены- поверхностные силы (градиент давления и молекулярную диффузию), пятый массовые силы (гравитацию). Для ньютоновских жидкостей тензор напряжений по гипотезе Стокса [54] представляется в виде ui ui 2 uк ij =, + (1.4) xi x j 3 xк ij где ij - символ Кронекера ( i = j ij = 1, i j ij = 0 ).

С учетом последнего соотношения уравнение изменения количества движения записывается в следующем виде u u 2 uк ui ui u j p i + i + gi. (1.5) + = + xi x j 3 xк ij t x j xi x j Проведя ряд преобразований с уравнениями (1.2) и (1.5) можно получить уравнение движения в форме Навье–Стокса [1].

Для несжимаемых жидкостей и при постоянной вязкости уравнение Навье–Стокса [4] записывается в виде 1 p 2ui ui ui +uj = + + gi. (1.6) xi x j x j t x j Уравнения баланса энергии и пассивного скаляра Уравнение баланса энергии записывается в виде qi ij ui c pT u i c pT + = + +S, (1.7) t x j xi x j V I III IV II первый член описывает нестационарность потока, второй конвективный перенос, третий- перенос за счет теплопроводности и T определяется законом Фурье для переноса энергии qi =, xi четвертый- диссипацию энергии, пятый- приток или убытие энергии за счет химических реакций, радиации и т.д.

Для расчета течений со смешением, теплообменом или химическими реакциями используется также уравнение переноса пассивного скаляра (температура, массовая доля), которое выглядит следующим образом ui + = +S. (1.8) D t xi xi xi Коэффициент диффузии D для массы записывается как соотношение вязкости и числа Шмидта D=. (1.9) Sc Для замыкания системы к уравнениям, полученным из упомянутых выше законов сохранения, следует добавить соотношения, устанавливающие связь между свойствами жидкости.

Примером такого соотношения может быть уравнение состояния, связывающее термодинамические параметры жидкости: давление, плотность и температуру p = RT. (1.10) За исключением некоторых случаев аналитическое решение данных уравнений невозможно, в частности, в случае турбулентных течений вследствие их нестационарного и случайного характера. Альтернативой аналитическому решению являются численные методы решения уравнений. В инженерных расчетах наиболее часто используется подход, предложенный Рейнольдсом, в котором мгновенные значения параметров потока представляются в виде суммы осредненной величины (например, статистически стационарного течения по времени) и ее пульсационной составляющей. Вследствие применения процедуры осреднения появляются новые незамкнутые корреляции (так называемые рейнольдсовые напряжения), которые требуют моделирования.

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

1.2. Двухфазные течения Двухфазные течения играют важную роль в технических приложениях. Большинство процессов, происходящих в камерах сгорания в тепловых машинах, в частности воздушно- реактивных двигателей, являются двухфазными (керосин–воздух).

Для описания двухфазных течений наиболее часто используют два метода: дискретно–траекторный Эйлер–Лагранж (Euler–Lagrange) и многожидкостный Эйлер–Эйлер (Euler–Euler).

В первом методе рассчитываются траектории и характеристики индивидуальных частиц в определенные промежутки времени. Т.е.

для частиц используется метод Лагранжа, а для жидкой фазы метод Эйлера. Во втором методе дискретная фаза (частицы) рассматривается как сплошная среда (континуум), а общее течение – как смесь двух или более жидкостей. На рис. 1.1 графически показано различие двух методов. Слева координатные системы для двух фаз совмещены (Эйлер–Эйлер), справа координатная система для твердой фазы привязана к частице (Лагранж–Эйлер).

Рис. 1.1. Метод Эйлера–Эйлера (слева), метод Лагранжа–Эйлера (справа) для двухфазных течений Кроме вышеперечисленного, моделирование двухфазных течений требует дополнительных моделей для описания столкновений частиц, коагуляции, распада, испарения и т.д. [12], [14], [49].

1.3. Горение Течения в камерах сгорания тепловых машин являются не только двухфазными, как было сказано выше, но и сопровождаются химическими реакциями (горением). В рамках данного учебного пособия рассмотрено только однофазное горение. Двухфазное горение подробно изложено в [63].

Определение горения в терминах макроскопической кинетики дано Франк–Каменецким в одной из первых работ в этой области [6]: «Горением, называется протекание химической реакции в условиях прогрессивного самоускорения, связанного с накоплением тепла или катализирующих продуктов реакции».

Макроскопическая теория кинематики горения изучает роль конвекции, диффузии и теплового переноса и их влияние на химические процессы. Горение можно условно разделить на три категории: предварительно смешанное (кинетическое), диффузионное и частично смешанное.

Химическая кинетика Химическая кинетика как элемент моделирования горения предполагает, что химическая реакция происходит без учета эффектов течения (конвекции, диффузии). Переход реагирующих веществ в продукты сгорания можно упрощенно описать одношаговой реакцией F F + OO = P P, (1.11) где количество молей: F - топлива;

O - окислителя;

P - продуктов сгорания. Однако в реальности механизм окисления топлива представляется сложнее и может описываться сотнями элементарных реакций. С учетом реальной химической системы для N компонент и M элементарных химических реакций можно записать следующую систему уравнений:

N N kj k = kj k, j = 1, M, (1.12) k =1 k = где k - символ для обозначения компоненты k и kj, kj молярные стехиометрические коэффициенты компоненты k в реакции j. Система (1.12) должна удовлетворять уравнению сохранения массы N kjWk = 0, j = 1, M, (1.13) k = где kj = kj kj и Wk обозначает молекулярный вес компоненты i k. Скорость kj образования компоненты k в реакции j записывается следующим образом:

i kj = r jWk kj, (1.14) i где r j - скорость протекания реакции j. Общая скорость k образования компоненты k записывается как сумма всех скоростей M элементарных реакций:

M m i i k = kj =Wk r j kj. (1.15) j =1 j = Из (1.15) и (1.13) можно записать N Ni M k = r j Wk kj. (1.16) j =1 k =1 j = Скорость протекания реакции r j находится как разность между скоростями протекания прямой и обратной реакции:

r j = r f, j rb, j, (1.17) kj kj Yk Yk N N где r f, j = K f, j k = ;

rb, j = Kb, j k =.

Wk Wk Скорости прямой и обратной реакции K f, j, Kb, j находятся из закона Аррениуса Ef,j nf, j K f, j = A f, jT exp. (1.18) RT В (1.18) три эмпирических параметра: предэкспоненциальный n коэффициент A f, j, температурный коэффициент T f, j и энергия активации E f, j. Скорость обратной реакции Kb, j находится обычно из скорости прямой реакции и константы равновесия S 0 H N kj p = K f, j a j j exp, k = (1.19) Kb, j RT RT R где pa соответствует окружающему давлению;

S 0, H 0 j j изменение энтропии и энтальпии во время реакции.

Задачей химической кинетики является нахождение наиболее важных компонент и реакций в уравнении (1.12), а также и определение эмпирических параметров в (1.18). Ряд различных схем решения данных задач приведен в [6]. Как упоминалось выше, горение делят на три вида. Рассмотрим каждый из них подробнее.

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

смеси и зона «бедной» смеси».

В отличие от кинетического горения, которое будет рассмотрено ниже, диффузионное пламя не может распространяться против течения. Кроме этого, к диффузионному пламени нельзя применить определение «толщины» пламени.

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

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

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

Простейшая модель диффузионного горения основана на введении так называемой массовой доли (mixture fraction), записываемой в следующем виде:

z z,O =, (1.20) z, F z,O где z обозначает массовую долю элемента.

Для массовой доли решается уравнение баланса, аналогичное уравнению для скаляра (1.19), без источникового члена:

ui + =. (1.21) t xi xi xi Простейшее описание диффузионного горения, представленного выше, было предложено Швабом и Зельдовичем [6]. Такой подход позволяет описать полную структуру диффузионного пламени с помощью только одной переменной.

Рис. 1.2. Одномерное описание диффузионного горения На рис. 1.2 представлена идеализированная одномерная структура диффузионного горения в значении одной переменной. В этом случае концентрация компонент, температура и плотность является функцией массовой доли: Yk = Yk ( ), T = T ( ), = ( ), в то время как массовая доля является функцией от времени и координат в пространстве = ( t, xi ). Для нахождения первых трех функций используется ряд различных моделей химической кинетики [6], [63]. Более детально моделирование диффузионного горения описано в [6].

Кинетическое горение В отличие от диффузионного горения, при кинетическом горении топливо и окислитель уже смешаны, перед тем как попасть в зону горения, поэтому иногда его называют предварительно смешанным горением. Схематически это показано на рис. 1.3 для идеализированного одномерного горения. Впервые такое представление было введено Зельдовичем и Франк-Каменецким [6].

Характерно увеличение температуры газа T вдоль оси x от минимального значения T0 (реагенты) до максимальной температуры Tmax. Фронт пламени, определенный как зона значительного изменения температуры, состоит из двух основных частей: предварительная зона, где диффузия тепла и масс наиболее интенсивна, пока химические реакции еще не начаты, и зона реакций, где скорость реакций резко возрастает, а потом падает;

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

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

Ниже изложены основные принципы, сформулированные Зельдовичем, Франк-Каменецким и Карманом (ZFK) для описания кинетического горения:

- зона реакций расположена в зоне высоких температур и имеет температуру близкую к Tmax ;

- толщина зоны реакций l приблизительно на порядок меньше толщины фронта пламени l f ;

- фронт пламени распространяется в направлении реагентов со скоростью s. Эта скорость пропорциональна квадратному корню скорости реакции (1.18) при T = Tmax и термической диффузии [6] E. (1.22) s~ exp cp RTmax Данное уравнение показывает природу распространения пламени, его связь с кинетикой выделения тепла и с переносом тепла из горячих слоев в холодные;

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

Частично смешанное горение Во многих камерах сгорания условия течения таковы, что горение невозможно классифицировать ни как диффузионное, ни как кинетическое в чистом виде. Такое смешанное горение называется частично смешанным. Петерс в [40] определил его так:

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

Рис. 1.4. Пример смешанного горения На рис. 1.4 схематично приведен один из примеров частично смешанного горения – так называемое струйное горение или тройное горение. Из центрального сопла подается топливо, снаружи подсасывается воздух. На границе стехиометрического смешения происходит диффузионное горение. При достаточно высоких скоростях истечения топлива из сопла возможно растяжение, а потом и разрушение структуры диффузионного горения. В этой области реагенты смешиваются друг с другом без горения и пламя стабилизируется («сдувается») ниже по течению струи. Стабилизация пламени происходит в точках, где скорость течения и скорость кинетического горения находится в равновесии и характеризуется длиной «отрыва». Таким образом, одновременно около этих точек образуется три зоны: диффузионного горения, «бедного» кинетического горения, «богатого» кинетического горения.

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

1.4. Турбулентность Большинство ньютоновских течений, представляющих интерес в технических приложениях, являются турбулентными.

Питер Брэдшоу во вступлении книги «Турбулентность» написал:

«один неоспоримый факт о турбулентности это то, что турбулентность наиболее сложное из движений жидкости».

Следуя Рейнольдсу, возможно охарактеризовать состояние течения числом Рейнольдса, определяемым по следующей формуле:

l0 u Re =. (1.23) Здесь l0 - характерная длина;

u0 - характерная скорость;

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

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

Большие вихри переходят в маленькие вихри;

в свою очередь, маленькие вихри, размеры, которых немного превышают длину молекулярного пробега, рассеивают энергию в тепло посредством молекулярной вязкости. Этот процесс впервые был математически описан Колмогоровым [33] и назван энергетическим каскадом.

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

Рис. 1.5. Генерация турбулентности за решеткой Итак, турбулентному течению присущи следующие свойства:

- случайный характер во времени и пространстве, - нестационарность, - трехмерность, - диссипативность, - вихревой характер.

Описание турбулентности возможно с помощью вероятностных законов с использованием функции плотности вероятности [40]. Однако в рамках данного пособия требуются только осредненные параметры (скорость, давление, плотность и т.д.) течения, поэтому в дальнейшем будут рассмотрены следующие методы описания турбулентности: DNS, LES, RANS и URANS.

Глава 2. МОДЕЛИРОВАНИЕ ТУРБУЛЕНТНЫХ ТЕЧЕНИЙ 2.1. Прямое численное моделирование (DNS) Прямое численное моделирование (DNS) предполагает решение полных нестационарных уравнений Навье–Стокса и уравнения неразрывности. Это означает, что не требуется дополнительного моделирования и происходит учет всех эффектов, присущих течению.

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

3 k =, (2.1) а временной шаг – колмогоровскому масштабу времени:

k =. (2.2) Ниже, в табл. 2.1, показана зависимость времени расчета от числа Рейнольдса для канала [38].

Табл. 2.1.

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

5103 5104 5105 5106 Re 200 68 дн. 444 дн. 610 лет Mflop/s 1 Tflops/s 13 дн. 88 дн. 122 года Очевидно, что DNS позволяет получать результаты при современном развитии вычислительной техники только при малых числах Рейнольдса. С практической точки зрения, статистика, полученная с DNS, может быть использована для тестирования и калибровки моделей, базирующихся на осредненных уравнениях Рейнольдса.

2.2. Моделирование крупных вихрей (LES) При моделировании крупных вихрей крупные вихри рассчитываются, а мельчайшие вихри подсеточного масштаба моделируются. Основной предпосылкой такого подхода является то, что наибольшие вихри несут максимум рейнольдсовых напряжений и должны быть рассчитаны. Мелкие же вихри содержат невысокие значения рейнольдсовых напряжений, кроме того, мелкомасштабная турбулентность близка с изотропной и имеет близкие к универсальным характеристики, в большей степени подающиеся моделированию.

Так как LES включает моделирование мельчайших вихрей, характерные размеры ячеек могут быть намного больше, чем колмогоровский масштаб длины, а временные шаги могут быть выбраны гораздо большими, чем они возможны в DNS. Таким образом, требования к вычислительным ресурсам LES гораздо ниже, чем DNS. С бурным развитием вычислительной техники сфера применения LES значительно возросла. Существует мнение, что LES в ближайшее время вытеснит RANS. Такая точка зрения является довольно спорной. До сих пор в LES не решена проблема пристеночных течений. Очевидно, что вблизи стенки вихри малы и анизотропны;

сеточные и временные шаги, требуемые для LES, падают до величин, характерных для DNS. Cуществующие решения, такие как анизотропные фильтры и динамические процедуры, пока не дают удовлетворительных результатов. Одним из решений данной проблемы является комбинирование LES и RANS – так называемое изолированное или присоединенное LES DES. В DES в пристеночной зоне используется RANS, а далее используется LES. Но и здесь существуют проблемы, например, сшивание зон, т.е. передача информации между зонами. Несмотря на существующие проблемы моделирование крупных вихрей является перспективным направлением в исследовании турбулентных течений в настоящее время.

2.3. Моделирование на базе осредненных уравнений Рейнольдса (RANS) Для технических приложений более оправданным является метод, базирующийся на решении осредненных уравнений Рейнольдса.

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

( xi, t ) = ( xi,t ) + ( xi t ). (2.3) В интегральной форме осреднение по времени можно записать так:

1t ( i ) = lim ( xi, t ) dt. (2.4) t t На рис. 2.1 приведен пример для скорости при стационарном и нестационарном течении.

Рис. 2.1. Осреднение для: (а) – стационарного течения;

(б) – нестационарного течения Применяя данное осреднение, возможно, получить осредненные уравнения сохранения массы, импульса и скаляра (энергии) [62].

Для несжимаемой жидкости они записываются в следующем виде:

u i = 0;

(2.5) xi ui uiu j u i ui u j 1 p + = + + gi ;

(2.6) xi x j xi t x j x j + u j D u j + S, = (2.7) x j x j t x j где uiu j - составляющие тензора напряжений Рейнольдса, а u j составляющие потока скаляра.

Составляющие тензора напряжений Рейнольдса являются дополнительными шестью неизвестными;

кроме того, добавляется три неизвестных в уравнении для скаляра (энергии). Таким образом, приведенная выше система уравнений является незамкнутой. Вопросы и проблемы нахождения вышеприведенных неизвестных рассмотрены в следующей главе.

При расчете двухфазных течений методом Эйлера решаются дополнительно к уравнениям (2.5) - (2.7) аналогичные уравнения для частиц. Полученные решения суммируются. При методе Лагранжа учет частиц происходит добавлением источниковых членов в уравнения (2.5) - (2.7) и решение проходит итерациями между двумя фазами.

Турбулентное течение с горением характеризуется значительными флуктуациями плотности. В этом случае осредненные уравнения, приведенные выше, будут иметь дополнительные неизвестные корреляции, включающие флуктуации плотности. Чтобы избежать этого, используют осреднение по плотности, предложенное Фавром в 1969 г. [62], записываемое следующим образом:

f f=. (2.8) Для скорости u = ui + ui ;

(2.9) для скаляра = i + i. (2.10) Осредненные по Фавру уравнения сохранения записываются в следующем виде:

ui + = 0;

(2.11) t xi ui u j ui u p i + uiu j + gi ;

+ = (2.12) xi x j t x j x j u j D u j +.

+ = (2.13) x j t x j x j Кроме проблемы замыкания осредненных уравнений (2.11) (2.13), при горении важную роль играет также взаимосвязь между турбулентной и химической составляющей горения (turbulence chemistry interaction). Для решения скалярного уравнения (2.13) требуется найти значение источникового члена. Учитывая нелинейный характер отношений между мгновенными скоростями реакций, температурой и химическими компонентами в уравнениях (1.14) - (1.19) (глава 1), значение химического источникового члена нельзя однозначно записать в виде функции значений скаляров:

(T, Y ), (2.14) так же, как и для плотности (T, Y ). (2.15) Источниковый член, учитывающий турбулентность, можно найти, используя функцию плотности вероятности (PDF), записываемой в общем виде как + ( ) P ( ) d, = (2.16) где – переменная, зависящая нелинейно от, например скорость реакции от температуры;

P ( ) – функция плотности вероятности.

Имея зависимость ( ) и функцию плотности вероятности P ( ), не составляет труда найти значение простым интегрированием (2.16). Однако функция плотности вероятности является неизвестной и требует моделирования. Для ее нахождения возможно использовать уравнение переноса [35], что сопряжено с аналогичными проблемами замыкания, как и осредненные уравнения Рейнольдса. Прямое же решение, так называемый метод Монте-Карло требует значительных компьютерных ресурсов.

Другим путем нахождения функции плотности вероятности является предположение о ее форме. В качестве распределения плотности вероятности обычно используется функция, реже функция Гаусса [35]. Подробно о взаимосвязи турбулентной и химической составляющих горения, в частности, о моделировании функции плотности вероятности для горения изложено в работах [6], [40], [63].

2.3.1. Предположение Буссинеска Модели турбулентности, используемые в инженерных приложениях в настоящее время, основаны на концепции вязкости и турбулентной диффузии. В 1877 г. Буссинеск [13] выдвинул предположение, что рейнольдсовы напряжения могут быть связаны со скоростью средней деформации через турбулентную вязкость.

Для тензора рейнольдсовых напряжений это дает ui u j uiu j = t k i.

+ (2.17) x j xi Данное уравнение не вводит модели турбулентности, а только характеризует структуру такой модели. При этом основной задачей является задание функции турбулентной вязкости t. В отличие от коэффициента молекулярной вязкости, коэффициент t определяется состоянием турбулентного течения и не связан со свойствами жидкости. Значение t может значительно изменяться от точки к точке в пространстве, в зависимости от характера течения. Согласно предположению Буссинеска, перенос количества движения происходит аналогично переносу за счет молекулярного движения.

По аналогии с кинетической теорией газов можно ожидать, что турбулентная вязкость с достаточной точностью представляется в виде t ut l, (2.18) где ut - характерная пульсационная составляющая скорости;

l характерный линейный масштаб турбулентности. Проблема состоит в оценке и определении ut и l.

Понятие турбулентной вязкости имеет ряд недостатков.

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

uiu j ij.

aij = (2.19) k Как следует из гипотезы Буссинеска, тензор анизотропии рейнольдсовых напряжений пропорционален тензору скоростей деформации осредненного течения:

u u j 2 t i + aij =. (2.20) x j xi k Вышеизложенное предположение не выполняется даже во многих простых течениях, например течение в трубе, вращающейся вокруг своей оси, не говоря уже о более сложных течениях. С другой стороны, во многих случаях, особенно при анализе течений, в которых основное влияние на осредненное движение оказывает лишь одна из компонент тензора рейнольдсовых напряжений – напряжение сдвига xy, нарушение гипотезы Буссинеска не приводит к заметным погрешностям.

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

2.3.2. Алгебраические модели Наиболее простыми моделями, определяющими турбулентную вязкость t, являются алгебраические модели, в которых связь между турбулентной вязкостью и параметрами осредненного потока задается алгебраическими соотношениями.

Алгебраическая модель для описания распределения t впервые была предложена Прандтлем в 1925 г. и известна как модель смешения [62].

Модель Прандтля записывается в виде vx t = lm, (2.21) x где lm - длина пути смешения, определяемая эмпирически.

В пограничном слое полагают, что lm = кy, (2.22) где к 0.39 - число Кармана;

y - расстояние от стенки.

Популярными алгебраическими моделями являются:

двухслойные модели Себеси–Смит [48], Балдвин–Ломакс [11], модель с половинным уравнением Джонсон–Кинг [27] и т.д.

Оценка применимости алгебраических моделей турбулентности детально обсуждена в работе Вилкокса [62] К достоинствам алгебраических моделей можно отнести:

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

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

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

2.3.3. Модели с одним дифференциальным уравнением Чтобы преодолеть ограниченность гипотезы пути смешения и алгебраических моделей вообще, были разработаны модели турбулентности, позволяющие учитывать влияние нелокальных эффектов (эффектов переноса) путем решения дифференциального уравнения для ut или l. С физической точки зрения, для величины ut наиболее подходящим оказывается масштаб k, где k кинетическая энергия турбулентных пульсаций. Если такой масштаб использовать для определения турбулентной вязкости в (2.18), то получается выражение Колмогорова–Прандтля:

* t = C k L, (2.23) где C - эмпирическая функция местного турбулентного числа Рейнольдса. Для масштаба энергии турбулентных пульсаций записывается уравнение переноса, модельный вид которого t k u j ui u j k k k + + +uj = + CD, (2.24) xi k x j t xi x j t x j x j L где C, CD, k - модельные константы. Значение линейного масштаба L определяется из простых эмпирических соотношений аналогично для пути смешения.

Существуют аналогичные модели с одним уравнением, где используются уравнения переноса для турбулентного трения, (Брэдшоу, Феррис, Атвел [15]), для турбулентной вязкости (Ни, Коважная [39]) и ряд других моделей [62].

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

2.3.4. Модели с двумя дифференциальными уравнениями Более универсальными моделями в инженерных расчетах турбулентных потоков являются модели с двумя дифференциальными уравнениями. Первая такая модель была предложена Колмогоровым (1942) [33]. Эта модель содержит уравнение переноса кинетической энергии турбулентности k и удельной скорости диссипации энергии.

k - модель Ниже приведена одна из распространенных моделей k типа, предложенная Вилкоксом [62]:

k ( ) k k u i *k + + * t + u j = ij ;

(2.25) t x j x j x j x j ui + ( +t ) ;

+ u j = ij (2.26) t xj k xj xj xj k t = ;

(2.27) u u j ij = uiuj = t i kij.

x j xi (2.28) Модельные константы:

* = 9 /100;

= 3 / 40;

= 5 / 9;

* = 1/ 2;

= 1/ 2.

k - модель Наиболее популярной моделью с двумя дифференциальными уравнениями является k - модель, предложенная Чоу (1945) [17] и получившая дальнейшее развитие в исследованиях Лаундера – Джонса (1972) [26].

k k k u i + u j = ij + + t ;

(2.29) k x j t x j x j x j 2 t ui + u j = c1 ij c 2 + + ;

(2.30) k xj xj t xj k xj k t = C ;

(2.31) ui u j ij = uiuj = t kij.

(2.32) x j xi Модельные константы:

c 1 = 1.44;

c 2 = 1.92;

C = 0.09;

k = 1.0;

= 1.3.

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

Основываясь на этом, Ментер (1993) предложил модель [37], сочетающую в себе указанные сильные стороны k и k моделей. Для этого k модель переформулируется в терминах k и, а затем в полученные модельные уравнения вводится весовая функция F1, обеспечивающая плавный переход от k модели в пристеночной области к k модели вдали от стенки. Таким образом, модель Ментера записывается путем суперпозиции моделей k и k, помноженных соответственно на весовую функцию F1 и (1 F1 ). Функция F1 конструируется таким образом, чтобы быть равной единице на верхней границе пограничного слоя и стремится к нулю при приближении к стенке. Кроме того, Ментер видоизменил стандартную связь между k, и турбулентной вязкостью t. В эту связь был введен специальный ограничитель (SST), обеспечивающий переход от нее к известной формуле Брэдшоу [15], согласно которой турбулентное напряжение пропорционально кинетической энергии турбулентности uiu j = 0.31k. Этот прием, получивший название SST (shear stress transport), в дальнейшем широко применялся и в других моделях турбулентности с двумя уравнениями, например Чена [16].

Ниже представлена базовая двухслойная модель Ментера:

k k k u i ( + k t ) ;

(2.33) * k + + u j = ij t x j x j x j x j ui + ( +t ) + + u j = ij t xj k xj xj xj (2.34) 1 k + 2(1 F )2 ;

xj xj k t = ;

(2.35) u u uj = t i + j k ij.

ij = ui (2.36) x j xi Обозначая обобщенным параметром 1 набор констант оригинальной модели k с индексами 1 и соответственно аналогичный набор констант трансформированной k модели, получаем = F11 + (1 F1 )2. (2.37) Используются следующие константы.

1. Модельные константы k модели (Вилкокс):

* = 0.09;

1 = 0.075;

k 1 = 0.5;

1 = 0.5;

= 0.41;

1 = 1 / * 1 2 / *.

2. Модельные константы стандартной k модели:

* = 0.09;

2 = 0.0828;

k 2 = 1;

2 = 0.856;

= 0.41;

2 = 2 / * 2 2 / *.

Весовая функция определяется следующим образом:

k 500 4 2 k arg1 = min max * ;

2 ;

;

(2.38) y y CDk y 1 k CDk = max 2 2,10 ;

(2.39) x j x j F1 = tanh(arg1 ). (2.40) Согласно предложению Брэдшоу, напряжение сдвига в пограничном слое пропорционально кинетической энергии = a1k (2.41) с постоянной a1 [15]. С другой стороны, в моделях с двумя уравнениями напряжение сдвига рассчитывается следующим образом:

= Ht, (2.42) u где =.

y Чтобы удовлетворить уравнению Брэдшоу в рамках концепции вихревой вязкости, коэффициент турбулентной вязкости переопределяют следующим образом:

a1k t =. (2.43) Чтобы расширить формулировку вихревой вязкости, приспособленную для свободных сдвиговых слоев для случаев, где предложение Брэдшоу не обязательно применимо, сделана модификация SST- модели для ограниченных стенками потоков.

Для этого применяется смесительная функция F2 :

a1k t =, (2.44) max ( a1w1F2 ) где F2 определяется подобно (2.40).

( ) F2 = tanh arg 2 ;

(2.45) arg 2 = max 2 k / 0,09 y;

2. (2.46) y Новые константы во внутреннем слое:

* = 0.09;

1 = 0.075;

k 1 = 0.85;

1 = 0.5;

= 0.41;

1 = 1 / * 1 2 / * ;

a1 = 0.31.

Во внешнем слое константы не менялись.

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

2.3.5. Нелинейные модели Одним из способов преодолеть недостатки моделей, базирующихся на предположении Буссинеска, является введение нелинейных членов в уравнение Буссинеска. Впервые такой подход был предложен Поупом [41], в дальнейшем он получил развитие в работах Специале [53].

Обычно в нелинейных моделях используют тензор анизотропии рейнольдсовых напряжений вместо обычного тензора рейнольдсовых напряжений aij. Кроме этого, в нелинейных моделях используются тензор деформации скоростей 1 u i u j 1 u i u j Sij = + и тензор завихренности ij =.

2 x j xi 2 x j xi Ниже представлена нелинейная модель Крафт, Лаундер, Шуга (CLS) [18].

Тензор анизотропии рейнольдсовых напряжений вычисляется следующим образом t t 2 1 S II s I + c2 t ( S S ) + c3 t 2 II I + a= S + c k 3 k k ( ) + c4 t2 S 2 S 2 + c5 t2 S 2 2 S IVI + k k + c6 t2 Sij S jk Ski + c7 t2 Sij jk ki. (2.47) Постоянный коэффициент в линейных моделях вычисляется c по следующей формуле:

0. 0. c = 1 exp. (2.48) exp ( 0.75max ( S, ) ) 1. 1 + 0.35 ( max ( S, ) ) Модельные константы:

2 2 c1 = 0.1, c2 = 0.1, c3 = 0.26, c4 = 10c, c5 = 0, c6 = 5c, c7 = 5c.

Аналогичной является нелинейная квадратичная модель Гибсон, Лаундер (GL) [23], более применимая для двухмерных течений.

Здесь тензор анизотропии находится по следующей формуле:

2 ( S S ) + c 3 t S 2 II s I,(2.49) k a = c 1 S c 2 t где ( ) 31 1 + c = ;

(2.50) ( ) 3 + 2 + 6 2 1 + 2 = ( 2 ) / 2, 2 = ( 3 S ) / 8.

2 (2.51) Модельные константы:

1 = 2 / 21, 2 = 1/ 7, 3 = 2 / 7.

Вышеперечисленные нелинейные модели, хотя и позволяют устранить ряд недостатков, присущих предположению Буссинеска, обладают сами таким существенным недостатком, как большое количество модельных констант. Например, CLS модель содержит 7 модельных констант, которые найдены экспериментально.

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

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

Ниже представлена полная форма транспортного уравнения рейнольдсовых напряжений uiuj uk uiu j u u i u j uk i = uiuk + t xk xk xk I II III u u p u u ( ) u j jk + puj ik i j + i + j uj uk + ui p xk x j xi xk V IV u u j 2 i. (2.52) xk xki VI Здесь первый член описывает нестационарность течения;

второй перенос рейнольдсовых напряжений за счет конвекции;

третий генерацию рейнольдсовых напряжений;

четвертый- перенос за счет диффузии;

пятый член- перекрестный член давления;

шестой перенос за счет диссипации.

Аналогично можно записать уравнение переноса для тензора анизотропии:

k Daij uiu j ul uiu j kul uiu j P = 1+ k x t xl k (2.53) ij ij Pij + + Cij, где первый член представляет собой перенос за счет конвекции и описывает нестационарность;

второй- перенос за счет диффузии;

третий и четвертый- генерацию анизотропии рейнольдсовых ' u i ' u i напряжений ( Pij = ui' uk u 'j uk, P = Pij / 2 );

четвертый xk xk перенос за счет диссипации;

пятый- перекрестный член давления;

шестой описывает кориолисовы силы.

Многие неоднородные течения, представляющие инженерный интерес, являются стационарными и удовлетворяют локальному равновесию тензора рейнольдсовых напряжений. В этом случае можно пренебречь конвекционным и диффузионным членами. Это и является основной идеей моделирования ARSM. В математической форме это записывается в виде k Daij uiuj ul uiuj kul = 0. (2.54) t xl k x Таким образом, в ARSM предположение об изотропности t в моделях, использующих предположение Буссинеска, заменяется на предположение о локальном равновесии uiu j, достоверность которого вполне очевидна для значительного числа течений.

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

С учетом вышеизложенного неявная алгебраическая форма уравнения записывается в виде uiu j P Pij ij ij 1 = + + Cij. (2.55) k Член генерации Pij, P = Pij / 2 и кориолисовы силы Cij не требуют моделирования, так как могут быть найдены непосредственно из тензора рейнольдсовых напряжений. Однако тензор скоростей диссипации ij и перераспределения давлений ij требуют моделирования. Моделирование этих членов изложено в следующей главе.

Подставляя уже смоделированные члены, можно получить следующую неявную форму:

7C + P ( ) С1 1 + aij = Sij + 2 aik kj ik akj 15 (2.56) 5 9C2 aik Skj + Sik akj ij aik Ski, 11 или P ( ) A3 + A4 aij = A1Sij + aik kj ik akj A2 aik Skj + Sik akj ij aik Ski.

11( C1 1) 5 9C 88 A1 =, A2 =, A3 =, A4 = Здесь 11( 7C2 + 1) 7C2 + 1 7C2 + 1 7C2 + модельные коэффициенты.

Значения модельных коэффициентов для некоторых моделей приведены в табл. 2.2.

Табл. 2.2.

Модельные константы для различных EARSM A1 A2 A3 A WJ, [60] 1.20 0 1.80 2. LRR, [34] 1.54 0.37 1.45 2. SSG, [53] 1.22 0.47 0.88 2. GS, [23] 1.22 0.47 5.36 1. Решение неявной формы представляет значительные численные трудности вследствие отсутствия диффузионного члена в уравнении переноса для рейнольдсовых напряжений, поэтому для решения этой проблемы используют явную алгебраическую форму. Наиболее общая форма aij в терминах Sij и ij записывается в виде десяти тензорно–независимых групп, в которых все комбинации тензоров более высокого порядка могут быть сокращены по теореме Кайли–Гамильтона. Различные способы сокращения количества тензоральных групп приведены в работах [52], [60], [65]. Согласно [60] тензор анизотропии рейнольдсовых напряжений записывается в виде:

1 a = 1S + 2 S 2 II s I + 3 2 II I + 4 ( S S ) + 3 ( ) + 5 S 2 S 2 + 6 S 2 2 S IVI + 3 (2.57) ( ) + 7 S 2 2 + 2 S 2 VI + 8 S S 2 S 2S 2 + ( ) ( ) + 9 S 2 2 S + 10 S 2 2 2 S 2.

u i u j u i u j, ij = нормированы через Здесь Sij = + 2 x j xi 2 x j xi k масштаб времени = турбулентный и для упрощения используются следующие обозначения: a = aij, S = Sij, = ij.

Коэффициенты i - функции пяти независимых инвариантов S и :

II s = S 2 = Sij Sij, II = 2 = ij ji, III = S 3 = Sij S jk Ski, IV = S 2 = Sij jk ki, (2.58) V = S 2 2 = Sij S jk kj jk.

Решая совместно уравнения (2.56) и (2.57) можно найти неизвестные коэффициенты i и сам тензор анизотропии рейнольдсовых напряжений.

Решение полных уравнений рейнольдсовых напряжений приведено в следующей главе.

2.3.6. Дифференциальные модели рейнольдсовых напряжений В предыдущей главе был приведен способ сведения уравнения переноса рейнольдсовых напряжений (2.52) к алгебраической форме (2.56). В моделях же рейнольдсовых напряжений это уравнение решается шесть раз для каждой компоненты тензора рейнольдсовых напряжений.

Ниже для удобства записано уравнение для переноса рейнольдсовых напряжений в условной форме uiuj + Cij = Pij Dij + ij ij. (2.59) t Здесь Cij описывает конвекцию;

Pij - генерация рейнольдсовых напряжений;

Dij - диффузия;

ij - член перераспределения давления;

ij - диссипация. Три последних члена требуют моделирования.

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

Турбулентная диффузия Dij может моделироваться, как в модели Дейли, Харлоу (DH) [19]:

uiuj k uiu j uk = Cs uk ul, (2.60) xl где эмпирическую константу Cs обычно принимают 0.22.

Моделирование члена перераспределения давления проходит с использованием уравнения Пуассона и делится на две части:

«быструю» и «медленную».

«Медленная» часть чаще всего записывается как линейная функция от анизотропии рейнольдсовых напряжений:

ijs ) ( = C1aij, (2.61) где C1 - модельная константа.

Для моделирования «быстрой» части наиболее часто используется модельное выражение, предложенное Лаундером и др. (Лаундер, Рис, Роди (LRR)) [34]:

ijr ) ( 9C + 6 4 aik Skj + Sik akj akm Skm ij + Sij + = 11 5 3 (2.62) ( ) 7C aik * * akj, +2 kj ik где C2 - модельная константа. В LRR модели приняты следующие модельные константы: C1 = 1.5 и C2 = 0.582 [34].

Диссипацию обычно моделируют с помощью уравнения переноса, аналогичного уравнению, используемому в k модели.

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

2.3.7. Моделирование процессов теплообмена и смешения Как уже упоминалось в главе 1, для расчета течений со смешением, теплообменом или химическими реакциями используется транспортное уравнение для скаляра.

Наиболее простым способом замыкания уравнения переноса является предположение о пропорциональности потока скаляра градиенту осредненного скаляра:

t u j =. (2.63) x j Более сложным является предположение о пропорциональности потока скаляра рейнольдсовым напряжениям, модель Дейли, Харлоу (DH) [19]:

ui = Ct1 t uiu j. (2.64) x j Более точной является модель Ким, Моин (KM) [31], предполагающая более сложное отношение:

uiuk uk u j ui = Ct1 t. (2.65) x j k И, наконец, модель Абе, Шуга (AS) [9], являющаяся комбинацией двух последних моделей:

uiu j uiuk uk u j ui = k t Ct1 + Ct 2. (2.66) x j k k Модельные константы:

Ct1 = 0.22, Ct 2 = 0.22.

Другим подходом для нахождения потока скаляра является решение транспортного уравнения, аналогичного уравнению переноса рейнольдсовых напряжений:

( ui ) + u j ( ui ) = Pi + i i + Di. (2.67) t x j Моделирование членов дифференциального уравнения для потока скаляра описано в [25].

Аналогично, как и в EARSM, возможно привести данное уравнение к алгебраической форме, как в модели Виркстром, Валлин, Йоханссон (WWJ) [61]:

k ( ) ui = 1 c 4 Bij uj uk, (2.68) xk где матрица находится по формуле B G Q1 I G ( cs S + c ) + ( cs S + c ) B= ;

зависимые 31 G GQ1 + Q 2 константы – по формулам cs = 1 c 2 c 3, c = 1 c 2 + c 3 ;

1 1 P G = 2c1 1 ;

параметр параметры, получаемые из r 23 2 инвариантов: Q1 = cs II s + c II, Q2 = cs III s + 2cs c IV.

Модельные константы:

c1 = 4.51, c 2 = 0.47, c 3 = 0.02, c 4 = 0.08.

2.3.8. Моделирование двухфазных течений Одной из характеристик двухфазных течений является число Стокса, определяемое по следующей формуле:

p St =, (2.69) k где p - время релаксации, характеризующее динамическую инерционность частиц, а k - временной интегральный или колмогоровский масштаб турбулентности. Число Стокса характеризует инерционность частиц по отношению к масштабу турбулентности. При большой инерционности частиц или при сильной турбулентности число Стокса стремится к бесконечности.

В этом случае движение частиц и жидкой фазы не коррелируются, т.е. частицы в основном движутся хаотично из одного вихря в другой. С другой стороны, при низких числах Стокса движение частиц коррелируется с движением жидкой фазы. Время релаксации p является функцией от числа Рейнольдса частицы Re p, вычисляемой по формуле W dp Re p =, (2.70) где W - разность скоростей между частицей и жидкой средой;

d p диаметр частицы.

Другой важной характеристикой является концентрация частиц. Различают три вида концентраций частиц: массовую, объемную и счетную. Массовую и объемную концентрацию определяют по следующим формулам:

V p CV = ;

(2.71) V m p Cm =, (2.72) m где индекс p относится к объему и массе частиц;

– к общему объему и массе. Счетная концентрация представляет собой количество частиц.

Метод Эйлера В методе Эйлера (раздел 1.2) осредненные уравнения для движения и теплообмена частиц записываются аналогично уравнениям для жидкой фазы (2.5-2.7). По аналогии требуется моделирование неизвестных корреляций скорости и температуры частиц, аналогичных рейнольдсовым напряжениям и переносу тепла в газе. Моделирование неизвестных корреляций аналогично, как и для жидкой среды. Расчет методом Эйлера проходит параллельно для двух фаз, и полученные поля суммируются. Метод Эйлера применяется к течениям с частицами, обладающими малой инерционностью (равновесные течения). Во всех других случаях в настоящее время используют метод Лагранжа, хорошо описывающий течения с малой концентрацией частиц. При высокой концентрации частиц и большой величине инерционности (неравновесные течения) добавляются модели, учитывающие соударения между частицами (коллизионные модели) [49], а также модели, учитывающие влияние частиц на жидкую среду.

Метод Лагранжа В методе Лагранжа записываются уравнения движения для одиночной частицы.

Уравнение изменения координат d xp = up. (2.73) dt Уравнение сил, которое представляет сумму сил, действующих на частицу du p = F A + F g + F vm + F Sa + F M + F P +..., (2.74) mp dt где F A - сила аэродинамического сопротивления;

F g - сила тяжести;

F vm - виртуальная массовая сила;

F Sa - сила Сэфмена;

F M - сила Магнуса;

F P - сила вследствие градиента давления.

Кроме этого могут быть добавлены члены, учитывающие электрические, магнитные и другие силы.

Уравнение вращения d p = Ti, (2.75) Ip dt где I p - момент инерции частицы;

Ti - сумма моментов, действующих на частицу.

Ниже приведено подробное описание сил, находящихся с правой стороны уравнения (2.74).

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

f FA = m p Cd Re p U rel. (2.76) pd 4 p Коэффициент сопротивления частицы Cd в случае несжимаемого потока является функцией от числа Рейнольдса. При малых числах Рейнольдса Re p 1 используется простая зависимость Cd = 24. При больших числах Рейнольдса Re p Re p 1000 используют зависимость, предложенную Шиллером и Нойманом (1933):

( ) 1.0 + 0.15Re0.687.

Cd = (2.77) p Re p Сила тяжести определяется по формуле Fg = m p g. (2.78) Виртуальная массовая сила представляет собой массу жидкой среды, которая увлекается при ускорении частицы. Виртуальная массовая сила пропорциональна отношению плотностей жидкой f среды и частиц = p :

d ( ) Fvm = m p Cvm u f up. (2.79) dt Коэффициент Cvm представляет собой долю объема жидкой среды по отношению к объему частиц. При низкой объемной доле дисперсной среды Cvm принимается равным 0.5.

Сила Сэфмена возникает вследствие неоднородного профиля скорости жидкой среды при обтекании частицы. Разница скоростей приводит к возникновению перепада давлений. Сила Сэфмена [44] определяется следующим образом 0. u FAS = 1.615d 2 U rel. (2.80) p y Сила Магнуса возникает вследствие вращения частицы. Сила Магнуса записывается в следующем виде U rel (U rel rel ), f d 2 cMa FM = (2.81) p rel здесь rel, U rel - скорость вращения и переносная скорость относительно потока.

Сила вследствие градиента давления записывается как f Du f f Du f FP = m p mp g. (2.82) p Dt p Dt Расчет течений методом Лагранжа проходит итеративно, в два этапа. В начале рассчитывается поле для жидкой фазы (2.5) – (2.7), потом для частиц (2.73) – (2.75), далее снова решаются уравнения для жидкой фазы с источниковым членом, вычисленным из второго этапа. Итерации продолжаются до полной сходимости обоих этапов.

2.3.9. Моделирование горения Как уже упоминалось в главе 1.3, структура диффузионного горения описывается уравнением сохранения пассивного скаляра (mixture fraction). При турбулентном горении осредненное по Фавру уравнение для переноса пассивного скаляра записывается в следующем виде:

ui D u.

+ = (2.83) j t xi xi xi Уравнение содержит неизвестные корреляции, моделирование которых рассматривалось в главе 2.3.7. Однако, в отличие от ламинарного диффузионного горения, турбулентное диффузионное горение не может представляться только значениями осредненных () () долей в виде Yk = Yk и Tk = Tk вследствие появления флуктуаций смешанных долей. Связь химической и турбулентной составляющих горения была рассмотрена в начале главы 2.3.

Значение переменной Yk может быть найдена из PDF интеграции Yk = Yk ( )P ( ) d. (2.84) Здесь для PDF функции P ( ) использована форма - функции.

Таким образом, к значению смешанной доли требуется добавить его вариации 2. Используя уравнение (2.13), транспортное уравнение для вариаций можно записать как 2 u j 2 D u 2 u + = j j x j t x j x j x j (2.85) 2.

Sc x j x j Уравнение (2.85) содержит больше неизвестных переменных, чем уравнение (2.83). Дополнительно к турбулентному потоку u появляются тройные корреляции u 2 и скорость j j диссипации скаляра. Наиболее простой способ x j x j моделирования первых двух членов – использование градиентного предположения, аналогично, как и для потока скаляра при моделировании процессов теплообмена и смешения (глава 2.3.7.). В этом случае перенос вариаций записывается как t u =. (2.86) j xi Тройные корреляции моделируются аналогично:

t u =. (2.87) j xi Член диссипации моделируется по следующей формуле [35]:

= 2 2. (2.88) Sc x j x j k Более сложные модели рассмотрены в [35], [47], [59].

Моделирование турбулентного предварительно-смешанного горения представляет собой более сложную задачу, чем моделирование диффузионного горения. В FASTEST-3D имеются следующие модели: BML и G–equation. В BML модели решается скалярное уравнение для прогрессивной переменной, локализующей фронт пламени (progress variable), которая может представлять собой функцию от плотности, температуры или концентрации. В G–equation модели решается скалярное уравнение для фронта пламени. Более подробно о моделях для кинетического горения изложено в работах [6], [40], [63]. При частично смешанном горении используются комбинированные модели на базе моделей для диффузионного и кинетического горения.

Ниже в качестве примера приведены осредненные по Фавру уравнения к модели с BML моделью горения [35]:

t k k k u ( ui k ) = uiu j i + + + k x j t xi x j x j (2.89) p ui c, u xi 2 t ui + ( ui ) = c1 ui u j c 2 + + k xj xj t xi xj k (2.90) p c 3 uic.

xj k Последние члены уравнений (2.88) и (2.89) являются добавочными членами BML модели. Дополнительная константа c 3 принимается обычно равной 1.44. Модель рейнольдсовых напряжений для горения (BML), а также уравнение для прогрессивной переменной приведены в приложении. Ряд других моделей горения рассмотрены в [35], [47], [59].

2.3.10. Моделирование пристеночных течений Пристеночные функции Приведенные в предыдущих главах турбулентные модели пригодны для полностью развитых турбулентных течений. Однако вблизи стенок местное турбулентное число Рейнольдса Ret является столь малым, что вязкие эффекты превалируют над турбулентными. Один из наиболее распространенных подходов к моделированию пристеночных течений связан с использованием метода пристеночных функций, который обладает двумя очевидными достоинствами: позволяет экономить вычислительные ресурсы и учитывать влияние различных факторов, в частности, шероховатости за счет введения эмпирических поправок.

Рис. 2.2. Турбулентный профиль скорости в пристеночной зоне Пристеночная область течения может быть разбита на три зоны:

1) вязкий подслой, в котором вязкие напряжения доминируют над рейнольдсовыми и имеет место линейная зависимость скорости uy u потока от расстояния до стенки: u + = y +, где u + =, y + =, а u u = w - динамическая скорость, 2) буферный слой, где вязкие и рейнольдсовы напряжения имеют один порядок. Соединяя профили для вязкого подслоя и логарифмического слоя, приближенно получают u + = 5 ln y + + 3.05. Часто буферным слоем пренебрегают, считая его частью вязкого подслоя. Объединенная зона лежит в диапазоне 0 y + 11.63 (рис. 2.2), 3) в логарифмическом слое y + 11.63 рейнольдсовы напряжения значительно превышают вязкие, а профиль скорости может быть представлен в форме логарифмического закона:

() u + = ln Ey + + B, (2.91) к где к 0.41- постоянная Кармана;

Е - постоянная, определяющая степень шероховатости (для гладкой стены E = 8.8 );

B 5.0 безразмерная константа.

Описанные участки обычно объединяются в одну внутреннюю область, которая занимает порядка 20% толщины турбулентного пограничного слоя и в которой генерируется около 80% всей энергии турбулентности. Одно из важных свойств внутренней области заключается в том, что профиль скорости слабо зависит от числа Рейнольдса, продольного градиента и прочих внешних условий. Именно это свойство послужило основой для построения универсальных соотношений (пристеночных функций), связывающих параметры течения с расстоянием от стенки. Наряду с универсальностью профиля скорости во внутренней области, метод пристеночных функций опирается на использование гипотезы о локальном равновесии турбулентных пульсаций, а также свойства локальной изотропности диссипирующих вихрей.

Подробное описание моделирования пограничного слоя приведено в работе Шлихтинга [45].

Аналогично строятся пристеночные функции для пассивного скаляра (Кадер, Яглом [29]) для описания процессов переноса тепла:

1) для вязкого подслоя y + 11:

t + = Pr y +, (2.92) где - эффективный перенос;

t - турбулентный перенос;

+ нормированный пассивный скаляр, 2) буферный слой 0 y + 11.63 :

t, (2.93) 3) логарифмический слой y + 11.63 :

k t + = ln y + + D. (2.94) kt В этой области скаляр + получается из наложения скоростного и скалярного поля:

( ) + = Prt U + + f, (2.95) где функция f найдена эмпирически согласно Кадеру и Яглому [29] f = 12.5Pr 2 / 3 + 2.12ln Pr + 1.5. (2.96) Эффективный перенос описывается следующим образом:

y+ eff =. (2.97) T+ Пристеночные эффекты Существует множество течений, где пристеночные функции дают неточные результаты. Одним из таких распространенных течений является течение с отрывом. Для корректировки этого предлагается в зоне пограничного слоя использовать различные демпфирующие функции, зависящие от пристеночных чисел Рейнольдса, а также включение члена, описывающего молекулярный перенос. Введение дополнительных членов, например, в уравнения для k по сравнению с исходными уравнениями (2.29), (2.30) обусловлено как потребностью более точного определения искомых функций в области малых значений Ret в непосредственной близости от стенки, так и тем, что скорость диссипации s принимает ненулевое значение на стенке, в то время как кинетическая энергия на стенке равна нулю. Это значит, что соотношение 2 / k в уравнении (2.30) стремится к бесконечности. Для устранения этого недостатка вводится понятие гомогенной диссипации = 0 +, где – гомогенная диссипация, которая находится из уравнения диссипации;

0 - диссипация на стенке (высчитывается по-разному в зависимости от модели) и – суммарная диссипация.

Ниже представлены низкорейнольдсовые (Low- Re) варианты некоторых рассмотренных ранее моделей.

Низкорейнольдсовая k модель Чена [16] k t k k u i + u j = ij + + ;

(2.98) k t x j x j x j x j u i + u j = c 1 f1 ij c 2 f 2 +E+ t x j k x j k (2.99) t + + ;

x j x j k = 0 +, 0 = 2 2 ;

(2.100) y + e y / E = 2 ;

(2.101) y k t = c f. (2.102) (Ret / 6) 0.0115 y + ;

f1 = 1;

f 2 = 1 0.22e Здесь f = 1 e демпфирующие функции.

Модельные константы:

c 1 = 1.35;

c 2 = 1.80;

c = 0.09;

k = 1.0;

= 1.3.

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

Ниже представлены низкорейнольдсовые варианты для моделей WJ и CLS.

Явная низкорейнольдсовая алгебраическая модель рейнольдсовых напряжений (Валлин & Йоханссон) [60] ( ) max(BII, II eq ) 2 S 2 1 II s I + f12 3 2 1 II I + 32 aij = 1 f ex 3 s s 2 ( ) B ( S S ) + + f1 4 1 f1 eq 2 max( II s, II s ) ( ) + f1 6 S 2 2 S IVI + f12 9 S 2 2 S. (2.103) Для демпфирующей функции используется функция Ван Дриста y+ f1 = 1 exp +, (2.104) A где A+ = 26.

Максимум-функция включена в уравнение (2.86). Чтобы избежать проблемы низкого значения II s в отрывных течениях, eq максимум-функция лимитируется значением II s для равновесных течений:

405c eq = 5.74, (2.105) II s 216c1 где c1 = 1.8.

Ниже приведена нелинейная низкорейнольдсовая модель (Craft, Launder & Suga) [18]:

t t 2 1 S II s I + c2 t ( S S ) + a= S + c k k ( ) + c3 t 2 II I + c4 t2 S 2 S 2 + (2.106) k + c5 t2 S 2 2 S IVI + k k + c6 t2 Sij S jk Ski + c7 t2 Sij jk ki ;

k t = c f ;

(2.107) 0. 0. c = 1 exp. (2.108) exp ( 0.75max ( S, ) ) 1. 1 + 0.35 ( max ( S, ) ) Демпфирующая функция f = 1 exp ( Ret / 90 ) ( Ret / 400 ).

1/ 2 (2.109) Местное число Рейнольдса Ret = k 2 /. (2.110) Модельные константы:

2 2 c1 = 0.1, c2 = 0.1, c3 = 0.26, c4 = 10c, c5 = 0, c6 = 5c, c7 = 5c.

Аналогично вводятся демпфирующие функции для уравнения переноса скаляра. Ниже приведена демпфирующая функция, используемая в модели Лэм-Брэмхорст (LB) [8] f T = (1 e0.255 Ret ) 2 (1 + ). (2.111) Ret К сожалению, ввод демпфирующих экспоненциальных функций создает определенные вычислительные трудности;

кроме этого, для низкорейнольдсовых моделей требуется очень мелкая сетка около стенок, что значительно замедляет скорость вычислений [64] (табл. П.2). Низкорейнольдсовые варианты существуют также для моделей рейнольдсовых напряжений;

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

Глава 3. ЧИСЛЕННЫЕ МЕТОДЫ К сожалению, решить систему дифференциальных уравнений, представленную в предыдущей главе, аналитически не представляется возможным. Поэтому решение осуществляется численно. Численный метод включает в себя: математическую модель, метод дискретизации, численную сетку, конечную аппроксимацию и метод решения. FASTEST–3D (Flow Analysis Solving Transport Equations Simulating Turbulence 3 Dimensional) имеет следующие особенности [20]:

- метод дискретизации конечных объемов, основанный на гексаэдрных контрольных объемах;

- декартовая система координат;

- неявная и полунеявная временная схема по времени;

- метод дискретизации в пространстве первого и второго порядка точности;

- строго неявная процедура (Strongly Implicit Procedure) для итерационного решения линеаризированной системы уравнений;

- параллелизация, основанная на пространственном разбиении (domain decomposition).

3.1. Метод конечных объемов В общем виде дифференциальные уравнения переноса величины можно представить как ui + = S, (3.1) Г t xi xi xi IV I II III где первый член описывает нестационарность;

второй- конвекцию;

третий- диффузию;

справа находится источниковый член.

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

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

( ) ( ) dV + t ui dV dV = S dV, (3.2) Г xi xi xi V V V V здесь - величина (скорость, скаляр);

Г - представляет собой коэффициент переноса (ламинарная/турбулентная вязкость, коэффициент диффузии и т. д.);

S - силы, не включенные в конвективную и диффузионную часть, а также дополнительные силы (сила тяжести, радиация и т.д.).

Используя теорему Остроградского–Гаусса, интеграл по объему можно записать через поверхностный интеграл:

divdV = ndF = S dV. (3.3) V Получаем ( ) dV + ui Г ni d = S dV, t (3.4) xi V V где – поверхность, ограничивающая объем V и ni, представляет собой компоненту вектора n нормального к поверхности и направленной наружу. Вышеизложенное уравнение должно решаться в каждом контрольном объеме.

В FASTEST-3D в качестве контрольного объема используется гексаэдр (рис. 3.1).

Рис. 3.1. Неортогональный гексаэдрный конечный объем Каждый гексаэдрный конечный объем содержит угловые точки и центральную точку. В центральной точке сохраняется значение величин для контрольного объема.

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

3.2. Преобразование координат С учетом неортогональности используемой численной сетки более удобным является использование в каждом контрольном объеме локальной координатной системы с трансформацией векторных величин из локальных координат в глобальные. На рис.

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

Локальные координаты обозначаются так (1, 2, 3 ), глобальные как ( x1, x2, x3 ). В [21] изложены преимущества использования декартовой систему координат в качестве базовой системы.

Рис. 3.2. Локальная система координат в контрольном объеме Матрица трансформирования [21] из глобальных координат в локальные записывается следующим способом x1 x x 1 2 x x x T = 2. (3.5) 1 2 x3 x x 1 2 Переменная по отношению к декартовой координатной системе может быть записана в локальных координатах как j =. (3.6) xi j xi Элементы обратной транспортной матрицы T 1 (локальные в глобальные координаты) записываются следующим образом:

1 x1 x2 x 1 2 T =, (3.7) x1 x2 x 3 x1 x2 x которая найдена из известного линейного алгебраического отношения (Tad )T, T= (3.8) J (Tad )T - дополнительная матрица, где J = det T - якобиан и найденная из матрицы T.

Таким образом, j1 x = ij = adj i (3.9) xi J i J и 1 = ij, (3.10) j xi J где ij - элементы матрицы, записываемой как x2 x3 x2 x3 x1 x2 x1 x x1 x3 x1 x 2 3 3 2 3 2 2 3 2 3 3 x x x x x1 x2 x1 x x1 x3 x1 x = 2 3 2 3.

3 1 1 3 1 3 3 1 3 1 1 x2 x3 x2 x3 x1 x2 x1 x x1 x3 x1 x 1 2 2 1 2 1 1 2 1 2 2 (3.11) Окончательно, подставляя полученное отношение для дифференциального оператора (3.10) в уравнение (3.2), получаем Г ( ) dV + ui ni d t ik ni d = S dV. (3.12) J k V V Используя дискретизационную транспортную матрицу и форму Гаусса (3.12), можно записать согласно [21] как ( ) Г t l bk = S dV, dV + ui Fnb bkj ui j l (3.13) V xi V V j где оператор Fnb - значение поверхностного интеграла на каждой стороне;

bkj - коэффициент переноса. Интеграция проходит по каждой стороне, в направлении j = 1 для “East” (индекс ’E’) и “West” (индекс ’W’);

в направлении j = 2 для “North” (индекс ’N’) и “South” (индекс ’S’);

в направлении j = 3 для “Top” (индекс ’T’) и “Bottom” (индекс ’B’) (рис. 3.3).

Рис. 3.3. Соседние стороны в одном контрольном объеме 3.3. Дискретизация диффузионного члена Для простоты рассмотрим только одну сторону контрольного объема – сторону ‘East’, изображенную на рис. 3.4. Аналогичная операция происходит и для других сторон. Соседние центральные точки обозначены буквами ‘P’ и ‘E’;

точки, лежащие на ребрах грани, обозначены буквами ‘te’, ‘ne’, ‘se’ и ‘be’.

Рис. 3.4. Сторона ‘East’ контрольного объема Поверхностный интеграл для диффузионного члена через сторону ‘East’ записывается следующим образом [21]:

() () () Г Г i bk = b1' 2 2 Fe1{bk 1 i 1 + b2 + b3 1 + v v неявная часть (3.14) Г 1 1 1 2 Г 1 3 1 3 1 b b + b b + b b, b1 b2 + b2b2 + b3b 3 2 + + v 1 1 2 2 3 3 v явная часть где разности i определяется через центральные точки и граничную поверхность:

1 = E P, 2 = ne se, (3.15) 3 = te be.

Неявная схема дискретизации всех членов для диффузии требует большую матрицу коэффициентов, поэтому все члены, содержащие E и P, находятся явно и добавляются в источниковый член (3.14). Значения ne, se, te, be находятся из трилинейной интерполяции между точками P, E, N, S, T, B, ne, se, te, be.

3.4. Дискретизация конвективного члена Аналогично диффузионному члену конвективный поток через сторону ‘East’ можно записать в виде ( )e e = m e e, i ( u )} = Fe1{bk 1 1 1 e b1 u1 + b2u2 + b3u3 (3.16) i где me обозначает конвективный перенос через грань ‘East’, ( ) i 1 1 me = b1 u1 + b2u2 + b3u3. (3.17) Из уравнения неразрывности это записывается как i i i i i i me + mn + mt m w m s mb = 0. (3.18) В FASTEST-3D значения скоростей и давления сохраняются в центральной точке контрольного объема;

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

3.5. Дискретизация временного члена В нестационарных течениях (или при решении псевдо нестационарным методом) в каждом транспортном уравнении (3.1) присутствует нестационарный член.

Нестационарный член аппроксимируется следующим образом:

dV ( ) P VP.

(3.19) t t VP Для аппроксимации по времени можно использовать явные методы. Одно из простейших предположений – явный метод Эйлера, который имеет первый порядок точности:

( P V )n ( P V )n1 = ANb1 n1 + S P1 AP1 n1. (3.20) n n n Nb P t Nb Здесь t обозначает временной шаг, а индексы n и n 1 относятся к последнему и предыдущему шагу. Первый член с левой стороны (3.20) определяется явно и добавляется в коэффициенты AP, второй же член добавляется в источниковый член S p. В общем случае дискретизация временного члена записывается в следующем виде:



Pages:   || 2 |
 




 
2013 www.netess.ru - «Бесплатная библиотека авторефератов кандидатских и докторских диссертаций»

Материалы этого сайта размещены для ознакомления, все права принадлежат их авторам.
Если Вы не согласны с тем, что Ваш материал размещён на этом сайте, пожалуйста, напишите нам, мы в течении 1-2 рабочих дней удалим его.