Л. В. Быков, А. М. Молчанов, Д. С. Янышев
основы
ВЫЧИСЛИТЕЛЬНОГО
ТЕПЛООБМЕНА И
ГИДРОДИНАМИКИ
URSS
Л. В. Быков, А. М. Молчанов, Д. С. Янышев
ОСНОВЫ
ВЫЧИСЛИТЕЛЬНОГО
ТЕПЛООБМЕНА
И ГИДРОДИНАМИКИ
Допущено
федеральным учебно-методическим объединением
в системе высшего образования по укрупненной группе
специальностей и направлений подготовки
«Авиационная и ракетно-космическая техника»
в качестве учебного пособия для студентов, обучающихся по основным
образовательным программам высшего образования по направлению
подготовки бакалавриата/специалитета/магистратуры 24.00.00
«Авиационная и ракетно-космическая техника»
Издание второе,
переработанное и дополненное
URSS
МОСКВА
ББК 22.317 22.353.3
Быков Леонид Владимирович,
Молчанов Александр Михайлович,
Янышев Дмитрий Сергеевич
Основы вычислительного теплообмена и гидродинамики:
Учебное пособие. Изд. 2-е, перераб. и доп. — М.: ЛЕНАНД, 2019. — 200 с.
Целью данного пособия является ознакомление читателей с основами моде-
лирования процессов теплопроводности и конвективного теплообмена. В пособии
излагаются основы метода конечных объемов для решения задач теплообмена
и гидродинамики, а также метода конечных элементов для решения задач тепло-
проводности. Приводится информация о моделировании турбулентности и методах
расчета химически реагирующих течении.
Рекомендуется студентам технических вузов, аспирантам и научным работ-
никам, специализирующимся в области современных проблем тепломассообмена.
Рецензенты:
главный конструктор по ВСО АО «ГН1Ш „Регион"»,
д-р техн, наук И. В. Гаранин;
проф. Института кибернетики НИ ТПУ, д-р физ.-мат. наук Г. Я. Мамонтов;
кафедра «Теоретическая и промышленная теплотехника» Национального
исследовательского Томского политехнического университета ( ГПУ)
ООО «ЛЕНАНД». 117312, Москва, пр-т Шестидесятилетия Октября, д. НА, стр. 11.
Формат 60x90/16. Печ. л. 12,5. Зак. № 131461
Отпечатано в АО «Т 8 Издательские Технологии».
109316, Москва, Волгоградский проспект, д. 42, корп. 5.
ISBN 978-5-9710-5803-8
11956 ID 238304
©ЛЕНАНД, 2018
НАУЧНАЯ И УЧЕБНАЯ ЛИТЕРАТУРА		
	1 URSS	E-mail: URSS@URSS.ru Каталог изданий в Интернете: httpy/URSS.ru Тел/факс (многоканальный): i	+7(499)724 25 45
Все права защищены. Никакая часть настоящей книги не может быть воспроизведена или
передана в какой бы то ни было форме и какими бы то ни было средствами, будь то элек-
тронные или механические, включая фотокопирование и запись на магнитный носитель,
а также размещение в Интернете, если на то нет письменного разрешения владельца.
Оглавление
Введение........................................6
Глава 1. Необходимое математическое введение....9
1.1.	Основные сведения
из векторного анализа.................11
1.2.	Основы тензорной нотации.........17
1.3.	Основные уравнения
механики жидкости и газа..............23
Глава 2. Основные методы численного решения
задач тепломассообмена........................37
Глава 3. Основы метода конечных объемов........43
3.1.	Простой пример...................45
3.2.	Уравнение энергии................49
3.3.	Дискретные аналоги
поверхностных интегралов..............52
3.4.	Дискретные аналоги
объемных интегралов...................60
3.5.	Граничные условия................61
3
Оглавление
3.6.	Производная по времени.............67
3.7.	Общий алгоритм решения задачи......74
Глава 4. Расчет поля течения.....................77
4.1.	Интегрирование уравнения движения...80
4.2.	Расчет поля давления...............82
4.3.	Подходы и проблемы.................88
Глава 5. Турбулентность: проблемы
моделирования и подходы к их решению..............91
5.1.	Феномен турбулентности.............94
5.2.	Явления отрыва.....................98
5.3.	Осредненное движение.
Уравнения Рейнольдса..................101
5.4.	Гипотеза Буссинеска...............106
5.5.	Путь смешения Л. Прандтля.
Алгебраические модели ..................109
5.6.	Модель к-Е.
Дифференциальные модели...................111
5.7.	О моделировании течений
вблизи стенки...........................114
5.8.	Общие данные о некоторых
моделях турбулентности....................119
5.9.	Более сложные модели турбулентности.122
5.10.	О применении моделей турбулентности.... 127
4
Оглавление
Глава 6. Особенвости расчёта химически
реагирующих течений..................129
6.1.	Основные положения..............132
6.2.	Основные уравнения..............135
6.3.	Жёсткие системы.................148
6.4.	Решение жёстких систем уравнений.160
6.5.	Метод расщепления для системы
уравнений переноса химических
компонентов.........................164
Глава 7. Особенности расчёта
многофазных течений.....................>....167
7.1.	Подход «Эйлер—Эйлер».............175
7.2.	Подход «Эйлер—Лагранж»...........182
7.3.	Массообмен между частицами
и непрерывной средой.................185
7.4.	Численный метод решения уравнений
переноса дискретной фазы.............188
7.5.	Общие рекомендации по численному
моделированию многофазных течений....191
Заключение....................................193
Список использованной литературы..............194
Введение
В настоящее время с интенсивным развитием компью-
терных технологий особое значение приобретает матема-
тическое моделирование различных физических процес-
сов. В задачах тепло- и массообмена численный экспери-
мент приобрел сейчас значение, сравнимое со значением
натурного эксперимента.
Целью данного пособия является ознакомление читате-
лей с основами моделирования процессов теплопроводно-
сти и конвективного теплообмена. Часть материала, пред-
ставленного в пособии, носит описательный характер и
служит для создания представления о существующих на
сегодняшний день математических моделях. Такой подход
к изложению обусловлен тем, что в настоящее время су-
ществует целый ряд коммерческих пакетов программ по
вычислительной гидродинамике и теплопередаче. Как пра-
вило, все они рассчитаны на продвинутых пользователей,
владеющих основными моделями, использование которых
требует знания основных подходов к моделированию те-
чений жидкостей и газов и применимости этих подходов
6
Введение
при решении конкретных задач с учётом имеющихся вы-
числительных ресурсов.
В связи с этим в данном пособии особое внимание уде-
лено сопоставлению русскоязычных и англоязычных тер-
минов и названий, поскольку большинство существующих
пакетов программ по вычислительной теплопередаче и
гидродинамике выпускаются исключительно на англий-
ском языке и не имеют локализации.
Для понимания сути математических выкладок от чита-
теля требуется знание некоторых основ векторного анали-
за. Основные сведения оттуда представлены ниже.
9-

Глава 1
Необходимое
математическое
введение

1.1. Основные сведения
из векторного анализа
Скалярное и векторное произведение векторов
Скалярное произведение двух векторов есть скаляр.
Пусть имеется вектор а с координатами ах, ay, az и век-
тор b с координатами bx, by, bz. Здесь и далее используется
декартова система координат. Тогда скалярное произведе-
ние этих векторов будет равно
а • b = axbx + ayby + azbz = |а| • |b| cos(a b).	(11)
Векторное произведение двух векторов является векто-
ром и вычисляется следующим образом:
а х b = (aybz - azby )i + (azbx - axbz) j +
+ (аД-«Д)к= Ч
(1-2)
Модуль векторного произведения равен произведению
модулей векторов на синус угла между ними
л
|а х b| = |а| • |b| sin(ab).
11
Глава 1. Необходимое математическое введение
Оператор набла (оператор Гамильтона)
Оператор набла (оператор Гамильтона) записывается
следующим образом:
V7 . д . д ,
V = —1 + — 1 +—к,
дх ду dz
(1.3)
где i, j, к — единичные векторы.
Если скалярно умножить оператор набла на векторную
величину, то мы получим дивергенцию данного вектора:
Va = diva =
Эа_ 5а даг
—- + —- + —
дх ду dz
(1.4)
ах, ау и аг — проекции вектора а на соответствующие оси
координат.
С точки зрения физики, дивергенция векторного по-
ля является показателем того, в какой степени данная
точка пространства является источником или стоком это-
го поля:
•	div(a) >0 — точка является источником поля а.
•	div(a) <0 — точка является стоком поля а.
•	div(a) = 0 — стоков и источников нет, либо они
компенсируют друг друга.
Если умножить оператор набла на скаляр, получается
градиент этого скаляра:
12
1.1. Основные сведения из векторного анализа
V^ = gradp = ^i+^j+^k.
дх ду dz
(1.5)
В различных отраслях физики используется понятие
градиента различных физических полей.
Например, градиент концентрации— нарастание или
уменьшение по какому-либо направлению концентрации
растворенного вещества, градиент температуры— уве-
личение или уменьшение по направлению температуры
среды и т. д.
Третьим важнейшим оператором векторного анализа
является ротор или вихрь. Как уже видно из названия, он
характеризует вихревую составляющую векторного поля,
показывая, насколько закручено поле в данной точке.
В русскоязычной литературе оператор ротора обозначается
как rot(a), в англоязычной — как curl(a).
Ротор есть вектор со следующими координатами (де-
картова система координат):
rota =
ду dz >
V dz дх
(1.6)
(дау да
__У х к
I Sr ду )
Для удобства запоминания можно представить ротор в
виде векторного произведения оператора набла и вектора
поля:
13
Глава 1. Необходимое математическое введение
	i	j	к
	д	д	д
rota = Vxa =			
	дх	ду	dz
	а X	а У	
(1-7)
(1.8)
Оператор Лапласа
Оператор Лапласа или лапласиан записывается сле-
дующим образом:
. д2 д2	д2
Л == "	"	  
дх2	ду2	dz2 ’
Можно показать, что оператор Лапласа есть скалярное
умножение оператора набла на самого себя, или, что есть
то же самое, дивергенция градиента:
Да = VVa = V2a = div(grada).	(1.9)
Наряду с лапласианом скалярной функции, существует
так же и лапласиан вектора. Он обозначается точно так же
и тоже представляет собой сумму вторых производных.
При этом, однако, лапласиан векторной функции является
вектором, а не скаляром.
Да = V2a = grad(diva) - rot(rota) =
_ Э2а д2а 52а
дх2 ду2 dz2
(110)
14
1.1. Основные сведения из векторного анализа
Как уже упоминалось ранее, при записи выражений
использована декартова система координат. Однако сле-
дует подчеркнуть, что при необходимости любое вектор-
ное выражение можно записать в произвольной системе
координат (в сферической, цилиндрической и т.д. — см.,
например, [12]).
Некоторые часто встречающиеся соотношения:
div($9 • a) = q> • diva + a • grad $9 =
= ^Va + a-V^,
дя дя дя
a- Via = aY — + a„— + a7 — =
7 x дх ydy z dz
, а	•
= grad — +rotaxa,
k 2 J
rot (grad$9) = 0,
(1.11)
Поток вектора
Поток вектора а через элементарную площадку dS оп-
ределяется как а • ndS или а • dS (dS = n • dS}, где n — еди-
ничный вектор нормали к площадке.
<D = Jad8.	(1.12)
s
15
Глава 1. Необходимое математическое введение
Теорема Остроградского—Гаусса
Рассмотрим векторное поле а, проходящее через объем
V, ограниченный поверхностью S. Тогда поток вектора а
через поверхность S будет равен интегралу дивергенции
этого вектора по объему V.
^adS = |divaJK.	(1.13)
S	V
Физически это можно интерпретировать следующим
образом: поток векторного поля через замкнутую поверх-
ность зависит от наличия в объеме, ограничивающем эту
поверхность источников или стоков рассматриваемого по-
ля. Если источников и стоков в данном объеме нет, или
они компенсируют друг друга, то поток вектора через
замкнутую поверхность равен нулю, т.е. сколько в объем
«втекает», столько из него и «вытекает».
Производные от формулы Остроградского—Гаусса
Можно показать (см. например [7]), что через поверхно-
стный интеграл можно выразить так же и объемные инте-
гралы от градиента и ротора. Приведем здесь эти формулы:
jgradpJK —	- dS,
V	S
JrotarfK = ^dSxa.
V	s
(1-14)
(1-15)
1.2. Основы тензорной нотации
В настоящее время в литературе очень широко исполь-
зуется тензорная запись уравнений механики жидкости и
газа (или механики сплошных сред). С целью ознакомле-
ния читателя с таким видом записи здесь приводятся ос-
новные ее правила. При этом теория тензорного исчисле-
ния здесь излагаться не будет. Кроме того, для простоты
восприятия мы ограничимся рассмотрением тензоров в де-
картовой системе координат, не затрагивая вопросов кри-
волинейных координат. Желающие могут подробно озна-
комиться с материалом по данным вопросам в специаль-
ной литературе, например, в [3,6].
Одной из причин применения тензорной формы записи
является то, что она более короткая. Достаточно часто это
помогает лучше уловить физический смысл математиче-
ского выражения.
Для начала рассмотрим самый простой пример. Возь-
мем два вектора а и Ь. Аналитически эти вектора могут
быть записаны следующим образом:
17
Глава 1. Необходимое математическое введение
ft	^2®2	^3®3 ’
b = bfti + ^2®2 + ^ЗеЗ ’
где ei, ei, ез — единичные координатные векторы.
Теперь, если мы захотим аналитически вычислить век-
тор с, который представляет собой сумму этих векторов, то
нам нужно будет записать:
с = (а, + Ьх) в) + (а2 + Ь2) е2 + (а3 + Ь3) е3.	(1.16)
Однако эта запись довольно длинная. Ее можно упро-
стить, если мы введем индекс I, который последовательно
пробегает значения от 1 до 3. Тогда вектора а и b можно
будет записать как а, и bi. А вектор с запишется как:
с^а^Ь,.	(1.17)
Точно так же мы можем записать сумму двух квадрат-
ных матриц. Допустим, мы хотим записать матрицу С,
элементы которой есть суммы соответствующих элемен-
тов матриц Л и В. Для этого воспользуемся двумя индек-
сами — z иу, которые так же пробегают значения от 1 до 3.
Тогда
С,=Л+г«'	(118)
Диапазон значений, которые пробегают индексы обыч-
но оговаривается в самом начале статьи (монографии,
18
1.2. Основы тензорной нотации
учебника) и более не оговаривается. Для трехмерных задач
это всегда будет от 1 до 3, для четырехмерных (задачи тео-
рии относительности) — от 1 до 4 и т.д. Для обозначения
индексов используют латинские буквы из середины алфа-
вита (/, j, к, I, т, р, q, г), однако обычно не принято прибе-
гать к буквам из начала (а, Ь, с...) и конца алфавита (и, v, w,
х, у, z). Так же в качестве индексов могут быть использова-
ны греческие буквы а, р, у.
Одним из самых главных в тензорной записи является
свертка или условие о суммировании, предложенное
Эйнштейном.
Правило Эйнштейна
Если в одночленном выражении имеются два одинако-
вых индекса, т. е. индекс повторяется, то этот индекс назы-
вается немым. Наличие немого индекса означает суммиро-
вание по всему диапазону, по которому пробегает индекс,
т. е. в нашем случае от 1 до 3. Результат этой операции на-
зывается сверткой; часто саму эту операцию называют
также сверткой, иногда свертыванием [3].
аи ~ аи + Л22 + азз-	(1-19)
То же самое можно было бы написать, используя знак
3
суммы:	.
i=i
19
Глава 1. Необходимое математическое введение
(1.20)
(1-21)
Используя свертку, можно так же записать операцию
скалярного умножения векторов:
afij = аД + а2Ь2 + а3Ь3.
Или операцию дивергенции вектора:
да, да, да, да,
дх, Sxj дх2 дх3
Следует также заметить, что повторение индекса более
одного раза не допускается правилом. Т.е. запись вида а»л
не имеет смысла.
Символ Кронекера
Этот символ очень часто используется в тензорном ис-
числении. Он записывается следующим образом:

1, при1 = у
0, при!*/
(1-22)
Фактически он представляет собой единичную матри-
цу, размерность которой зависит от того диапазона, кото-
рый пробегают индексы i и j.
Символ Леви-Чивиты
Символ Леви-Чивиты, или абсолютно антисиммет-
ричный единичный тензор записывается следующим
образом:
20
1.2. Основы тензорной нотации
1, если элемент имеет чётную
перестановку индексов
eijk = < -I, если элемент имеет нечётную	(1.23)
перестановку индексов
О, при i = j\i = kJ = к
Чтобы лучше понять, что значит четные и нечетные пе-
рестановки, рассмотрим всевозможные перестановки ком-
бинации чисел 1,2,3. Всего таких перестановок возможно
3!=6.
1,2,3 — четная перестановка;
1,3,2 — нечетная перестановка;
3,1,2 — четная перестановка;
3,2,1 — нечетная перестановка;
2,3,1 — четная перестановка;
2,1,3 — нечетная перестановка.
Для наглядности распишем символ Леви-Чивиты:
ГО О ОЛ <0 О
-Г|
О
(1-24)
21
Глава 1. Необходимое математическое введение
Основные формулы векторной алгебры и анализа в тен-
зорных обозначениях
Векторные обозначения	Тензорные обозначения
1). Вектор	
а	
2). Модуль вектора	
|а|	у1а{а{
3). Скалярное произведение	
а-Ь	aibi
4). Векторное произведение	
axb	eiklak^l
5). Градиент	
grad^	дф дхр
6). Дивергенция	
diva	daq dxq
7). Ротор	
rota	дак eijk я J dxj
8). Лапласиан	
Чгф = Ьф	д дф axjsxj
1.3.	Основные уравнения механики
жидкости и газа
Здесь мы кратко рассмотрим вид уравнений гидроди-
намики и теплообмена, а также приведем их вывод.
Уравнение неразрывности (continuity equation)
Его также называют уравнением сохранения массы.
Оно выводится из очень простых предположений.
Рассмотрим объем V, ограниченный некоторой произ-
вольной поверхностью S (см. Рис. 1.1)
Рис. 1.1. К выводу уравнения неразрывности.
23
Глава 1. Необходимое математическое введение
Каким образом может измениться масса данного объе-
ма? Очевидно, что только за счет втекающего и вытекаю-
щего из него потока жидкости. Тогда скорость изменения
массы в объеме будет равна потоку (расходу) жидкости
через данный объем.
Известно, что расход жидкости с плотностью р, проте-
кающей со скоростью U через сечение F вычисляется по
формуле:
G = pUF.	(1.25)
Тогда поток через малую площадку dS будет равен
pV dS. Интегрируя по всей поверхности S, получим:
-^y = fpUdS.	(1.26)
S
Учитывая, что т = J pdV, получим:
V
-f^</r = 4>pUdS.	(1.27)
V С)Т	S
Теперь, используя теорему Остроградского—Гаусса
(формула (1.13)), получим:
t^-dV = [div(pU)</K,
(1-28)
+ div(pU) dV = 0.
24
1.3. Основные уравнения механики жидкости и газа
Интеграл по произвольному объему от непрерывной
функции может быть равен нулю, только если сама функ-
ция тождественно равна нулю.
Отсюда получаем уравнение неразрывности:
^ + div(pU) = 0.	(1.29)
В тензорных обозначениях оно запишется следующим
образом:
^+/-(рЦ) = 0-	(1.30)
дт oxj	'	'
В несжимаемой жидкости плотность постоянна, т.е.
— = 0 и уравнение (1.29) можно будет записать в виде:
дт
divU = 0 или = 0.	(1.31)
Уравнение движения (momentum equation)
Уравнение движение сплошной среды фактически
представляет собой второй закон Ньютона записанный для
элементарного объема.
Пусть элементарный объем представляет собой парал-
лелепипед массы dm с гранями dx, dy и dz. Запишем для
него второй закон Ньютона:
25
Глава 1. Необходимое математическое введение
. </U JTZrfU „
dm = pdV-------= R,
dr dr
(1.32)
где R — вектор результирующей силы, действующей на
объем.
Теперь перед нами встает вопрос — а какие именно си-
лы действуют на элементарный объем жидкости? Начнем с
того, что все силы можно сразу разделить на объемные и
поверхностные. К объемным силам относят силы тяготе-
ния, электрические силы, магнитные силы и т.п. Не вдава-
ясь в природу объемной силы, будем обозначать ее J.
Поверхностные силы можно представить в виде тензо-
ра напряжений'.
(1.33)
Здесь напряжения, с одинаковыми индексами	)
суть нормальные напряжения (т.е. они направлены по нор-
мали к поверхности). Все остальные напряжения — каса-
тельные (они, соответственно, направлены по касательной).
Причем можно показать, что
Рассмотрим вывод уравнения движения в проекции на
осьх.
26
1.3. Основные уравнения механики жидкости и газа
Рис. 1.2. К выводу уравнения движения.
На Рис. 1.2 представлены все силы, действующие на
элементарный объем. Запишем выражение для равнодей-
ствующей этих сил:
Rx=JxpdV +
(1-34)
Произведя соответствующие выкладки для других осей,
расписав полную производную в выражении (1.32) как
27
Глава 1. Необходимое математическое введение
^ = ™+(U
dr дт '
аи тт аи
дт х дх
аи
аи
dz
(135)
и разделив все на dV = dx х dy х dz, получим уравнения
движения в напряжениях:
/тт Г7\ТТ	т (
+ (U-V)C7y = pJT + —— + —2-+—s-
V ’ x) и x У дх ду dz )
(1.36)
Заметим также, что при использовании тензорной но-
тации можно было бы вывести это уравнение намного бы-
стрее и записать гораздо компактнее.
Результирующая поверхностных сил Rs, действующих
на некоторый объем, есть интеграл от тензора (1.33) по по-
верхности этого объема:
Rs=\°adSr	(1.37)
s
Теорема Остроградского—Гаусса (формула (1.13))
применима и в данном случае. Преобразовав поверхност-
ный интеграл в объемный, получим:
28
1.3. Основные уравнения механики жидкости и газа
(1.38)
С учетом (1.35) и (1.38) перепишем второй закон Нью-
тона (формула (1.32)):
V
1 дх.
(1.39)
J Q X
V ил]
V
Отсюда приходим к окончательной формуле:
дт J дх.
дет
___ч_
дх.
(1.40)
Очевидно, что данное равенство полностью совпадает с
уравнениями (1.36).
Однако пользоваться уравнениями (1.36) невозможно
пока напряжения, действующие на рассматриваемый
элементарный объем неизвестны. Эти напряжения мож-
но выразить через скорости деформации среды и через
давление.
Жидкости, у которых напряжения зависят от дефор-
маций линейно, называют ньютоновскими. Касательные
напряжения в ньютоновской жидкости выражаются сле-
дующим образом (и здесь — коэффициент динамической
вязкости):
29
Глава 1. Необходимое математическое введение
	(ьих	аи;	
<Xxy=<XyX=fi	1 ду	дх	>
	( dUx dz	+—- dx >	p
	'dUz	dU + —y- dz ,	•
(1-41)
Нормальные напряжения вызывают деформацию жид-
кости не только в направлении их действия, но и перпен-
дикулярных, приводя к сдвиговым и объемным деформа-
циям. Наглядной моделью такого явления может служить
растяжение резинового стержня, уменьшающегося при
этом в диаметре [17].
Исследования показали, что нормальные напряжения
можно представить в следующей форме:
(1-42)
При использовании тензорной нотации можно было бы
записать выражение для напряжений намного короче:
30
1.3. Основные уравнения механики жидкости и газа

дх, дх, 3 дхк 7	7
J	1	К S
(1.43)
Первый член в этом выражении представляет собой
вязкие напряжения, второй — напряжения, связанные с
давлением.
Подставив выражения для напряжений в (1.36), по-
лучим:
р\ — + (и • V)U | = рЗ - VP + Д (pU) + -grad(pdiv U).	(1.44)
I дт	)	3
В тензорной нотации оно будет выглядеть как:
dut TT dut
p —L+u,—- =
П дт J dx.
. дР д
= pJt----+—
dxj
dUk dU. 2 dU,
—-+—-—o,t—-
dx, dxk 3 J dx.
(1-45)
А
Если мы рассматриваем течение несжимаемой жидко-
сти, divU = 0. Тогда уравнение движения можно записать
в виде:
р[ —+ (UV)u| = pJ-VP + A(pU)
V дт	J
или
(dUk тт диЛ	т дР а Г (диЛ
дт	7 dxj j дхк dxj dxj
(1.46)
31
Глава 1. Необходимое математическое введение
Если сложить уравнения (1.45) и (1.46) с уравнением
неразрывности в соответствующей форме ((1.30) и (1.31))
и провести соответствующие преобразования, то можно
получить уравнения движения в т.н. консервативной
форме:
дт	dxj
Т др д Г (dUk dUj 2 с диД
= pj.--+--- и —-+—-—8ik—-
dxk дх,	 дх,	дхк 3 к дх, t
*	J	\ J	* *	* У
'ги,,	=
dr dxj
(1-48)
Такая форма записи часто оказывается более удобной
при численном решении.
Уравнение энергии (energy equation)
Уравнение энергии фактически представляет собой
закон сохранения тепла, записанный для элементарного
объема.
Чтобы вывести это уравнение, будем рассуждать
следующим образом. Каким образом может измениться
32
1.3. Основные уравнения механики жидкости и газа
энтальпия h {Дж/кг} объема жидкости (напомним, что
элементарное изменение энтальпии равно dh = cpdT)? Са-
мое очевидное — из-за потока тепла за счет теплопровод-
ности. По закону Фурье поток тепла через единичную по-
верхность (плотность теплового потока) определяется как
(А, — коэффициент теплопроводности):
q = -2VT.	(1.49)
Тогда изменение энтальпии в объеме будет равно сум-
марному тепловому потоку через поверхность этого объема:
= dg[	(1.50)
у dr	s
Воспользовавшись теоремой Остроградского—Гаусса
и подставив выражение для плотности теплового потока,
получим:
£(^£1^7 .(2VT).	(1.51)
d т
п	- «	dh dh
Если жидкии объем неподвижен, то — = — , если
dr дт
dh dh тт
жидкость движется, то — = — + U • VA.
dr дт
Что же еще может изменить энтальпию жидкого объе-
ма? В жидкости может возникнуть объемное тепловыделе-
33
Глава 1. Необходимое математическое введение
ние (например, вследствие химической реакции). Обозна-
чим объемное тепловыделение через qv. Тогда уравнение
энергии запишется в виде:
+ и • V (^,7)1 - V • (AVT) + .	(1.52)
Другим важным фактором, влияющим на энтальпию
влияние на энтальпию напрямую следует из первого нача-
ла термодинамики (dQ = dh-dplр, где dQ — количество
теплоты, измеряемое в Дж/кг).
Третьим фактором является наличие работы сил тре-
ния. Часть энергии движущейся среды из-за действия этих
сил будет переходить в тепло. Часто работу сил трения
обозначают как //Ф, где Ф — диссипативная функция, ко-
торая записывается как:
dUz dUУ
34
1.3. Основные уравнения механики жидкости и газа
В тензорных обозначениях она записывается следую-
щим образом:
dU. dUk 2 диД2
<дхк dxt 3 dxt>
(1-54)
Тогда наиболее общее уравнение энергии примет вид:
^^ + U-V(pCp7’) =
= V(2VT) + gv+^ + Z^.
ат
(1.55)
Или то же самое в тензорной нотации:
д(рсрТ)_
I С/ L *	—
дт	дхк
(1.56)
Если жидкость несжимаема и скорости ее течения не-
велики, то два последних члена уравнения можно не учи-
тывать и тогда уравнение (1.56) превращается в уравне-
ние (1.52).

Глава 2
Основные методы
численного решения
задач тепломассообмена

Как известно, процессы тепло- и массообмена с точки
зрения механики сплошных сред могут быть описаны сис-
темой нелинейных дифференциальных уравнений в част-
ных производных. Эту систему уравнений принято назы-
вать уравнениями Навье—Стокса.
Для примера еще раз запишем систему уравнений для
неизотермического (с теплообменом) течения несжимае-
мой жидкости в декартовых прямоугольных координатах.
Она будет состоять из уравнений неразрывности (2.1),
движения (2.2) и энергии (2.3).
VU = 0,	(2.1)
— + (U-V)U = J-—VP+vAU,	(2.2)
дт	р
^+(UV)T = aAr+-^-.	(2.3)
5т	рср
В данных уравнениях J есть результирующий вектор
массовых сил, v — кинематическая вязкость среды (v = /z/p),
а — коэффициент температуропроводности (а = Х/рср).
При этом коэффициенты v и а считаются постоянными.
В тензорных обозначениях она будет выглядеть как:
39
Глава 2. Численное решение задач тепломассообмена
dUk тт dUk	т \ дР	д (дик}
дт	1 дх, р дхк	дх, дх,
дТ гт дТ д (дт} qv
дт dxj дх} J рср
(2.5)
(2.6)
В наиболее общем случае система уравнений Навье—
Стокса включает в себя уравнения неразрывности, движе-
ния, энергии и диффузии. Если в потоке происходят хими-
ческие реакции, задача усложняется введением уравнений
модели протекания этих реакций.
Уравнения Навье—Стокса могут быть решены в общем
виде лишь в некоторых случаях и при ряде допущений.
Общего аналитического решения системы этих уравнений
пока не получено. При этом численные методы решения
уравнений Навье—Стокса развиты довольно хорошо и на
сегодняшний день нашли широкое применение в различ-
ных областях науки и техники. Численное моделирование
уравнений Навье—Стокса является неотъемлемой частью
процесса проектирования летательных аппаратов, судов,
различных двигательных установок, объектов ракетной
техники, автомобилей и т.д.
В настоящее время развиты три основных подхода к
численному решению уравнений Навье—Стокса.
Первый из них носит название Метода конечных разно-
стей. На английском — Finite Difference Method (FDM). Его
суть заключается в прямой замене производных, входящих в
40
Глава 2. Численное решение задач тепломассообмена
исходные уравнения, их дискретными (разностными) анало-
гами. Решение ищется в узлах сетки, на которую разбивается
расчетная область. Достоинством метода является относи-
тельная простота реализации, при этом, однако с точки зре-
ния физического смысла этот метод не очень нагляден. Дру-
гим недостатком этого метода являются особые требования
к построению сетки, что часто усложняет процесс решения.
Второй называется Методом конечных объемов или
методом контрольного объема. В англоязычной литературе
он называется Finite Volumes Method (FVM). Основа мето-
да заключается в том, что расчетная область с помощью
сетки разбивается на совокупность конечных объемов. Уз-
лы, в которых ищется решение, находятся в центрах этих
объемов. Для каждого объема должны выполняться законы
сохранения массы, количества движения и энергии.
То есть, например, изменение во времени массы среды в
контрольном объеме может происходить только за счет
внешнего потока массы, входящего в объем, или за счет
потока массы из данного объема выходящего. Более под-
робно мы познакомимся с этим методом несколько позд-
нее. Метод конечных объемов применяется во многих вы-
числительных гидродинамических (CFD) пакетах, таких
как FlowVision, Flow3d, PHOENICS и ряде других.
Третий метод решения — Метод Конечных Элементов
(МКЭ). В англоязычной литературе его называют Finite
41
Глава 2. Численное решение задач тепломассообмена
Elements Method (FEM). Суть метода состоит в прибли-
женном решении вариационной задачи. Для формулировки
этой задачи напомним понятие функционала. Оператор
/Дх)] называется функционалом, заданным на некотором
множестве функций, если каждой функции Дх) ставится в
соответствие определенное числовое значение /Дх)] [1].
Иными словами, функционал является как бы «функцией
он функции». Часто функционалы имеют вид интегралов.
Вариационная задача состоит в отыскании такой функции
Дх), которой бы соответствовало минимальное значение
функционала /Дх)]. Вид этого функционала различен для
различных задач и подбирается специально.
В настоящее время Метод Конечных Элементов нашел
широкое применение при решении задач теплопроводно-
сти в твердых телах и при расчетах на прочность. Однако
он может быть применен и к задачам течения жидкостей
и газов (см. например [2]). Известны также методы, кото-
рые сочетают в себе черты метода конечных объемов и
метода конечных элементов [4, 5]. Сочетание этих мето-
дов позволяет использовать более широкий ряд расчет-
ных сеток (тетраэдрические сетки, пирамиды, призмы,
полиэдры), что необходимо при решении задач со слож-
ной геометрией.
Этот подход используют CFD пакеты Ansys CFX, Ansys
Fluent, Star-CD, Star-CCM+.
Глава 3
Основы метода
конечных объемов

3.1. Простой пример
В качестве относительно простого примера рассмотрим
основные принципы метода на уравнении неразрывности
для сжимаемой жидкости.
^ + div(pU) = 0.	(3.1)
от
Выделим контрольный объем жидкости, ограниченный
параллелепипедом (см. Рис. 3.1). Часто задача может быть
сведена к двумерной постановке и тогда контрольным объ-
емом будет уже прямоугольник (см. Рис. 3.2).
Проинтегрируем уравнение (3.1) по рассматриваемому
объему
J^ZK + Jdiv(pU)c/K.	(3.2)
г ОТ у
Воспользовавшись теоремой Остроградского—Гаусса
(формула (1.13)), получим:
f^ + JpU-dS = O.	(3.3)
V Ul s
45
Глава 3. Основы метода конечных объемов
Рис. 3.1. Контрольный объем в трехмерной
постановке задачи.
Возможно, вы заметили, что когда мы выводили урав-
нения течения среды в п. 1.3, то мы тоже пользовались
теоремой Остроградского—Гаусса, однако — в обратном
направлении с тем, чтобы перейти от интегральных харак-
теристик к дифференциальным.
Уравнения вида (3.3) являются основными в методе ко-
нечных объемов. Фактически они представляют собой за-
коны сохранения для контрольного объема. Уравнение
46
3.1. Простой пример
(3.3) собственно и является уравнением закона сохранения
массы. Его физический смысл весьма прозрачен — масса
контрольного объема может измениться только за счет по-
тока жидкости, втекающей через его грани.
Если предположить, что профиль скорости на гранях
равномерный, а изменение плотности во всех точках объе-
ма происходит одинаково, то мы получим следующее дис-
кретное (алгебраическое) уравнение сохранения массы для
контрольного объема:
47
Глава 3. Основы метода конечных объемов
дт х	(3.4)
где Sxy, Syz, Sxz — площади соответствующих граней объе-
ма. Индексы in и out обозначают вход и выход соответст-
венно. Такая аппроксимация достаточно часто использует-
ся, однако она далеко не единственная.
Также здесь встает проблема, как выразить величины
атт™ атт™ атт™ ат т out ATT°ut at Tout
pu x >pUy >pU z 9 pu x >pU у , pUz через значения
pU в узловых точках (которые, как уже говорилось ранее,
находятся в середине каждого контрольного объема).
Более подробно все эти вопросы мы рассмотрим далее.
3.2.	Уравнение энергии
Рассмотрим уравнение (2.3) и выведем для него дис-
кретный аналог. При этом будем считать, что поле скоро-
стей в потоке нам заранее известно. Тогда задача сведется
лишь к отысканию поля температур. С задачей решения
всей системы уравнений (2.1) - (2.3) мы ознакомимся не-
сколько позднее.
Будем рассматривать задачу в двумерной постановке.
Введем сетку с шагами 8х и 8у, представленную на Рис. 3.3.
Контрольный объем здесь обозначен буквой Р и заштрихо-
ван. Стрелками показаны единичные вектора нормали к
его граням. Буквы N, NE, Е, ЕЕ и т.д. обозначают «стороны
света» относительно рассматриваемого объема. N — North
(север), Е — East (восток), S — South (Юг), W — West (за-
пад). Такое обозначение сторон контрольного объема для
ячеек прямоугольной формы является общепринятым и
часто используется в литературе (например, в [4,5]).
Интегрируя (2.3) по контрольному объему, получим:
J^JK + j(UV)TJK = JaATJK + J-^WK. (3.5)
V V	V	¥ РСр
49
Глава 3. Основы метода конечных объемов
Рис. 3.3. Контрольный объем в двумерной постановке
в декартовых координатах.
Применив теорему Остроградского—Гаусса (формула
(1.13)), а также вспомнив, что оператор Лапласа есть ди-
вергенция градиента (формула (1.9)), будем иметь:
+ ^7XJ dS = a^gradT dS + J-^Ч/Г.	(3.6)
V	S	S	V РСр
Перепишем его следующим образом:
J^JK = a^gradr-dS-^TUdS + J-^4/K.	(3.7)
V^T S	3	V РСр
Слева стоит член, отвечающий за изменение темпера-
туры в объеме по времени. Первый член в правой части
50
3.2. Уравнение энергии
представляет собой тепловой поток за счет теплопровод-
ности, второй — конвективный тепловой поток, а третий —
объемное тепловыделение.
Таким образом, перед нами стоит задача получения дис-
кретных аналогов объемных и поверхностных интегралов
уравнения (3.7). Кроме того нужно получить еще и дискрет-
ный аналог производной по времени и градиента.
Начнем с рассмотрения поверхностных интегралов.
3.3. Дискретные аналоги
поверхностных интегралов
Как известно интеграл по замкнутому контуру можно
представить в виде суммы интегралов по каждой из частей
контура:
( \
5	* V*
Jfds
(3.8)
В нашем случае мы имеем 4 стороны контрольного
объема. Таким образом, интеграл по замкнутому контуру
будет равен сумме интегралов по каждой из сторон кон-
трольного объема. Надо заметить, что в нашем случае f
может представлять собой либо TU, либо grad Г.
Далее мы будем рассматривать только одну из сторон
контрольного объема. Выражения для остальных сторон
записываются аналогично. Главное, на что здесь следует
обратить внимание — это на знак интеграла по той или
иной стороне. Этот знак определяется скалярным произве-
дением вектора на нормаль п к элементарной площадке dS
(см. формулы (1.1), (112) и комментарии к ним). Как из-
52
3.3. Дискретные аналоги поверхностных интегралов
вестно, косинус нуля равен 1. Это соответствует случаю,
когда векторы сонаправлены. Если векторы направлены
противоположно друг другу, то угол между ними состав-
ляет 180°, а косинус этого угла равен -1. Таким образом,
составляя дискретный аналог уравнений гидродинамики и
теплообмена, следует обращать внимание на то, как отно-
сительно рассматриваемой стороны контрольного объема
направлен вектор скорости. При этом направление вектора
скорости выбирается в самом начале. Обычно предполагают,
что проекции вектора скорости совпадают с направлениями
осей координат. Если в процессе вычисления скорость в ка-
кой-либо расчетной точке получается отрицательной, зна-
чит в этой точке направление скорости противоположно на-
правлению, которое было выбрано изначально.
Рассмотрим «восточную» сторону контрольного объе-
ма (см. Рис. 3.1).
Самая простая и очевидная аппроксимация интеграла
вида (3.8) основывается на предположении постоянства f
по всей поверхности рассматриваемой стороны:
\fdS = f. \dS = f,S'.	(3 9)
S.	S.
Эта аппроксимация имеет второй порядок точности [4].
Для повышения порядка точности можно использовать
формулу Симпсона (см., например, [7]):
53
Глава 3. Основы метода конечных объемов
(3.10)
5.	6
Здесь fne т/ifse — значения в угловых точках контроль-
ного объема. Формула (3.10) дает четвертый порядок точ-
ности [4].
Здесь встает закономерный вопрос об определении ве-
личин в формуле (3.9) (об использовании формулы (3.10)
мы поговорим ниже). Рассмотрим некоторые из вариантов.
При этом вспомним, что f может представлять собой кон-
вективный /с=фи или диффузионный fd = &аЛф пото-
ки (для уравнения энергии ф = Т).
Схема «против потока» (английское название —
upwind differencing scheme, UDS)
Она состоит в следующем:
<рр если (U • п)е > 0
<рЕ если (U-n) <0‘
(3.U)
(U п)е здесь показывает, совпадают или нет направле-
ния нормали к стороне контрольного объема и вектор ско-
рости ((U n)e >0 — вектора совпадают по направлению,
(U-n)e<0 — вектора направлены в противоположные
стороны). Схема называется схемой «против потока»
54
3.3. Дискретные аналоги поверхностных интегралов
именно поэтому: мы всегда выбираем значение, располо-
женное с «подветренной» стороны рассматриваемой части
контрольного объема (см. Рис. 3.4).
о 43	
	
и	
UH п => <ре = <Рр
U Н П => <Ре = Ve
Рис. 3.4 Схема «против потока».
Для дискретизации по схеме «против потока» диффу-
зионных членов будем иметь следующее выражение:
г д<р\
<8х)е
(РЕ~(РР
——— ,если
*Е ” *Р
——— ,если
(U-n)e>0
(Un)e<0.
(3-12)
55
Глава 3. Основы метода конечных объемов
Схема против потока хороша тем, что никогда не при-
водит к «нефизичным» решениям. Т.е., скажем, при тече-
нии холодного газа около горячей стенки использование
этой схемы никогда не приведет к такому решению, при
котором температура газа оказалась бы больше температу-
ры стенки. Она «гасит» колебания решения (вид таких ко-
лебаний показан на Рис. 3.5). Подробнее о свойствах этой
схемы см., например, [5].
Рис. 3.5. «Нефизичные» колебания решения.
Схема с линейной интерполяцией
Если предположить, что между двумя соседними узла-
ми (скажем, Р и Е) (р (в нашем случае — температура) из-
меняется по линейному закону, то можно будет записать
следующее выражение:
<Ре = <РЕ1А<РР^~1е\	(3.13)
56
3.3. Дискретные аналоги поверхностных интегралов
где 1е представляет собой интерполяционный коэффи-
циент:
1 Хе~ХР
1е=—----(3.14)
хЕ~хР
Отсюда можно получить формулу для определения
градиента, по сути совпадающую с формулой (3.12):
(3.15)
<^Х)е ХЕ ХР
Схема с линейной интерполяцией также имеет второй
порядок точности [4].
Схема с квадратичной интерполяцией
(QUICK — Quadratic Upwind Interpolation
for Convective Kinematics)
Данная схема является логическим продолжением схе-
мы с линейной интерполяцией. Она интерполирует пара-
метры между узлами с помощью парабол. Из курса чис-
ленных методов известно, что для построения полинома
степени п необходимо л+1 точек. Таким образом, для ин-
терполяции параболами требуется три точки. Исходя из
природы конвекции, третья точка выбирается с противопо-
точной стороны [4]. Тогда мы получим:
<Ре = <Ри + g\ (<Pd - <Ри) + ?2 (<Ри ~<PuuY (3.16)
57
Глава 3. Основы метода конечных объемов
Здесь D обозначает значение по потоку, U — значение
против потока. Например, для потока, текущего из W ъ Р,
D= W, U = Р, UU = Е. gi и g2 — интерполяционные коэф-
фициенты, зависящие от координат.
(3.17)
Дифференцируя по хе, получим:
' ду'
<дх)
= g\(<PD-<Pu)^g'2(<Pu-<Puu)^
(3.18)
gi _ хи 2хе + Хцу
(,XD ~ Xu)(XD ~ Х1
1 _ ХР ~ 2-Хе + Хц
6 2 ~ (
\ХР
(3.19)
и хии)
Схема с квадратичной интерполяцией несколько точ-
нее, чем схема с линейной интерполяцией. Она имеет
ошибку третьего порядка [4].
Схемы более высоких порядков
Очевидно, что возможны интерполяции и полиномами
более высоких порядков. Для использования формулы
58
3.3. Дискретные аналоги поверхностных интегралов
(3.10)	требуется интерполяция, по крайней мере, полино-
мом третьей степени:
<р(х) = а0+ ахх + а2х2 + а3х3.	(3.20)
Для подобной интерполяции, как и для интерполяций
полиномами более низких порядков, можно использовать
интерполяционный полином Лагранжа (см., например, [7]).
3.4.	Дискретные аналоги
объемных интегралов
Самым распространенным способом здесь является
предположение о том, что среднее по объему значение по-
дынтегральной функции равняется ее значению в узле Р.
Математически это может быть выражено следующим
образом:
\qdV = qV*qpV.	(3.21)
v
Формулы с интерполяцией применяются гораздо реже,
особенно в случаях, когда задача рассматривается в трех-
мерной постановке.
3.5.	Граничные условия
Для получения решения любого уравнения в частных
производных требуется задание условий на границе рас-
сматриваемой области.
При решении задач конвективного теплообмена выде-
ляют несколько типов граничных условий (boundary con-
ditions):
Граничное условие на входе (inlet). Обычно на входе
задается температура жидкости.
Граничные условия на стенке (wall). Здесь может быть
задана температура стенки (в русскоязычной литературе
такие граничные условия называют граничными условия-
ми первого рода), тепловой поток от стенки (граничные
условия второго рода) или коэффициент теплоотдачи (гра-
ничные условия третьего рода).
Граничные условия симметрии (symmetry). Такие усло-
вия применяются для задач, в которых существует сим-
метрия. Расчетная область ограничивается линией (плос-
костью) симметрии. Угол наклона касательной к профилю
температуры (а также всех других рассчитываемых пара-
61
Глава 3. Основы метода конечных объемов
метров) на линии (плоскости) симметрии равен нулю (см.
Рис. 3.6). Таким образом удается в два и более раз сокра-
тить объем рассматриваемой области, что снижает время,
требующееся на расчет.
Рис. 3.6. Течение в трубе — один из примеров
симметричного течения.
Периодические граничные условия (periodic boundary
conditions). Эти условия сходны с граничными условиями
симметрии. Их еще иногда называют условиями повторе-
ния. Они применяются для расчетов таких течений, где
существует много повторяющихся участков. Вместо того,
чтобы рассматривать всю картину течения, мы рассматри-
ваем только один участок, что в разы экономит вычисли-
тельные затраты. Примером может послужить расчет тече-
ния в кожухотрубном теплообменнике (линейная перио-
62
3.5. Граничные условия
дичность) или в лопаточной машине (радиальная перио-
дичность). Пример течения с периодичностью показан на
Рис. 3.7.
Рис. 3.7. Течение в кожухотрубном теплообменнике.
Вместо рассмотрения всего теплообменника
мы рассматриваем лишь обтекание одной трубки.
Граничных условий на выходе (outlet) по температуре
не задается.
Рассмотрим подробнее граничные условия на стенке.
Если задана температура стенки, то не возникает осо-
бенных трудностей, и нет необходимости в дополнитель-
ных уравнениях. Считается, что жидкость на поверхности
стенки имеет температуру стенки.
63
Глава 3. Основы метода конечных объемов
Если задан тепловой поток, то необходимо составить
для пристеночного контрольного объема, как и для всех
других, дискретный аналог уравнения энергии.
Будем считать контрольный объем около стенки поло-
винным (см. Рис. 3.8).
Рис. 3.8. Обычный и половинный контрольные объемы.
В данном случае при записи в дискретном виде уравне-
ния энергии для пристеночного объема изменятся только
выражения для поверхностных интегралов со стороны
стенки.
В уравнении (3.7) присутствует два поверхностных инте-
грала. Один из них
a J gradT dS
\ Стенка	у
есть поток тепла за
счет теплопроводности от стенки, а другой
j TUdS
УСтенка	J
поток тепла за счет конвекции. Поскольку на поверхности
стенки скорость среды равна нулю, то и | 7TJ - dS = 0.
Стенка
64
3.5. Граничные условия
Теперь, используя закон Фурье, выразим поток тепла за
счет теплопроводности. Вспомним общую формулировку
закона:
q = -2gradT.	(3.22)
Знак «минус» здесь означает, что поток тепла направ-
лен в сторону уменьшения температуры, т.е. тепло перете-
кает из более горячей области в более холодную.
Если рассматривать только поток тепла от стенки qw,
можно получить следующее равенство:
Под у здесь подразумевается ось координат перпенди-
кулярная стенке. Сопоставляя выражения для потока за
счет теплопроводности в уравнении (3.7) и формулу (3.23),
получим:
a J gradT dS =
Стенка
,	\	(3.24)
= а- [	= f dS = -^-5.
ОпеиюЛ / Р'СР Стенка РСР
Знак «минус» здесь исчез вследствие того, что направ-
ление теплового потока (мы положили его направленным
от стенки в поток) противоположно направлению нормали
к поверхности пристеночного контрольного объема. Если
направление теплового потока будет направлено от потока
65
Глава 3. Основы метода конечных объемов
в стенку (то есть стенка будет более холодной, нежели по-
ток), то значение qw нужно будет взять со знаком «минус».
Примечание: при выводе формулы (3.24) мы исходили
из постоянства qw по всей пристеночной грани контроль-
ного объема. Если тепловой поток от стенки не постоянен
по ее поверхности, а распределен по какому-либо закону,
то необходимо либо вычислять интеграл от qw по поверх-
ности, либо считать, что для каждого контрольного объема
^w=const. Во втором случае погрешность, вызванная таким
допущением, будет зависеть от закона изменения qw и от
выбранных размеров контрольных объемов.
3.6.	Производная по времени
Время в физике может рассматриваться как своеобраз-
ная «четвертая координата», с той лишь разницей, что бу-
дущее и прошлое четко разграничены между собой. Про-
изошедшее в настоящий момент может повлиять только на
будущие события, но не на прошлые (так называемый
принцип причинности).
Исходя из приведенной аналогии, при численном ре-
шении задач время процесса можно разбить на некоторое
количество шагов, создав наряду с пространственной сет-
кой временную.
Существует два основных метода решения нестацио-
нарных задач. Метод, в котором неизвестные величины на
текущем временном шаге выражаются через величины
предыдущего временного шага, которые уже известны, на-
зывается явным. В неявном же методе величины на теку-
щем временном шаге выражаются друг через друга.
Явный метод более прост, однако на него накладыва-
ются дополнительные ограничения по временному шагу, в
то время как никаких ограничений (кроме физических) на
временной шаг в неявном методе не накладываются.
67
Глава 3. Основы метода конечных объемов
Явный метод Эйлера (Euler explicit method)
Название метода отражает тот факт, что при его ис-
пользовании получаются явные зависимости искомой ве-
личины от величин уже известных. То есть, нам не требу-
ется производить никаких дополнительных действий для
получения решения.
Рассмотрим явный метод на примере одномерного
уравнения энергии в предположении, что внутренние ис-
точники тепла отсутствуют (qv = 0):
дТ тт дТ д2Т
дт х дх дх2
(3.25)
Введем равномерную сетку с шагом дх и зададимся
шагом по времени Дт. Используя схему с линейной интер-
поляцией, получим (верхним индексом будем отмечать
шаги по времени, нижним — номера ячеек сетки):
(3.26)
После алгебраических преобразований получим:
+f Fo- —1	+ f Fo+—)
L 2 J 1 v 2 J '
(3.27)
»л-1
68
3.6. Производная по времени
Здесь Fo = -—-у — сеточное число Фурье (cell Fourier
\Sx)
number), Со = —— — сеточное число Куранта (cell
ох
Courant number). Эти два безразмерных параметра являют-
ся одними из наиболее важных при решении нестационар-
ных задач гидродинамики и теплообмена.
Физически число Фурье представляет собой отношение
временного шага Дт к времени, за которое тепловая волна
распространится на расстояние дх.
Число же Куранта есть отношение временного шага Дт
к характерному времени конвекции 8 x/Ux.
Таким образом, эти числа отражают, насколько велик
шаг времени по отношению к характерному времени про-
текания процессов в системе. Если мы выберем слишком
большой временной шаг, то, очевидно мы просто «пропус-
тим» (т.е. не будем учитывать) какие-то процессы в систе-
ме. Это может привести к нефизичности решения.
Пределы значений для Fo и Со могут быть получены в
общем виде, однако мы здесь приводить эти выкладки не
будем, а ограничимся следующими простыми рассужде-
ниями. Как видно из (3.27), температура на временном ша-
ге п в точке i зависит от температур на предыдущем вре-
менном шаге в точках i- 1, i, i + 1. Если коэффициенты
69
Глава 3. Основы метода конечных объемов
перед этими температурами будут иметь произвольный
знак, то это будет означать, что рост температуры в точках
i и i + 1 на предыдущем временном шаге может привести к
падению температуры в точке i на текущем временном ша-
ге, что противоречит принципам термодинамики.
Таким образом, мы пришли к искомым ограничениям
по значению чисел Fo и Со для явной схемы при использо-
вании схемы с линейной интерполяции:
Fo< |,
Co<2Fo.
Решая систему этих неравенств, получим:
Дг<-^—
2а
(3.28)
(3.29)
То есть мы получили, что для случая конвективного те-
плообмена должен быть ограничен не только шаг по вре-
мени, но и шаг пространственной сетки. В случае, если
конвекция отсутствует (Цх= 0), ограничение на шаг сетки
снимается (из второго неравенства в (3.29) получим, что
Зх < оо э что выполняется всегда). Разумеется, это не озна-
чает, что ограничения на шаг сетки не существует вовсе,
однако они определяются другими параметрами.
70
3.6. Производная по времени
Если преобразовать второе неравенство в (3.29), можно
получить безразмерный параметр, называемый числом се-
_ Ux8x
точным числом Пекле Ре = —~— . Таким образом, для
задач конвективного теплообмена при использовании яв-
ного метода и схемы с линейной интерполяцией должно
выполняться условие Ре < 2.
При использовании других схем дискретизации (метода
«против потока», методы с квадратичной интерполяцией
и т.д.) условия, налагаемые на шаги по времени и по про-
странственной сетке должны так же исходить из положи-
тельности коэффициентов при температурах в различных
точках.
Неявный метод Эйлера (Euler implicit method)
В неявной схеме искомая величина выражается через
величины на том же временном шаге. Для схемы с линей-
ной интерполяцией получим:
Т”=Т”-'
|д Т^-2Т”+Т^
26 х
Дт. (3.30)
2
После преобразований, получим:
(i+2Fo)7--+( v-F»]7;:,+(-v-F°F
k х J k х J
= Г-1. (3.31)
71
Глава 3. Основы метода конечных объемов
Как уже было упомянуто выше, данный метод не нала-
гает ограничений на временной шаг (ограничения на вре-
менной шаг налагаются чисто физические — если времен-
ной шаг будет больше характерного времени протекания
какого-либо процесса, то мы его попросту «пропустим»).
Однако он требует больших вычислительных затрат. Эти
затраты особенно ощутимы в задачах с большим количест-
вом ячеек в сетках. К примеру, стандартные задачи на
внешнее обтекание могут потребовать расчетных сеток с
несколькими миллионами ячеек. Это означает, что при ис-
пользовании неявного метода потребуется решить систему
из нескольких миля ионов уравнений! Даже при нынешнем
уровне развития вычислительной техники это потребует
значительного времени.
Существуют и другие методы получения дискретных
аналогов производной по времени. Мы не будем приво-
дить их здесь. Интересующиеся могут найти их в [4,5, 8].
Метод установления
Этот метод часто применяется для решения стационар-
ных задач. Он заключается в том, что стационарное реше-
ние получается из решения нестационарной задачи. Счита-
ется, что система, у которой существует стационарное со-
стояние, стремится прийти в него. Таким образом, если мы
будем решать нестационарную задачу на достаточно
72
3.6. Производная по времени
большом промежутке времени, то, в конце концов, полу-
чим стационарное решение. Начальные же условия задачи
при этом могут быть практически любыми.
Метод установления хорош в первую очередь тем, что
позволяет проверить, есть ли у системы стационарное со-
стояние как таковое. Также его часто используют в пакетах
прикладных программ, поскольку иначе пришлось бы со-
ставлять отдельные модули для решения стационарных и
нестационарных задач, а метод установления позволяет
использовать один модуль для решения задач обоих типов.
Обычно для решения задач методом установления при-
меняют неявные схемы. Они позволяют задать достаточно
большие шаги по времени, что приводит к более быстрому
получению искомого стационарного решению.
3.7.	Общий алгоритм
решения задачи
Мы рассмотрели здесь основные вопросы, связанные с
решением задач вычислительной механики сплошной сре-
ды с помощью метода конечных объемов. Теперь сведем
все в единую картину.
В общем и целом составление системы приближенных
алгебраических уравнений из исходных уравнений в част-
ных производных представляет сборку из отдельных гото-
вых «деталей». Надо лишь выбрать, какие «детали» лучше
подойдут для данной конкретной задачи, и учесть требова-
ния, которые при этом предъявляются.
Приведем общий алгоритм, которому надо следовать
при составлении численной схемы:
Постановка задачи. Начинать следует с выбора тех
уравнений, которые мы будем решать. Что нам нужно най-
ти и что нам известно при этом? Необходимо также опре-
делиться с точной геометрией расчетной области, типом
задачи (стационарная или нестационарная), граничными и
начальными условиями.
74
3.7. Общий алгоритм решения задачи
Выбор численной схемы. Сюда входит выбор схемы
аппроксимации («против потока», схемы с интерполяци-
ей и т.п.), выбор метода получения дискретного аналога
производной по времени (если задача нестационарная
или выбран метод установления решения стационарной
задачи).
Построение сетки и выбор шагов по времени. После
выбора схемы следует учесть все ограничения, которые
накладываются ею на размер шага сетки и временного ша-
га. Следует определиться, будет ли сетка равномерной или
нет. Если сетка равномерна, то следует выбрать шаг. Если
неравномерна, то следует собрать в единый массив коор-
динаты всех расчетных узлов.
Построение численной схемы. Ячейки сетки разде-
ляются на группы (ячейки внутри расчетной области и
ячейки на её границах). Для каждой группы ячеек запи-
сываются проинтегрированные по контрольному объему
(т.е. вида (3.7)) уравнения. В соответствии с выбранной
схемой аппроксимации следует расписать все члены, вхо-
дящие в уравнение. При рассмотрении поверхностных
интегралов особое внимание следует обращать на знак
(подробнее см. п. 3.3).
Решение полученной системы уравнений. В результа-
те построения численной схемы обычно получают систему
линейных алгебраических уравнений. Решение может
75
Глава 3. Основы метода конечных объемов
осуществляться с помощью различных методов. Основные
методы решения систем линейных уравнении рассматри-
ваются в курсе линейной алгебры (см. например, [7]).
Проверка адекватности полученного решения. Ре-
зультаты расчетов должны быть проверены с точки зрения
соблюдения законов сохранения для области в целом, а
также с точки зрения их общей физичности. Например, ес-
ли в ходе решения задачи о течении воды в канале с тепло-
обменом расход воды на выходе больше расхода на входе,
или температура воды превышает температуру насыщения
(кипения), то об адекватности полученных результатов
расчета необходимо задуматься.
Глава 4
Расчет поля течения

В разделе 3 мы рассмотрели общие принципы метода
конечных объемов. При этом мы считали, что поле скоро-
стей нам известно. Однако на практике такое встречается
редко. Чаще всего поле скоростей заранее неизвестно и его
требуется найти.
В данном разделе мы рассмотрим вопросы, касающиеся
решения полной системы уравнений Навье—Стокса. Пока
будем считать, что течение у нас ламинарное. Вопросы,
связанные с проявлением эффекта турбулентности будут
рассмотрены позднее.
4.1. Интегрирование
уравнения движения
Мы уже рассматривали интегрирование по контроль-
ному объему уравнений неразрывности и энергии. Теперь
настала очередь уравнения движения:
V	V
= J JdV --fvPdV + vj АШГ.
У Р V	У
(4.1)
Воспользовавшись формулой (1.14), преобразуем инте-
грал от градиента давления. Будем иметь:
V	У
= J JdV - - JPdS + V J АШГ.
У Р S	У
(4-2)
Данное уравнение — векторное. Его можно свести к
трем обычным (скалярным) уравнениям, спроецировав на
оси координат. Запишем лишь одну такую проекцию:
80
4.1. Интегрирование уравнения движения
J	jv + j div(C/xU) dV =
у дт	у
(4.3)
v Р s v
Отсюда сразу видно, что мы можем преобразовать
некоторые интегралы по объему к поверхностным инте-
гралам. Применив теорему Остроградского—Гаусса,
получим:
+ dS =
V	S
(4.4)
V
Ps
S
Чтобы получить второе и третье уравнение, надо х за-
менить соответственно на у и z.
В тензорной нотации то же самое запишется в виде:
^\u,dv+\ufj,-ds^
V	S
=	-dSt.
i Pi	idx,
(4.5)
Если бы мы сразу воспользовались тензорной нотаци-
ей, то нам бы не пришлось записывать проекцию уравне-
ния, а мы бы сразу получили окончательный результат.
Обычно в литературе тензорная нотация применяется
совместно с обычной.
4.2.	Расчет поля давления
Когда мы рассматриваем течение сжимаемой жидкости,
то уравнение неразрывности имеет ввд (3.1). В это уравне-
ние входит плотность жидкости р. Зная уравнение состоя-
ния, можно легко найти поле давления из поля плотности.
Чаще всего уравнением состояния является уравнение
Менделеева—Клапейрона:
Р = pRT.	(4.6)
В системе же уравнений течения несжимаемой жидко-
сти (2.1) - (2.3) давление присутствует только в виде сво-
его градиента в уравнении (2.2), из которого мы будем на-
ходить составляющие скорости. При этом уравнение для
явного определения поля давления отсутствует и способ
нахождения давления далеко не очевиден.
К данной проблеме существует несколько подходов.
•	Исключение давления из системы уравнений.
•	Получение уравнения типа уравнения Пуассона для
давления.
82
4.2. Расчет поля давления
Методы исключения давления
из системы уравнений
Вспомним известную формулу векторного анализа
(1.11) и запишем ее для вектора скорости:
(UV)U = grad
+ rotUxU.
(4.7)
Подставив (4.7) в (2.2) и применив операцию rot, полу-
чим уравнение для вихря '.
д	.Г—JZ - rot (U х rot U) = v • Д (rot U) + rot J.	(4.8)
Это уравнение вместе с уравнением (2.1) составляют
систему для нахождения вектора скорости. Их можно ре-
шить, используя метод конечных объемов и применяя из-
вестные соотношения векторного анализа (см. п. 1).
К недостаткам этого подхода можно отнести некото-
рую сложность в понимании физического смысла полу-
ченных уравнений (что, впрочем, субъективно).
Кроме того, можно модифицировать этот подход вве-
дением функции тока. Этот прием широко используется
для описания двумерных течений.
Функцию тока у/ определяют следующим образом:
Ux^; Ur=-?V
ду	дх
(4.9)
83
Глава 4. Расчет поля течения
Очевидно, что при таком рассмотрении уравнение не-
разрывности удовлетворяется автоматически (проверьте
это, подставив (4.9) в (2.1)).
Для ротора скорости в двумерном случае мы имеем
следующее выражение:
диу дих
rotU = a> = —(4.10)
дх ду
Принимая во внимание формулу (4.9), можно пока-
зать, что:
й? = -Дул	(4.11)
Окончательно получим следующую форму уравнения
движения:
да) (да) ду/ да)	_
— + -—=v-A®+rotJ. (4.12)
дт ^Зх ду ду дх)	v ’
Преимуществом введения переменных «функция тока —
вихрь» является то, что вместо трех уравнений (уравнение
неразрывности и две проекции векторного уравнения дви-
жения), мы решаем два — уравнение (4.11) и (4.12). Если
подставить (4.11) в (4.12), то можно вообще свести задачу
к одному уравнению четвертого порядка. После получения
решения для введенных переменных, с их помощью вы-
числяются так называемые естественные переменные (ско-
рость, давление и т.п.).
84
4.2. Расчет поля давления
Серьезным недостатком этого подхода является неоче-
видность граничных условий. Он более интересен со своей
математической стороны, однако несет мало физического
смысла. Поэтому метод конечных объемов в данном слу-
чае применить тоже несколько затруднительно. Опять же,
восстановление значений естественных переменных по по-
лю значений функции тока и вихря — тоже не совсем про-
стая задача.
Желающих подробнее рассмотреть методы решения за-
дач с помощью введения переменных «функция тока —
вихрь» мы отсылаем к [9] и больше останавливаться на
них не будем.
Уравнение Пуассона для давления
Используя уравнения неразрывности и движения, мож-
но получить уравнение, из которого можно явным образом
найти поле давления.
Для этого применим к обеим частям уравнения (2.2)
операцию дивергенции:
div(gradP) = р• divf J - — - (U • V)U + vAU |. (4.13)
\ от	)
Произведем некоторые преобразования, учитывая фор-
мулы (1.11), а также тот факт, что:
85
Глава 4. Расчет поля течения
div
(4.14)
Получим следующее выражение:
ДР = p(div J —— (divU) -
дт	(4.15)
-div[(U • V)u] + vdiv(AU)).
Примем во внимание тот факт, что по формуле (1.10)
AU = grad(divU)-rot(rotU). Но дивергенция скорости со-
гласно уравнению неразрывности (2.1) равна нулю. Ди-
вергенция любого ротора так же равна нулю (см. формулы
(1.11)). Таким образом, мы приходим к следующему урав-
нению:
\Р = р {div J-div [(U-V) и]}.	(4.16)
Это и есть искомое уравнение для давления. Из него в
частности видно, что для несжимаемой жидкости поле
давления от вязкости не зависит.
Запишем его так же и в тензорной нотации:
д (дРУ
дх} \дхи
dJt
dxt
д (ди(и/
dxt >
(4-17)
Таким образом, система уравнений, описывающих не-
изотермическое течение несжимаемой жидкости состоит
86
4.2. Расчет поля давления
из уравнений (4.16), (2.2) и (2.3), причем уравнение (4.16)
заменило здесь уравнение неразрывности (2.1).
Проинтегрируем уравнение (4.16) по контрольному
объему:
^gradP dS = p |j-dS-|(U-V)lbdS
S’	Ls s
(4.18)
Тут следует сделать важное замечание: для получения
корректных результатов способ, которым получаются дис-
кретные аналоги интегралов в уравнении (4.18) и (4.5)
должен быть одним и тем же. Проще говоря, если для ин-
теграла по градиенту давления в (4.18) выбран способ с
линейной интерполяцией, то и в уравнении (4.5) интеграл
от давления должен быть аппроксимирован с помощью
линейной интерполяции (но никак не с помощью квадра-
тичной или кубической интерполяции или какого-либо
другого метода).
Все это обуславливается тем, что мы получили уравне-
ние для давления путем преобразования уравнения движе-
ния и применения уравнения неразрывности. Таким обра-
зом, напрямую уравнения неразрывности мы не решаем,
оно выполняется автоматически. Однако если при числен-
ном решении уравнений мы будем применять для каждого
из них неодинаковые способы, то уравнение неразрывно-
сти выполняться не будет.
4.3.	Подходы и проблемы
Как мы уже показали выше, решение каждого из урав-
нений системы (2.1) - (2.3) или (4.16), (2.2) и (2.3) сводит-
ся к решению системы линейных алгебраических уравне-
ний. И уравнений этих столько, сколько ячеек в расчет-
ной сетке.
Каждую из систем этих уравнений можно решать как
совместно со всеми другими, так и по отдельности. Когда
мы решаем уравнения по отдельности (т.н. раздельный
подход, segregated approach), мы рассматриваем только ка-
кое-то одно уравнение системы (2.1) - (2.3). Так мы посту-
пали в п. 3, считая, что поле скорости нам известно. Этим
мы лианеризовали уравнение энергии. Проще говоря, свели
нелинейное уравнение к линейному.
Раздельное решение требует меньше памяти, хотя при
этом может возникнуть ряд проблем.
Во-первых, при раздельном решении требуется приме-
нять особые приемы для того, чтобы решения, получаемые
для поля давления и для поля скорости, соответствовали
друг другу (в англоязычной терминологии это называется
88
4.3. Подходы и проблемы
термином pressure-velocity coupling, т.е. — «стыковка ско-
рости и давления»). Для этой «стыковки» часто вводят для
каждого уравнения системы свою собственную особую
сетку (хотя есть подходы, где для всех уравнений исполь-
зуется одна сетка) и специальные формулы для корректи-
ровки давления и скорости. Два популярных метода со-
кращенно называются SIMPLE (Semi-Implicit Method for
Pressure-Linked Equations) и SIMPLER (SIMPLE Revised).
Подробно эти методы описаны в [5].
При совместном решении (coupled approach) все урав-
нения решаются совместно. Это усложняет задачу, по-
скольку все уравнения в системе (2.1) - (2.3), кроме урав-
нения (2.1) в этом случае нелинейны. Для их решения тре-
буется больше ресурсов, однако при использовании этого
подхода часто удается добиться более быстрой сходимости
решения, чем при раздельном решении.

Глава 5
Турбулентность:
проблемы
моделирования
и подходы
к их решению

До этого мы не касались вопроса режима течения. Мы
лишь упомянули, что уравнения, представленные в п. 4 в
такой форме могут быть успешно применены в основном
лишь для ламинарных потоков. Данная глава посвящена
проблеме турбулентности и основным подходам к ее мо-
делированию.
Все дифференциальные уравнения здесь записываются
только в тензорной нотации, поэтому читателю рекомен-
дуется ещё раз просмотреть п. 1.2.
5.1.	Феномен турбулентности
Начать рассмотрение феномена турбулентности проще
всего именно с того самого исторического примера, кото-
рый и положил начало её изучения. Речь идет об опыте
О. Рейнольдса 1883 года.
Экспериментальная установка показана на Рис. 5.1 и
Рис. 5.2. Она состоит из резервуара А с водой, от которого
отходит стеклянная труба В с краном С на конце, и сосуда
D с водным раствором краски, которая может по трубке
вводиться тонкой струйкой внутрь стеклянной трубы В.
Рис. 5.1. Опыт Рейнольдса; ламинарный режим течения.
94
5.1. Феномен турбулентности
Рис. 5.2. Опыт Рейнольдса; турбулентный режим течения.
Первый случай движения жидкости. Если немного
приоткрыть кран С и дать возможность воде протекать в
трубе с небольшой скоростью, а затем с помощью крана Е
впустить краску в поток воды, то увидим, что введенная в
трубу краска не будет перемешиваться с потоком воды.
Струйка краски будет отчетливо видимой вдоль всей стек-
лянной трубы, что указывает на слоистый характер тече-
ния жидкости и на отсутствие перемешивания. Такой ре-
жим движения называется ламинарным (Рис. 5.1).
Второй случай движения жидкости. При постепенном
увеличении скорости течения воды в трубе путем откры-
тия крана С картина течения вначале не меняется, но затем
при определенной скорости течения наступает быстрое ее
изменение. Струйка краски по выходе из трубки начинает
колебаться, затем размывается и перемешивается с пото-
ком воды, причем становятся заметными вихреобразова-
95
Глава 5. Турбулентность: проблемы моделирования
ния и вращательное движение жидкости. Такое течение
называется турбулентным (Рис. 5.2).
Если бы мы замеряли скорость в каком-либо сечении
трубки в первом и втором случае, то увидели бы, что в
первом случае она оставалась бы величиной постоянной, а
во втором — совершала бы нерегулярные колебания около
некоторого среднего значения.
Рейнольдс обнаружил, что переход течения от лами-
нарного режима к турбулентному наблюдается при неко-
тором значении скорости, которую принято называть кри-
тической. Значение критической скорости зависит от соот-
ношения диаметра трубки, плотности среды и вязкости.
Обобщая данные опытов, Рейнольдс ввел безразмерный
критерий, названный впоследствии в его честь:
П pud
Re = ——.	(5.1)
Р
Физический смысл критерия Рейнольдса заключается
в том, что он представляет собой соотношение сил инер-
ции и сил вязкости. Если силы инерции превышают на
какой-то порядок (определяемый критическим числом
Рейнольдса) силы вязкости, течение переходит к турбу-
лентному режиму.
Опыт показал, что критическое число Рейнольдса для
течения в круглой трубе Re,,, примерно равно 2300.
96
5.1. Феномен турбулентности
Таким образом, критерий подобия Рейнольдса позволя-
ет судить о режиме течения жидкости в трубе. При
Re < Re«p течение является ламинарным, а при Re > Re«p
течение является турбулентным. Точнее говоря, вполне
развитое турбулентное течение в трубах устанавливается
лишь при Re примерно равному 4000, а при Re =
= 2300.. .4000 имеет место переходная область.
Так в чем же феномен турбулентности? Почему после
того, как силы инерции во сколько-то раз превышают силы
вязкости, в потоке начинают образовываться вихри, а все
параметры колеблются?
Все дело в том, что когда силы вязкости малы по отно-
шению к силам инерции, течение теряет устойчивость.
Ведь на самом деле случайные и хаотические процессы
происходят в жидкости постоянно, даже в условиях, когда
она никуда не движется или течет в ламинарном режиме
(вспомните, что такое броуновское движение). Однако, ко-
гда силы вязкости достаточно высоки, любые флюктуации
очень быстро угасают. Вязкость не дает им развиться. На-
против, когда вязкость мала, любое случайное возмущение
не угасает, а наоборот может даже усиливаться. Так про-
исходит вихреобразование, а течение становится турбу-
лентным.
5.2.	Явления отрыва
Турбулентность может возникнуть также в случае
внешнего обтекания тел. Рассмотрим, например, обтекание
цилиндра.
При Re < 5 поток обтекает цилиндр в ламинарном ре-
жиме, полностью огибая его (см. Рис. 5.3).
Рве. 53. Обтекание цилиндра.
Безотрывный режим. Re = 1,54.
Начиная с Re ~ 5, за цилиндром происходит отрыв тече-
ния и формируются вихри. По мере увеличения числа Re
98
5.2. Явления отрыва
они вытягиваются (Рис. 5.4 и Рис. 5.5). При числе Re порядка
100 вихри начинают отрываться и двигаться вниз по потоку,
образуется так называемая дорожка Кармана (Рис. 5.6).
Рис. 5.4. Отрывное обтекание цилиндра, Re = 13,1.
Рис. 5.5. Отрывное обтекание цилиндра, Re - 26.
99
Глава 5. Турбулентность: проблемы моделирования
Рис. 5.6. Дорожка Кармана при обтекании цилиндра. Re <= 140.
Цилиндр — гладкое тело, однако уже при малых числах
Рейнольдса обтекающий их поток отрывается, и образуются
вихри. А что будет, если мы рассмотрим какое-нибудь бо-
лее плохообтекаемое тело, например квадратный выступ на
стенке? Очевидно, что отрыв на нем будет происходить да-
же при очень малых числах Re (Рис. 5.7). Это используется
в технике, когда нужно специально турбулизовать поток.
Рис. 5.7. Обтекание прямоугольного выступа на пластине.
5.3.	Осредненное движение.
Уравнения Рейнольдса
Как же можно описать турбулентное движение? Вооб-
ще, строго говоря, никаких новых уравнений не нужно.
Турбулентное течение вполне может быть описано обыч-
ной системой уравнений Навье—Стокса, поскольку при их
выводе не налагалось никаких ограничений на соотноше-
ние между силами инерции и силами вязкости. Однако,
если мы каким-то обычным способом попытаемся числен-
но решить эти уравнения для турбулентного режима тече-
ния (этот подход называется прямым численным модели-
рованием, Direct Numerical Simulation, DNS), то у нас, ско-
рее всего ничего не получится. Для прямого решения
уравнений Навье—Стокса в случае турбулентного режима
течения нужны очень качественные сетки с очень малым
шагом. Шаг по времени также должен быть очень мал. Все
это потребует очень мощных вычислительных ресурсов с
большим объемом памяти и быстродействием. И потреб-
ность в вычислительных ресурсах растет примерно про-
101
Глава 5. Турбулентность: проблемы моделирования
порционально кубу из числа Рейнольдса в рассматривае-
мом течении (подробнее о DNS см. [10, 11]). Таким обра-
зом, в инженерных и научных расчетах, часто возникает
необходимость как-то упростить систему уравнений На-
вье—Стокса, снизив тем самым потребности в вычисли-
тельных мощностях.
Для количественного описания развитого турбулентно-
го течения Рейнольдс предложил следующий, получивший
широкое применение прием. Регистрируя во времени па-
раметры потока можно предположить, что значение каж-
дого из них можно представить в виде суммы осредненной
(обозначается чертой сверху) и пульсационной (обознача-
ется штрихом) составляющих (см. Рис. 5.8). Например,
возьмем скорость:



Рис. 5.8. Пульсации осевой скорости.
102
5.3. Осредненное движение. Уравнения Рейнольдса
UX=UX+U'X; Uy=Uy+U'y; UZ=UZ+U'Z. (5.2)
То есть мы будем рассматриваем ее как некую не ме-
няющуюся (или меняющуюся по какому-то определенному
закону) величину к которой добавлены случайные пульса-
ции. Таким образом, мы рассматриваем турбулентное тече-
ние, как некий случайный процесс, применяя к нему приемы,
используемые в теории вероятностей и математической ста-
тистике.
Указанное выше осреднение имеет следующие свойства:
(р = Ф,	(5.3)
(p + W = (p + W,	(5.4)
const = const (p,	(5.5)
<p-y = <p	(5.6)
(p' = <p-<p = Q,	(5.7)
(5.8)
^дх) дх
Представив в виде (5.2) скорость и температуру, подста-
вив их в уравнения системы (2.1) - (2.3) и выполнив неко-
торые преобразования (здесь и далее подразумевается, что
жидкость несжимаема, т.е. р = const, получим уравнения
Рейнольдса (RANS — Reynolds Averaged Navier-Stokes):
103
Глава 5. Турбулентность: проблемы моделирования
РСР
dUt
дх.
О,
~dUt ! д(Ц-Ё7,)~| др
dxt
(5.9)
дт
д
дх.
dUt
dxj
(5.Ю)
(5.И)
Данная система уравнений содержит девять неизвест-
ных членов (шесть членов вида -p-UjUj и три члена вида
-рсрйуг).
Если посмотреть на структуру уравнений, то становит-
ся понятен физический смысл этих величин. Член
- р UPj представляет собой тензор турбулентных напря-
жений трения или тензор рейнольдсовских напряжений
(Reynolds stress tensor). В самом деле, если расписать его,
мы получим девять величин:
	[ й* и X	U'U' Х	у	U'xU'z	
	U' U'r	и* У	U',U'Z	(5.12)
	U'ZU'X	U'zU'y	J	
104
5.3. Осредненное движение. Уравнения Рейнольдса
Следует заметить, что на самом деле эти дополнитель-
ные напряжения (равно как и дополнительный тепловой
поток) появились только потому, что мы представили ско-
рость (и температуру) в виде (5.2). Чисто физически, кроме
сил трения, обусловленных молекулярной вязкостью, дру-
гих сил трения в потоке нет.
Напряжений в этом тензоре (5.12) девять, однако, оче-
видно, что UjUj = UjUi, поэтому неизвестных величин
только шесть.
Член -pcpU'jT' представляет собой дополнительный
поток тепла за счет турбулентного переноса.
Чтобы вычислить эти члены, требуются дополнитель-
ные уравнения. Эти уравнения обычно называют моделью
турбулентности.
Выбор модели турбулентности для конкретной задачи —
вопрос непростой. Необходимо учитывать свойства каж-
дой модели и имеющиеся в распоряжении исследователя
вычислительные мощности.
В последующих параграфах мы изложим самые про-
стые концепции и подходы к моделированию турбулент-
ности и приведем иерархию существующих на сегодняш-
ний день моделей.
5.4.	Гипотеза Буссинеска
Сама по себе гипотеза Буссинеска (Boussinesq) не яв-
ляется моделью турбулентности. Она лишь утверждает,
что турбулентные напряжения, как и обычные напряжения
трения пропорциональны градиенту скорости. Если в слу-
чае обычных напряжений трения коэффициентом пропор-
циональности был коэффициент динамической вязкости р
(dynamic viscosity), то в случае напряжений Рейнольдса
это коэффициент турбулентной вязкости рч (eddy viscosi-
ty). То есть влияние турбулентности на течение учитыва-
ется с помощью дополнительной вязкости, возникающей в
потоке.
В случае течений, когда две проекции скорости малы
по сравнению с третьей (например, течение в трубе или
течение в пограничном слое) гипотеза Буссинеска записы-
вается очень просто (ст в данном случае — полное напря-
жение трения; a — т.н. эффективная вязкость):

дЦх
ду
ay
(5.13)
106
5.4. Гипотеза Буссинеска
Ддя общего случая течения несжимаемой жидкости ги-
потеза Буссинеска записывается следующим образом:
------	(ди. ди
-pU'U^ft- -T-L+-T^
дх, дх.

(5-14)
В данном уравнении $у — символ Кронекера (см.
U'U'-
(1.22)),	а к-—1—*- — кинетическая энергия турбу-
лентности.
Подставив (5.14) в (5.10), после преобразований по-
лучим:
р
ди, t =
дт dxj
дР д ,
т- + 7— '
dxt дх}
dUt
дх.
(5.15)
(- 2	А
Здесь Р = Р + ~zPk I — давление с учетом турбулент-
\	«э	J
ных пульсаций.
Гипотезу Буссинеска можно применить и к уравнению
энергии. В этом случае она будет выглядеть следующим
образом:

(5.16)
107
Глава 5. Турбулентность: проблемы моделирования
То есть дополнительный турбулентный тепловой поток
выражается через дополнительную турбулентную теплопро-
водность среды. В этом случае уравнение (5.11) примет вид:
р4^+^;^1 = /“к1+;1т)^1+‘г'' (5Л7)
дт J dxj dxj	дх}
По аналогии с обычным числом Прандтля (которое ус-
танавливает связь между молекулярными коэффициентами
вязкости и теплопроводности), можно ввести турбулентное
число Прандтля, устанавливающее связь между соответст-
вующими турбулентными коэффициентами:
Ргг=-^.
А-
Для многих течений жидкостей и газов (например — те-
чения в трубах, течения в соплах, струи) предполагают, что
Ргг = 0.9. При заданном турбулентном числе Прандтля
необходимо найти только один из коэффициентов, а вто-
рой тогда получается автоматически. Обычно ищут турбу-
лентную вязкость, а турбулентную теплопроводность вы-
числяют через Ргг.
Подходы, использующие гипотезу Буссинеска, имеют
один существенный недостаток. Коэффициент турбулент-
ной вязкости есть скаляр. Таким образом, предполагается,
что турбулентность изотропна, т.е. не зависит от направ-
ления, что верно далеко не всегда.
5.5.	Путь смешения Л. Прандтля.
Алгебраические модели
Модель для описания распределения турбулентной ки-
нематической вязкости
у -El
vt -
Р
впервые была предло-
жена Л. Прандтлем в 1925г. и известна как модель пути
смешения.
Рассматривая осредненные сдвиговые течения без гра-
диента давления, Прандтль постулировал, что характерный
масштаб пульсаций скорости равен градиенту осредненной
скорости, умноженному на характерный масштаб длины
1т, который он назвал путем смешения.
Понятие пути смешения в теории турбулентности
сродни понятию свободного пробега молекулы в кинети-
ческой теории газов. Путь смешения есть некоторая длина,
на которой вихрь теряет свою индивидуальность, смеши-
ваясь с окружающим потоком. Это расстояние также равно
среднему расстоянию пульсаций.
В свете гипотезы пути смешения можно получить сле-
дующее выражение для турбулентной вязкости:
109
Глава 5. Турбулентность: проблемы моделирования
dU

ду
(5.19)
Длина пути смешения определяется эмпирически. Ус-
пех предложенной Прандтлем модели был предопределен
тем обстоятельством, что для многих простых типов тече-
ний со сдвигом 1т может быть выражена относительно не-
сложными формулами. При рассмотрении течения в по-
граничном слое полагают
1т=ку.	(5.20)
Здесь к = 0.4 —постоянная Кармана.
С помощью модели Прандтля было получено много
важнейших для теплообмена и гидродинамики результа-
тов. Однако сейчас она мало применяется, поскольку при-
годна лишь для простых типов течений (скажем, отрывные
течения она описывает плохо).
Модель пути смешения Прандтля является алгебраиче-
ской моделью. Алгебраические модели отличаются тем,
что для получения значения турбулентной вязкости не
требуется решать дополнительных дифференциальных
уравнений. Моделей этой группы для разных случаев соз-
дано достаточно много. Хороший их обзор дан в [10].
5.6.	Модель к-е.
Дифференциальные модели
Модели с дифференциальными уравнениями являются
более сложными, однако они и более универсальны, их
можно применять к широкому кругу задач.
Самой известной и популярной на сегодняшний день
дифференциальной моделью является модель к-с. Для ее
построения вводятся два параметра: К	"— — кине-
dU\ dU\
тическая энергия турбулентности и е ~v~q^	—
диссипация турбулентной энергии. Фактически уравнения
модели турбулентности представляют собой закон сохра-
нения турбулентной энергии. Из анализа размерностей к и
е следует, что турбулентная кинематическая вязкость
должна выражаться как:
к2
Ут=<у—>	(5.21)
£
111
Глава 5. Турбулентность: проблемы моделирования
где сц — некоторая константа. Экспериментально уста-
новлено, что сц = 0.09.
Уравнение для к является точным. Оно получается при
подстановке значения к в уравнения Навье—Стокса (2.2).
Уравнение для е также выводится из них с последующими
преобразованиями и гипотезами относительно величины
отдельных членов. Запишем уравнения модели без вывода:
дк — дк
— + U,----
дт 1 dxj
o-f = 1.3; cfl=1.44; ce2=1.92
Кроме к—е существует целый ряд других дифференци-
альных моделей. Например, семейство моделей к—(О.
Впервые модель такого вида предложена знаменитым
русским математиком Колмогоровым в 1942 году. Эта мо-
дель содержит уравнения переноса кинетической энергии
турбулентности к и удельной (в единице объема) скорости
112
5.6.	Модель к-е. Дифференциальные модели
TZ	2
диссипации энергии со. Иногда со определяют как осред-
ненный квадрат пульсаций завихренности. Ее размер-
ность — (время)-2. Она связывается с к и е следующим со-
отношением:
<5-24>
где cd — некоторая константа. Обычно полагается, что
cD = ср = 0.09 .
При этом далеко не все дифференциальные модели со-
стоят из двух уравнений. К примеру, распространенная в
расчетах дозвуковой аэродинамики модель Спалларта и
Алмареса содержит только одно уравнение, а модель Пола
Дурбина v — f— четыре уравнения. У каждой из моделей
есть свои преимущества и недостатки, свои требования к
сетке и особенности решения. Достаточно подробно об
этих особенностях излагается в [10], мы же несколько
позже приведем таблицу с кратким описанием наиболее
популярных моделей.
5.7.	О моделировании
течений вблизи стенки
Оказалось, что далеко не все модели турбулентности
адекватно отражают процессы, происходящие в присте-
ночной области. Причиной этого является тот факт, что
вблизи стенки силы вязкости часто преобладают над сила-
ми инерции, течение там не полностью турбулентно.
По признаку того, способна ли модель адекватно смо-
делировать пристеночный слой, все их делят на высоко-
рейнольдсовкие (High-Reynolds) и низкорейнольдсовские
(Low-Reynolds). Низкорейнольдсовкие модели без привле-
чения каких-либо дополнительных приемов позволяют
адекватно моделировать пристеночный слой (при этом им,
правда, требуется гораздо более качественная сетка). Вы-
сокорейнольдсовские модели не обеспечивают адекватно-
го моделирования пристеночной области течения.
Самым частым приемом для моделирования присте-
ночных слоев при помощи высокорейнольдсовских мо-
делей является использование метода пристеночных
функций.
114
5.7. О моделировании течений вблизи стенки
Известно, что пристеночная область течения может
быть разбита на три зоны (см. Рис. 5.9):
1)	Вязкий подслой, в котором вязкие напряжения до-
минируют над рейнольдсовыми и имеет место линейная
зависимость скорости потока от расстояния от стенки:
+	+	+ U	+ и.у
и = у , где м = — — безразмерная скорость, у =------
ит	V
безразмерное расстояние от стенки, ыг
— динами-
ческая скорость.
2)	Буферный слой, где вязкие и рейнольдсовы напря-
жения имеют один порядок. «Сшивая» профили скорости
для вязкого подслоя и логарифмического слоя, прибли-
женно получают:
и+ = 51пу+ +3.05.
(5.25)
3)	В логарифмическом слое рейнольдсовы напряжения
намного превышают вязкие, а профиль скорости может
быть представлен в форме логарифмического закона:
u+=-ln(£y+),	(5.26)
к v '
где к — постоянная Кармана, ^-постоянная, определяю-
щая степень шероховатости (для гладкой стенки экспери-
ментально установлено Е = 8.8).
115
Глава S. Турбулентность: проблемы моделирования
Рис. 5.9. Структура пристеночного слоя.
Описанные участки обычно объединяются в одну внут-
реннюю область, которая занимает порядка 20% толщины
пограничного слоя и в которой генерируется около 80%
всей энергии турбулентности. Одно из важных свойств
внутренней области заключается в том, что профиль ско-
рости слабо зависит от числа Рейнольдса, продольного
градиента давления и прочих внешних условий (которые,
тем не менее, могут вызвать уменьшение толщины внут-
ренней области или даже полное ее вырождение). Именно
это свойство послужило основой для построения универ-
сальных соотношений (пристеночных функций), связы-
116
5.7. О моделировании течений вблизи стенки
вающих параметры течения с расстоянием от стенки. На-
ряду с универсальностью профиля скорости во внутренней
области, метод пристеночных функций опирается на ис-
пользование гипотезы о локальном равновесии энергии
турбулентных пульсаций, а также свойства локальной изо-
тропности диссипирующих вихрей.
Следует также сделать замечание по поводу требова-
ний к сетке. Обычно они выражаются через безразмерное
расстояние от стенки у+ ближайшего узла. Для низкорей-
нольдсовских моделей у+ ~ 1, для высокорейнольдсов-
ских оно обычно находится в пределах от 30 до 300. При-
менение высокорейнольдсовских моделей за пределами
указанного диапазона может приводить к существенным
ошибкам в определении характеристик трения и теплооб-
мена на стенке.
К сожалению, точное значение у+ может быть получе-
но только после вычисления напряжений на стенке (т.е.
уже после решения задачи). Поэтому сетки, после первых
расчетов иногда приходится перестраивать заново, чтобы
_ _ +
подобрать требуемое значение расстояния^ .
Чтобы получить хотя бы примерные сведения о том, ка-
кого размера следует задавать ячейки у стенки, можно вос-
пользоваться прикидочными формулами. Например, из-
вестно, что для круглой трубы при стационарном течении
117
Глава 5. Турбулентность: проблемы моделирования
(5-27)
где £ — коэффициент трения, вычисляемый по эмпириче-
ским зависимостям (см. например [15]), I — длина канала, „
d—диаметр канала.
5.8.	Общие данные
о некоторых моделях турбулентности
Приведем данные о самых популярных и известных
моделях турбулентности:
Название модели	Тип	Low Re или High Re	Краткое описание
Модель Себеси—Смита (Cebeci-Smith)	алг.	Low Re	Двухслойная модель, разделяет по- ток на две области — внешний по- ток и пристеночный слой. Применя- ется для высокоскоростных безот- рывных потоков. Очень полезна для предварительных расчетов, когда не важна детальная физика процесса.
Модель Болдуина— Ломакса (Baldwin-Lomax)	алг.	Low Re	Сходна с моделью Себеси—Смита по своим характеристикам
Модель Спаларта— Аллмареса (Spalart- Allmaras, SA)	диф.	Low Re	Модель содержит одно дифферен- циальное уравнение. Создавалась для задач внешней дозвуковой аэ- родинамики. Опыт эксплуатации модели SA показал, что ее реальные возможности заметно шире, чем предполагалось при ее создании. Более того, после введения в нее поправок на кривизну линий тока и вращение, границы ее применимо- сти модели заметно расширились.
119
Глава 5. Турбулентность: проблемы моделирования
Окончание таблицы
Название модели	Тип	Low Re или High Re	Краткое описание
уг-92	диф.	Low Re	Эта модель содержит одно уравне- ние. Она обеспечивает вполне удовлетворительное описание не только большинства канонических сдвиговых течений (плоская и осе- симметричная струя, слои смеше- ния в несжимаемой и сжимаемой жидкости, пограничный слой на плоской пластине при отсутствии и при наличии шероховатости по- верхности и др.), но и ряда более сложных течений, представляющих практический интерес.
к-е	Диф.	High Re	Хорошо описывает полностью раз- витую турбулентность. Одна из самых популярных моделей, вклю- чена во все коммерческие пакеты по вычислительной гидродинамике. Имеет целый ряд модификаций, в том числе и низкорейнольдсовские модификации.
Модель к-а> Колмогорова	Диф.	High Re	Исторически самая первая модель с двумя дифференциальными урав- нениями. Не содержит членов, от- ражающих влияние молекулярной вязкости на турбулентность. Сейчас практически не применяется.
Модель к-а> Саффмена— Вилкокса	Диф.	Low Re	Хорошо описывает пристеночные течения. Хуже дело обстоит со сво- бодной развитой турбулентностью.
120
5.8. Общие данные о некоторых моделях турбулентности
Название модели	Тип	Low Re или High Re	Краткое описание
Модель перено- са касательных напряжений (Shear Stress Transport, SST) Ментера	диф.	Low Re	Является комбинацией моделей к-8 и к-а). Для пристеночного слоя используется к - о, для внеш- него региона — к - в. В настоящее время эта модель является очень популярной и входит во многие пакеты вычислительной гидроди- намики.
V2-/ модель	диф.	Low Re	Данная модель сходна со стандарт- ной к —8 моделью, однако в ней сделана попытка учесть неизотроп- ные эффекты, возникающие около стенки. Для этого добавлено урав- нение для квадрата пульсаций по- перечной скорости V2 (что факти- чески представляет собой пульса- цию линии тока) и специального релаксационного фактора / Таким образом, модель содержит четыре уравнения.
5.9.	Более сложные модели
турбулентности
Как упоминалось выше, модели, основанные на гипоте-
зе Буссинеска не всегда применимы. Поэтому было разра-
ботано несколько более продвинутых подходов к модели-
рованию турбулентных течений. Кратко рассмотрим самые
важные из них.
Модели переноса рейнольдсовых напряжений
Точные уравнения для каждого компонента тензора
(5.12) могут быть получены из уравнения (2.2). Затем вы-
двигаются гипотезы о возможном значении отдельных
членов в этих уравнениях (модели для переноса рейнольд-
совых напряжений отличаются друг от друга как раз этими
гипотезами). К шести уравнениям для рейнольдсовых на-
пряжений добавляется также уравнение для диссипации
турбулентной энергии вида (5.23). Таким образом, модель
состоит из семи дифференциальных уравнений.
Модели для рейнольдсовых напряжений рекомендуется
применять для моделирования сложных отрывных тече-
122
5.9. Более сложные модели турбулентности
ний, решения задач о взаимодействии струй с препятст-
виями, моделирования течений в каналах сложной формы.
DNS, LES и DES
Все подходы к моделированию турбулентности, опи-
санные выше, основаны на одном и том же принципе —
осреднении Рейнольдса. Если внимательно рассмотреть
суть математических операций данного осреднения, то
становится ясно, что в конечном итоге операция осредне-
ния по Рейнольдсу сводится к тому, что в уравнение дви-
жения в том или ином виде добавляется дополнительная
вязкость. Таким образом, задача о турбулентном течении
заменяется задачей о ламинарном (точнее будет сказать —
«квазиламинарном») течении. Главным недостатком тако-
го подхода является то, что модели, основанные на осред-
нении Рейнольдса не в состоянии воспроизвести многие
эффекты, возникающие при рождении, взаимодействии и
распаде вихревых структур.
Указанных недостатков лишен подход, называемый
прямым численным моделированием (DNS).
Как уже упоминалось выше, под прямым численным
моделированием понимается подход, когда для математи-
ческого моделирования турбулентных течений никаких
специальных моделей турбулентности не применяется, а
все эффекты, связанные с возникновением и эволюцией
123
Глава 5. Турбулентность: проблемы моделирования
турбулентных вихрей, получаются напрямую из решений
системы уравнений Навье—Стокса.
Прямое численное моделирование (DNS), начиная с
80-х годов, достаточно быстро прогрессирует, хотя дос-
тижимые расчетные числа Рейнольдса пока еще остаются
слишком низкими, чтобы интересовать практикующих
инженеров. К настоящему времени получены данные
DNS для ряда двумерных и трехмерных течений, в том
числе с отрывом потока, и список приложений продолжа-
ет расти [10].
Более простой моделью является моделирование
крупных вихрей (Large Eddy Simulation, LES). В этом
подходе крупные вихри рассчитываются, а мельчайшие
вихри подсеточного масштаба (Sub-Grid Scale, SGS) мо-
делируются (см. Рис. 5.10). Основной предпосылкой тако-
го подхода является то, что наибольшие вихри, которые
находятся под прямым воздействием граничных условий,
несут максимум рейнольдсовых напряжений и должны
быть рассчитаны.
Мелкомасштабная турбулентность является слабой, со-
держащей меньше рейнольдсовых напряжений. Также она
близка к изотропной и имеет характеристики, близкие к
универсальным. Поэтому она в большей мере поддается
упрощенному моделированию.
124
5.9. Более сложные модели турбулентности
U
-----DNS
-----LES
Рис. 5.10. Схематичное представление масштабных отличий между
LES и DNS (в подходе LES напрямую моделируются только круп-
ные вихревые структуры, а в подходе DNS — вихревые структуры
всех масштабов).
Основной прием, используемый в LES в чем-то сходен
с осреднением, использованным Рейнольдсом (формула
(5.2)). Любая величина представляется здесь в следующем
виде:
<р = ф + ф.	(5.28)
Здесь ф — отфильтрованная, крупномасштабная часть
величины (моделируется напрямую), ф — подсеточная
часть величины (моделируется упрощенно).
Основное различие между вариантами LES состоит в
разном подходе к моделированию подсеточных масшта-
бов. Подробнее о LES можно прочесть в [11].
Метод моделирования отошедших вихрей (Detached
Eddy Simulation, DES) является логическим развитием LES.
125
Глава 5. Турбулентность: проблемы моделирования
Суть метода заключается в том, что расчетная область раз-
деляется на слой около стенки и остальную область. При-
стеночный слой моделируется с помощью более простых
моделей турбулентности (например, одной из моделей с
турбулентной вязкостью), а вся остальная область — с по-
мощью LES.
5.10.	О применении моделей
турбулентности
В данном разделе мы познакомились с подходами к
численному моделированию турбулентных течений.
Как мы увидели, наиболее распространенным подхо-
дом является замена исходных уравнений, описывающих
течение на уравнения для осреднённых характеристик с
добавлением дополнительных («турбулентных») членов,
описывающих турбулентный перенос (вещества, количест-
ва движения, энергии). При этом опыт показывает, что эти
дополнительные члены часто существенно превосходят по
величине перенос, осуществляемый за счёт молекулярного
взаимодействия. Таким образом, турбулентный перенос
при превышении критического числа Рейнольдса стано-
вится одним из основных механизмов переноса вещества,
количества движения и энергии в потоке.
Определение характеристик турбулентного переноса
осуществляется с помощью моделей турбулентности. При
этом, как мы уже убедились, универсальных моделей тур-
булентности на сегодняшний день не существует. Для раз-
127
Глава 5. Турбулентность: проблемы моделирования
ного типа течений разные модели дают разные результаты,
которые могут быть как близки к наблюдаемым в действи-
тельности, так и нет.
В случае, когда результаты расчетов являются неудов-
летворительными (плохо отражают реальную картину те-
чения), проблема очень часто кроется именно в неудачном
выборе модели турбулентности.
Отсюда можно сделать вывод, что правильный выбор
модели турбулентности при решении задачи о турбулент-
ном течении является ключом к получению адекватных
результатов.
Глава 6
Особенности
расчёта химически
реагирующих течений

Очень часто течение газов и жидкостей при высоких
температурах сопровождается химическими реакциями.
Эти реакции могут быть самыми разнообразными — от
горения топлива и взаимодействия среды со стенкой, до
реакций диссоциации.
В данном разделе мы рассмотрим лишь гомогенные
химически реагирующие течения, не касаясь многофазных
(гетерогенных) систем.
6.1. Основные положения
Представим течение многокомпонентной газовой сме-
си, в которой происходят химические реакции. Отметим,
что, в отличие от многофазной смеси, все химические ком-
поненты смешаны на молекулярном уровне.
Введем основные параметры, характеризующие газо-
вую смесь.
Компоненты обычно обозначаются большими латин-
скими буквами: А, В, С и т.д. Nc — общее число компо-
нентов газовой смеси; А — массовая концентрация ком-
понента I, т.е. масса компонента I в единице объема,
размерность [кг/м3].
Очевидно, что сумма массовых концентраций всех
компонентов равна плотности газовой смеси:
А
Z А=А	(6.1)
1=Л,В,С,...
Массовая доля I-го компонента Cj определяется по
формуле:
132
6.1. Основные положения
C/=J-	(6-2)
Из (6.1) следует, что
[7] — мольная концентрация компонента I — это ко-
личество молей (киломолей) данного вещества в единице
объема. Чаще всего используется размерность [кмоль/м3].
Очевидно, что
= = (64)
где Wj —молекулярная масса компонента I [кг/кмоль].
Суммируем (6.4) по всем компонентам
(6-5)
и вводим понятие мольной доли компонента I:
(6.6)
Из формул (6.4) - (6.6) следует связь между мольной и
массовой долями компонента I
133
Глава 6. Расчёт химически реагирующих течений
/ Nc
Cj/W.
Обратную связь можно получить, если домножить
уравнение (6.4) на W[, потом просуммировать по всем ком-
понентам и разделить полученные уравнения друг на друга:
с А _	№
Известно, что объемная доля компонента I равна его
мольной доле.
6.2. Основные уравнения
Наряду с уравнениями системы (2.1) - (2.3) для много-
компонентной смеси решаются также уравнения нераз-
рывности (переноса) для каждого компонента. Это уравне-
ние для компонента I имеет вид (здесь используется тен-
зорная форма записи):
д(рс,) t аршА)
dt dxj
+ 51/,
(6.9)
(I) (II) (Ш) (IV)
где Гleff — эффективный коэффициент диффузии для
компонента I:
Гл<#=Г/+7“’	(6.10)
ЙС у»
где Г, = pDj — молекулярный коэффициент диффузии
компонента I;
Dj —коэффициент диффузии;
Рт — коэффициент турбулентной вязкости (о нём мы го-
ворили в предыдущем разделе);
135
Глава 6. Расчёт химически реагирующих течений
Sc
_ Нт
pDT
— турбулентное число Шмидта, обычно пола-
гают, что Scr = 0.9;
DT — коэффициент турбулентной диффузии;
S, — источник компонента I за счет всех химических
реакций, в которых он участвует (скорость образования
компонента I [кг /(м3 с)]).
Физический смысл членов, входящих в уравнение (6.9):
Изменение концентрации компонента I (член I) про-
исходит за счет конвективного подвода этого компонента
(Я), за счет диффузионного подвода (Ш) и за счет образо-
вания компонента в результате химических реакций (IV).
Вопросы решения уравнений, содержащих конвекцию
и диффузию, достаточно подробно освещены в предыду-
щих разделах. Поэтому основное внимание уделим рас-
смотрению в уравнении источникового члена и его влия-
нию на выбор способа решения системы уравнений (6.9).
Если в химически реагирующей системе протекают ре-
акции, число которых равно К, то скорость образования
компонента I складывается из скоростей всех реакций, в
которых он участвует.
Обычно химическую реакцию под номером к представ-
ляют в виде:
136
6.2. Основные уравнения
Nc	Nc
£ y'Jo £ v*/,	(6.U)
где Nc — число компонентов, — стехиометрический
коэффициент компонента I в к-ой прямой реакции (слева
направо), у£ — стехиометрический коэффициент компо-
нента I в к -ой обратной реакции (справа налево).
Скорость химической реакции Rk определяется как из-
менение мольной (молярной) концентрации одного из реа-
гирующих веществ за единицу времени (Т.е. размерность
Rk равна кмоль/(м3 с). Строго говоря, в этом определении
следует нормировать изменение мольной концентрации
реагирующего вещества на его стехиометрический коэф-
фициент).
С учетом (6.4) массовая скорость образования компо-
нента I, таким образом, определяется по формуле:
к
(6.12)
*=1
Проиллюстрируем эти определения на примере.
Рассмотрим взаимодействие смеси окиси углерода СО
и водорода Щ с кислородом Ог. Примитивное школьное
представление этого взаимодействия выглядит так:
Реакция 1: СО+-^О2 -»СО2,	(6.13)
137
Глава 6. Расчёт химически реагирующих течений
Реакция 2: Н2 +|о2 Н20.	(6.14)
В этом случае скорость образования СОг равна *$с02 =
= V, скорость образования СО равна Sco = -ffcoRi и
имеет отрицательный знак, т.е. компонент СО2 расходует-
ся. А вот кислород расходуется в двух операциях и его
скорость образования равна
«о. -	-I*,] =	(Я. +«,).
\ Z Z J L
Вернемся к скоростям реакций.
В 1865 году Н. Н. Бекетовым и в 1867 году Гульдбер-
гом и Вааге был сформулирован закон действующих масс:
Скорость химической реакции в каждый момент
времени пропорциональна концентрациям реагентов,
возведенным в некоторые степени.
Для элементарных реакций показатель степени при
значении концентрации каждого вещества равен его сте-
хиометрическому коэффициенту. Следует отметить, что
для более сложных реакций это правило не соблюдается.
( Nr	Nc	\
я* = Fk П -Вк п
\	I~A,B,...	I~A,B,...	J
(6.15)
138
6.2. Основные уравнения
где Fk и Вк — коэффициенты скоростей прямой и обратной
А?ой реакции, соответственно (размерность [(м3/кмоль)"-1/с]),
п — порядок реакции. Например, для реакции
Н+ОН+М = Н2О+М
(6.16)
порядок прямой реакции равен 3, а обратной — 2. Ско-
рость этой реакции определяется по формуле:
/? = F[H][0H][M]-5[H20][M]	(6.17)
Здесь под М подразумевается любой химический ком-
понент газовой смеси. Этот компонент выступает в роли
катализатора реакции, т.к. он не изменяется в процессе ре-
акции, но обеспечивает ее протекание. Для такого компо-
нента существует специальный термин: «третье тело».
Коэффициенты скоростей реакций рассчитываются по
формулам:
F,=4TAexp(-^),	(6.18)
В, = 4Гаехр^-^,	(6.19)
где Ак — предэкспоненциальный множитель; Д — без-
размерный температурный показатель степени; Тк — тем-
пература активации; Т —текущая температура.
139
Глава б. Расчёт химически реагирующих течений
Константы Ак, 0к, Тк имеют разное значение для пря-
мой и обратной реакций.
Существует два способа задания этих констант:
1) задаются все константы для прямой и обратной ре-
акции,
2) задаются константы для прямой реакции, а скорость
обратной реакции определяется по формуле:
B = FIKC	(6.20)
(индекс «£» для упрощения записи опущен, но подразуме-
вается). Константа равновесия Кс определяется с помощью
свободной энергии Гиббса:
( р V Г ЛЬ	_ -|
кс = ехр	,	(6.21)
Nc
где:	•' = £(>'.'Р.=10![Па];
i=l
Rv = 8314.41 [Дж/кмольхК]; Rs — газовая постоянная
компонентам; gs=hs-ssT —энергия Гиббса.
Для определения теплоемкостей, энтальпий и энтропий
компонентов обычно используются полиномы вида:
^2 = а1+а2Т+а3Т2 + а4Г3+а5Т4,	(6.22)
R
h\T) Т Т2 Т3 Т* Ь
-^ = *1+*2у + «зу + «4-^+а5Т+т’	6'23)
140
6.2. Основные уравнения
^-^ = aIln(T)+a2T+
, а д	(624)
+а,— + ал — + а,—+Ь,,
3 2	4 3	5 4	2
где СЦТ) — теплоемкость компонента при постоянном
давлении при температуре Т; R — газовая постоянная ком-
понента,
h° (Г) = {й° (Т) - Л° (0)} + h° (0) =	(6.25)
= {Л"(Г)-»’(Г,)}+Л’(Г0),
Т0 =298.15К —стандартная температура; Л°(0) —хими-
ческая энергия компонента при 0 К; h°(T0) — энтальпия
компонента при стандартной температуре То, которую
принимают равной
A°(r0) = Az№(T0),	(6.26)
гдеАу.№(Т) — энтальпия образования (теплота образова-
ния) компонента из элементарных веществ при температу-
ре Т и стандартных условиях. В этих формулах индекс «у»,
относящийся к номеру компонента, для простоты записи
опущен.
При использовании для определения термодинамиче-
ских свойств полиномов вида (6.22) - (6.24) энергия Гиб-
бса задается формулой
141
Глава 6. Расчёт химически реагирующих течений
— = -^—— = а, Г1-1п(Т)"1-а2 —-
R,T RST R, ML V 7J u 2	(6.27)
T2 T3 T4 bls ,
-a,,----a.,-----a., — ч—--D,,
3,. g	4,j	V 20 T 3,2
Упрощенная формула получается без учёта изменения
теплоёмкости реакции (так называемое первое приближе-
ние Улиха):
де° = ДЯ^-ТД5®и,	(6.28)
ехр
ехр
-ая;°98
. RVT
(6.29)
*4»
Выбор набора элементарных реакций зависит от кон-
кретной задачи.
А что же такое элементарные реакции?
Элементарные реакции — это химические реакции, ко-
торые не могут быть представлены более простыми хими-
ческими превращениями. Элементарные реакции — со-
ставные части сложной реакции. Иногда вместо термина
«элементарная реакция» пользуются терминами «элемен-
тарная стадия» или просто «стадия» (сложной реакции).
Дело в том, что реакции, представленные формулами
(6.13) и (6.14), являются сложными реакциями, и скорость
их протекания аналитически определить практически не-
возможно. Для математического описания процесса горе-
ния СО и Нз в кислороде необходимо представить этот
142
6.2. Основные уравнения
процесс в виде системы элементарных реакций. Возможен
набор различных систем элементарных реакций (цепочек
реакций). Один из самых простых представлен в Табл. 6.1.
Эти реакции могут идти как вправо, так и влево. Как
уже указывалось, М в реакциях 5-8 и 10 — это так назы-
ваемое третье тело; им может быть любой химический
компонент, входящий в смесь. Третье тело не изменяется в
результате реакции, оно только способствует распаду мо-
лекулы при соударении с ним (реакции диссоциации), либо
соединению атомов или молекул в более сложные молеку-
лы (реакции рекомбинации).
Табл. 6.1. Элементарные реакции процесса
горения СО и Hj в кислороде
Цепные реакции	
Реакция 1	Н+О2=ОН+О
Реакция!	О+Н2=ОН+Н
Реакция 3	ОН+Н2=Н2О+Н
Реакция 4	ОН+ОН=Н2О+О
Реакции рекомбинации — диссоциации	
Реакция 5	н+н+м=н2+м
Реакция 6	Н+ОН+М=Н2О+М
Реакция?	н+о+м=он+м
Реакция 8	о+о+м=о2+м
Реакции с углеродом	
Реакция 9	со+он=со2+н
Реакция 10	со+о+м=со2+м
143
Глава 6. Расчёт химически реагирующих течений
Для оценки влияния химических реакций имеет смысл
ввести так называемое число Дамкелера (Damkohler
number), которое для каждой реакции определяется как
_ скорость химической реакции
скорость конвективного переноса	<Х».ЗО)
или
Da = ^~,	(6.31)
^chan
где ^gas — характерное газодинамическое время (масштаб
времени переноса вещества за счет конвекции); гсАеот —
характерное время протекания химической реакции.
Если рассматривать все реакции, то число Дамкелера
это отношение члена (IV) в уравнении (6.9) к члену (II).
Рассмотрим более подробно, как происходит окисление
водорода в системе реакций, представленной в таблице 6.1.
Предположим в системе появился 1 (один) атом сво-
бодного радикала О.
Один из механизмов следующий:
радикал О взаимодействует с молекулой Нг через реак-
цию 2, и образуются один радикал ОН и один радикал Н;
радикал Н взаимодействует с молекулой Ог через реак-
цию 1, и образуются один радикал ОН и один радикал О
(т.е. имеется уже два радикала ОН);
144
6.2. Основные уравнения
два радикала ОН взаимодействуют между собой через
реакцию 4 , и образуются одна молекула воды НгО и один
радикал О.
Общий результат такой цепочки:
О+Н2+О2 -> Н2О+2О,	(6.32)
т.е. цепочка замкнулась, вместо одного радикала О полу-
чилось 2 радикала, молекулы Щ и Ог превратились в моле-
кулу воды НгО.
Можно рассмотреть и другие варианты цепочек.
Основной эффект реакций 1-4 это лавинообразное уве-
личение числа свободных радикалов о, Н, ОН и превраще-
ние Нг и Ог в НгО, т.е. горение водорода в кислороде.
Реакции 1-4 называются цепными. ЦЕПНЫЕ РЕАК-
ЦИИ — химические реакции, идущие путем последова-
тельности одних и тех же элементарных стадий, на каждой
из которых возникает одна или несколько активных частиц
(атомов, свободных радикалов, ионов, ион-радикалов).
В частности, по цепному механизму протекают реакции
горения.
Эти реакции протекают с очень большой скоростью и
считаются быстрыми. Для них TChem очень мало.
Отток свободных радикалов осуществляется за счет
медленных реакций рекомбинации 5-8. Для них характер-
145
Глава 6. Расчёт химически реагирующих течений
ное время протекания химических реакций tchem значи-
тельно больше, чем для реакций 1-4.
Резонно задать вопрос, откуда в смеси На и Ог появля-
ются первые свободные радикалы?
Они могут возникать в результате реакций диссоциа-
ции (реакции 4-8, идущие справа налево). Реакции диссо-
циации протекают только при достаточно высоких темпе-
ратурах.
Таким образом, для воспламенения смеси На и Оа необ-
ходима некоторая начальная вспышка, которая создаст
достаточное количество свободных радикалов для ини-
циирования цепных реакций горения.
Справка из энциклопедии. Свободные радикалы в
химии— частицы (как правило, неустойчивые), содержа-
щие один или несколько неспаренных электронов, оксидан-
ты, в целом частицы (или интермедиаты) электронейтраль-
ны. По другому определению свободный радикал— вид
молекулы или атома, способный к независимому существо-
ванию (то есть обладающий относительной стабильностью)
и имеющий один или два неспаренных электрона.
Как было уже указано, характерные времена протека-
ния химических реакций тскет для цепных реакций и для
реакций рекомбинации-диссоциации очень сильно отли-
чаются друг от друга. Т.е. справедливо:
146
6.2. Основные уравнения
Da .. > Da .,
chain	recomb*
(6.33)
Dachain — число Дамкелера для цепных реакций, Darecomb —
число Дамкелера для реакций рекомбинации-диссоциации.
Учитывать необходимо и те, и другие реакции. Это по-
рождает серьезную математическую проблему, возникаю-
щую при решении системы уравнений переноса (6.9), т.к.
эта система является жесткой.
В следующем пункте мы рассмотрим понятие жестко-
сти систем уравнений.
6.3. Жёсткие системы
Рассмотрим пример (пример 1). Пусть имеется система
уравнений
(6.34)
с начальными условиями
y,(0) = y2(0) = l.
(6.35)
Зададим значения числовых коэффициентов равны-
ми: Л, = -1 и	= —10б
Система легко решается аналитически, и ее решение
имеет вид:
У1 (г) = ет; у2 (г) = е"10‘т.	(6.36)
Нам интересно проверить решение системы численны-
ми методами, поэтому используем для решения системы
(6.34) явный метод Эйлера.
Из теории численных методов известно, что для сходи-
мости решения конечно-разностных уравнений к точному
148
6.3. Жёсткие системы
решению необходима так называемая устойчивость разно-
стной схемы.
Метод Эйлера применительно к системе (6.34) дает
следующее численное решение:
т	(6.37)
где
i = (1,2), т — шаг конечно-разностной сетки,
у" — значение У,- на п -ом шаге по т ,
— значение на (п +1) -ом шаге по т.
Конечно-разностная сетка задается соотношением:
тп=пт.	(6.38)
Таким образом,
у! = (1+Ц),	= (1 + Ц)’..у; = (1+ гД)".	(6.39)
Условием устойчивости разностной схемы (6.37) явля-
ется:
|Дт+1|<1, |Лг+1|<1;	(6.40)
т.е. т<2-10~6.
Предположим, что нам необходимо определить реше-
ние системы при t = 3. Точное решение (25) равно:
149
Глава 6. Расчёт химически реагирующих течений
yt (3) = 0,049787 , у2 (3) = е 310 . Значение у2 выходит за
рамки разрядности любой вычислительной машины и фак-
тически равно нулю.
Решение конечно-разностного уравнения (6.37) с уче-
том (6.38) и (6.39) имеет вид:
Л(«) = (1+А’-)<"')-	<6.41)
При шаге г = 5 • 10~7 и f = 3 получаем:
у, (3) = (1 -10"6	= 0.04978, у2 (3) = 0.
Для более раннего момента времени t = 10^ при том же
шаге г = 510-7 получаем:
точное решение —
у, (10*4) = 0,9999, у2 (10-4) = 3.72 10-44;
численное решение —
у, (10^) = 0,9999, у2 (Ю-4) = 6.22 Ю^1.
Мы видим, что для У\ выбранный шаг дает хорошее
совпадение для обоих моментов времени; для у2 относи-
тельное совпадение плохое, но самом деле это не имеет
большего значения, т.к. в обоих случаях решение фактиче-
ски равно нулю.
150
6.3. Жёсткие системы
Для получения более точного численного решения вто-
рого уравнения необходимо еще больше уменьшить шаг
интегрирования. Например, если т = 110-’, то численное
решение при t = 10-4 гораздо ближе к точному решению:
у2(10'4) = 3.54-10'44.
С другой стороны для получения хорошей точности
при решении первого уравнения совсем необязательно за-
давать такой маленький шаг как т = 5 10~7. При шаге в
20000 раз больше, т.е. при т = 1102 получаем численное
решение: у, (3) = 0,049, у, (10-4) = 0,9999 — очень близкое
к точному решению.
Систематизируем полученные данные в Табл. 6.2 и
Табл. 6.3.
Табл. 6.2. Сравнение точного и приближенного решения
первого уравнения при различных шагах по времени
Время т = 3				
Шаг т	У1 точное решение	У1 численное решение	У2 точное решение	У2 численное решение
5-Ю"7	0,049787	0,049787	0	0
110“2		0,04904 (погрешность 1.5%)		решение невозможно
110^		0,049787		0
151
Глава 6. Расчёт химически реагирующих течений
Табл. 6.3. Сравнение точного и приближенного решения
второго уравнения при различных шагах по времени
Время / = !•]			Ю"4	
Шаг т	У1 точное решение	У1 численное решение	У1 точное решение	У1 численное решение
5 10"7	0,9999	0,9999	3.72-10"**	6.22-10"61
110"2		0,9999		решение невозможно
110"’		0,9999		3.54 10м4
Основной вывод, который можно сделать, исходя из
полученных результатов, состоит в том, что шаги Т опти-
мальные для первого и второго уравнения отличаются друг
от друга на несколько порядков. Это приводит к основа-
тельным проблемам. Решать первое уравнение с шагом
110"’ непозволительная роскошь с точки зрения расходо-
вания вычислительных ресурсов, а использовать шаг 1 • 10~2
для второго уравнения невозможно с точки зрения устой-
чивости.
Вообще говоря, трудности численного решения подоб-
ных систем, получивших название жестких, связаны с вы-
бором шага интегрирования. Дело в том, что характерные
времена исследуемых процессов могут различаться в 106
раз. Следовательно, при численном решении системы не-
обходимо выбирать шаг по самому быстрому процессу.
152
6.3. Жёсткие системы
И в таком случае затраты машинного времени для иссле-
дования самых медленных процессов будут неоправданно
велики.
Для решения рассматриваемых задач имеются следую-
щие подходы:
1.	Численно решать систему с шагом, выбранным по
условию (6.40), т.е. с учетом характерных времен всех
процессов, описываемых данной системой.
2.	Решать систему ОДУ с различными шагами, соответ-
ствующими физическим процессам, имеющим существен-
но различные характерные времена. В этом случае необхо-
димо задавать условия перехода к другому шагу по времени
при интегрировании.
3.	«Пренебречь» быстропротекающими процессами и
численно рассматривать лишь медленные, проводя интег-
рирование с шагом, превышающим характерные времена
быстрых процессов. В этом случае придется конструиро-
вать численные методы, позволяющие проводить расчеты
с шагом большим, чем следует из условия (6.40).
Недостатки подхода 1 мы уже показали.
Второй подход приемлем для только что рассмотрен-
ной задачи, когда правая часть первого уравнения зависит
только от первой функции, а правая часть второго — толь-
ко от второй. В общем случае в правой части каждого
уравнения могут быть все искомые функции.
153
Глава 6. Расчёт химически реагирующих течений
Систему уравнений (6.34) и начальных условий (6.35)
удобно представить в матричном виде:
^- = Ау, y(Q) = ye;	(6.42)
где У~
— вектор искомых функции, А — матрица
\Уы)
коэффициентов, У о — начальное значение вектора У,
N — число искомых функций.
В рассмотренном примере матрица коэффициентов
имеет диагональный вид
г-1	0
-io6,
(6.43)
и каждое уравнение системы можно решать отдельно.
В общем случае необходимо решать всю систему со-
вместно. Поэтому нам представляется наиболее оптималь-
ным использовать третий подход.
Но сначала давайте подробнее исследуем понятие же-
сткости систем обыкновенных дифференциальных урав-
нений.
Описанная выше ситуация возникает из-за большого
разброса собственных значений матрицы системы (6.42):
154
6.3. Жёсткие системы
= Ю6. Компонента с большим (по модулю) собствен-
ным значением вынуждает выбирать мелкий шаг и, одно-
временно, быстро перестает влиять на решение. Класс
дифференциальных уравнений с таким поведением выде-
ляется в теории численных методов понятием жестких
уравнений.
Точнее, система линейных автономных дифференци-
альных уравнений (6.42) называется жесткой, если, во-пер-
вых, все собственные значения Л. матрицы А имеют отри-
цательную вещественную часть (т. е. система (6.42) экспо-
ненциально устойчива).
(6.44)
и, во-вторых,
max|Re(/L)
£ _ 1<|<У । v '7
min|Re(A)
1«<лН v 17
(6.45)
Число S при этом называют коэффициентом жесткости
системы (6.42).
Значок » («значительно превосходит») на практике
обычно означает, что S > 100, хотя в химической кинетике
часто встречаются задачи с коэффициентом жесткости
» ю6 и более.
155
Глава 6. Расчёт химически реагирующих течений
Более подробно с понятиями жёсткости и устойчивости
систем дифференциальных уравнений можно ознакомить-
ся в [6,23].
Попробуем применить к системе (6.42) неявный метод
Эйлера.
уЛ + 1 - уЛ
----^ = Л/+1,	(6.46)
т
т.е.
=(E-.4r)V,	(6.47)
где Е —единичная матрица.
Отсюда получаем:
у2 =((£-Лт)’’У,
v 7	(6.48)
Уо> •••>
уп =((Е~АтУ1)п Уо-
дыдущем примере, матрица Е Ат -
Пример 2. Для системы (6.34), рассмотренной в пре-
+ т О
О	1 + 10бг
, а
обратная ей:
156
6.3. Жёсткие системы
О
1
1 + 1Обт>
(6.49)
Тогда решение на и-ом шаге интегрирования получает-
ся равным:
(6.50)
Если взять т = 1 • 10~2 (это максимальное значение ша-
га, использованного в предыдущем примере), то получим
результаты, представленные в Табл. 6.4.
Табл. 6.4. Сравнение численного решения
неявным методом Эйлера с аналитическим
Время г = 3				
Шаг Т	У1 точное решение	У1 численное решение	Уг точное решение	У1 численное решение
110’2	0,049787	0,05053 (погрешность 1.5%)	0	0
Основным преимуществом, полученным в результате
использования неявного метода Эйлера, является возмож-
157
Глава 6. Расчёт химически реагирующих течений
ность совместно решать всю систему с шагом интегриро-
вания намного превосходящим шаг, полученным из крите-
рия устойчивости (6.40). При этом появляется возмож-
ность ориентироваться только на медленные процессы,
проводя интегрирование с шагом, превышающим харак-
терные времена быстрых процессов.
Рассмотрим еще один пример, менее жесткий, в кото-
ром вторая функция не столь быстро уходит в нуль.
Пример 3. Пусть матрица коэффициентов и начальный
вектор равны
<-1	99 "
0 -100>
(6.51)
, Уо~
Собственные значения матрицы равны (-1, -100). Они
имеют отрицательную вещественную часть, т.е. решение
экспоненциально устойчиво. Коэффициент жесткости сис-
темы равен 100.
Точное решение системы:
У =
^2е"г + 100еЧООг
(6.52)
Результаты расчета для 2-х моментов времени пред-
ставлены в таблицах 6.5 и 6.6.
158
6.3. Жёсткие системы
Табл. 6.5. Сравнение численного решения
неявным методом Эйлера с аналитическим
Время г = 0.1				
Шаг т	Ух точное решение	Ух численное решение	Уг точное решение	Уг численное решение
0.01	1.814215	1.81 (погрешность 0.25%)	4.54 Ю’5	9.76-Ю"4
0.1		1.73 (погрешность 5%)		9.1-102
Табл. 6.6. Сравнение численного решения
неявным методом Эйлера с аналитическим
Время т = 1				
Шаг Т	Ух точное решение	Ух численное решение	Уг точное решение	Уг численное решение
0.01	0.7357589	0.739 (погрешность 0.5%)	3.78-10"14	7.9 Ю~31
0.1		0.771 (погрешность 5%)		3.86 Ю-11
Таким образом, мы получили неплохую точность для
медленных процессов.
6.4. Решение жёстких
систем уравнений
В предыдущем параграфе было введено понятие жест-
кости систем линейных обыкновенных дифференциальных
уравнений (ОДУ) и много внимания было уделено спосо-
бам их решения.
Спрашивается, какое это имеет отношение к расчету
химически реагирующих течений? Самое прямое. Дело в
том, что система уравнений неразрывности химических
компонентов (6.9), во-первых, также является жесткой, а,
во-вторых, для ее решения применимы все те идеи и мето-
ды, которые были изложены в п. 6.3.
Рассмотрим для начала систему (6.9) без учета конвек-
ции и диффузии:
— Oj.	(6.53)
dt
В векторной форме это уравнение можно представить
в виде:
/(<")’	(654)
160
6.4. Решение жёстких систем уравнений
где С
С2
— вектор массовых долей компонентов;
С,
%/р '
S2/p
^Nc / Р,
— вектор источников.
К уравнению (6.54) можно применить те же методы,
что и к линейным ОДУ, только вместо матрицы постоян-
ных коэффициентов необходимо использовать матрицу
Якоби:

Л.
(6.55)
VNc
dYNJ
Неявный метод Эйлера применительно к системе
(6.54) имеет вид:
-J1+1 __ п
-----У- = /" + л(у+1 - у).	(6.56)
161
Глава 6. Расчёт химически реагирующих течений
Откуда
(Е-тАп\У 1~уЯ =fn.	(6.57)
\ > т
Решение (6.57) основано на следующих действиях, вы-
полняемых на каждом n-ом шаге интегрирования:
1.	Вычисляется матрица производных (6.55) в точке
/•
2.	Следующая точка у"+1 находится из матричного
уравнения (6.57).
Приведем классический пример — систему Робертсона.
Рассмотрим систему трех уравнений:
dyx
-TL = -aiy1+a2y2y3,
ат
'	- «2У2 Уз - азУ2 >	(6.58)
ат
Начальные условия:
?1(0) = 1, л(0) = 0, Л(о) = о.
Система (6.58) представляет модель химического взаи-
модействия трех веществ: вещество «1» медленно превра-
щается в вещество «2»: «1» -> «2» (со скоростью ах = 0.1),
162
6.4. Решение жёстких систем уравнений
вещество «2» превращается очень быстро в вещество «3»:
«2»	«3» (а3 = 103).
И, наконец, вещество «2» при каталитическом воздей-
ствии вещества «3», превращается в вещество «1»
(а2 =102): «2» + «3» -> «1» + «3»
Матрица Якоби имеет вид:
«1
О
«2^3
~а2Уз~аЗ
аз
а1У2
~а2У2
о ,
(6.59)
Рис. 6.1. Результат решения системы (6.S8).
6.5. Метод расщепления для системы
уравнений переноса химических
компонентов
Основную систему (6.9) можно представить в вектор-
ной форме
^ = £(С)+/(С),	(6.60)
ОТ
тце L(C) — оператор, учитывающий диффузию и кон-
векцию:
£(<?) =
_а_Г 8С) a(/>“jC)
dxj kff dxj} dxj
(6.61)
/(С) — вектор источников (см. п. 6.4).
В системе (6.60) лучше отдельно рассматривать про-
цессы переноса (конвекция и диффузия) и процесс образо-
вания компонента в результате химических реакций. По-
этому для ее решения имеет смысл применить метод рас-
щепления по физическим процессам. Он состоит в сле-
дующем. На каждом шаге по времени вместо системы
уравнений (6.60) решается последовательность уравнений:
164
6.5. Метод расщепления
при выполнении условий:
(6.62)
(6.^3)
(6.64)
В работе [24] показано, что для сходимости такой сис-
темы (т.е. аппроксимации и устойчивости) достаточно,
чтобы сходилась каждая из систем (6.62) и (6.63).
Схема расщепления позволяет существенно сократить
вычислительные ресурсы.
В каждое уравнение системы (6.62) входит только один
химический компонент; таким образом, можно решать ка-
ждое уравнение отдельно. При численном решении это
приводит к системе, состоящей из Nnodes уравнений (Nnodes —
число узлов сетки).
Уравнение (6.63) решается для каждого узла сетки для
всех компонентов. Число неизвестных в нем равно Nc.
Если же решать сразу всю систему (6.60), то число не-
известных будет равно *NC- Решение таких систем
представляет значительную сложность с точки зрения вы-
числительных ресурсов.

Глава 7
Особенности расчёта
многофазных течений

Многофазные потоки состоят из двух и более отдель-
ных частей или «фаз», которые могут быть жидкостями
или твердыми частицами. Каждой из фаз присущи свои
характерные свойства.
В предыдущих разделах уже рассматривались много-
компонентные потоки (или газовые смеси). В чем отличие
многофазных и многокомпонентных потоков? В много-
компонентных смесях химические компоненты смешаны
на молекулярном уровне, т.е. различные микрочастицы
могут находиться в непосредственной близости. Много-
фазные течения, как правило, характеризуются достаточ-
но большими скоплениями однородных компонентов,
размер которых на несколько порядков превышает разме-
ры молекул.
Чаще всего мы можем разделить многофазные потоки
на 2 принципиально отличающиеся сущности: непрерыв-
ная среда и дискретная среда. Непрерывной средой в дан-
ном случае называют такую субстанцию, в которой мож-
но переместиться из одной точки в любую другую, не по-
кидая данную субстанцию. В дискретной среде это не-
возможно.
С этой точки зрения эти потоки можно разбить на не-
сколько классов.
169
Глава 7. Особенности расчёта многофазных течений
1.	Непрерывная среда — газ, дискретная среда — твер-
дые или жидкие частицы.
Пример — капли или кристаллики воды в непрерывной
среде воздуха. Для сравнения: вода почти всегда может
находиться в воздухе как газообразная фаза, и в этом слу-
чае мы имеем дело с многокомпонентной газовой смесью,
а не с многофазной средой.
Еще один пример — течение частиц AI2O3 в продуктах
сгорания твердотопливного ракетного двигателя. При этом
AI2O3 может находиться как в твердом, так и в жидком со-
стоянии.
2.	Непрерывная среда — несжимаемая жидкость, дис-
кретная среда — твердые частицы.
Пример — гидротранспорт: перенос плотно распреде-
ленных твердых частиц в потоке жидкости.
3.	Непрерывная среда — несжимаемая жидкость, дис-
кретная среда — газообразные частицы или частицы дру-
гой жидкости (случай, когда жидкости не растворяются
друг в друге).
Пример первого варианта — пузырьки газа в жидкости,
возникающие, в частности, при начале кипения. Второй
случай — масляные или бензиновые капли в воде.
Существуют и более сложные многофазные системы,
например, когда присутствуют 3 и более фаз. Примеры та-
ких систем будут рассмотрены позже.
170
Глава 7. Особенности расчёта многофазных течений
Большой интерес представляет случай, когда непре-
рывная и дискретная среды как бы меняются ролями. На-
пример, при течении воды в нагретом канале непрерывной
средой вначале является жидкость, а дискретной — пу-
зырьки пара; далее, по мере вскипания и испарения все
большей части воды уже пар становится непрерывной сре-
дой, а капли воды — дискретной.
По сравнению с теорией однофазной жидкости, (мо-
дель которой основана, по сути, на решении уравнений
Навье—Стокса) описание динамики многофазных пото-
ков представляет собой более сложную проблему. И хотя
довольно просто получить уравнения сохранения массы,
импульса и энергии для произвольной смеси, общий ана-
лог уравнений Навье—Стокса для многофазных течений
не найден. Тем не менее, используя правильную процеду-
ру усреднения, вполне возможно получить набор «урав-
нений для многофазных потоков», которые в принципе
достаточно адекватно описывают динамику любой мно-
гофазной системы и подчиняются только очень общим
допущениям. Недостатком такого подхода является усло-
вие, что эта система уравнений включает в себя большое
количество новых неизвестных переменных, и поэтому
для ее решения необходимо введение дополнительных
соотношений и применение дополнительных материаль-
ных законов.
171
Глава 7. Особенности расчёта многофазных течений
При моделировании многофазных потоков использу-
ются различные подходы. «Гомогенные модели» или «од-
нородные модели» являются простейшими. В них предпо-
лагается, что скорость и температура дискретных частиц в
каждой точке потока совпадают с соответствующими ло-
кальными свойствами непрерывной среды. В этом случае
описание течения сводится к решению уравнений Навье—
Стокса для непрерывной среды. Однако на практике ба-
ланс количества движения и энергии между разными фа-
зами наблюдается довольно редко.
Более обоснованные физические модели многофазных
потоков были разработаны, главным образом, с использо-
ванием двух различных подходов. В рамках подхода «Эй-
лер—Эйлер» (или сокращенно просто «метода Эйлера»)
все фазы формально рассматриваются как жидкости, кото-
рые подчиняются нормальным однофазным уравнениям
движения. Различные фазы обрабатываются математиче-
ски как взаимопроникающие континуумы. Поскольку объ-
ем фазы не может быть занят другими фазами, вводится
понятие фазовой объемной доли. Эти объемные доли счи-
таются непрерывными функциями пространства и време-
ни, а их сумма равна единице. Для каждой фазы получает-
ся набор уравнений сохранения, и эти уравнения имеют
сходную структуру для всех фаз. В уравнения входят но-
вые неизвестные величины, описывающие межфазное
172
Глава?.
взаимодействие, поэтому для замывания системы исполь-
зуются дополнительные соотношения, получаемые на ос-
нове экспериментальных данных или молекулярно-
кинетической теории.
Преимуществом метода Эйлера является его общность
и универсальность: в принципе он может быть применен к
любой многофазной системе, независимо от количества и
характера фаз. Недостатком метода является то, что он
достаточно часто приводит необходимости использовать
очень сложный набор основных уравнений и замыкающих
соотношений.
Другим распространенным подходом к описанию мно-
гофазных потоков является так называемая модель «Эй-
лер—Лагранж» (или сокращенно просто «метод Лагран-
жа»), который основан на анализе траекторий частиц.
В рамках этого подхода только жидкая фаза рассматрива-
ется как непрерывная, а движение дисперсной фазы частиц
получается путем интегрирования уравнения движения от-
дельных частиц вдоль их траекторий (в практических при-
ложениях «частица» может представлять собой одну физи-
ческую частицу или группу частиц). Жидкая фаза рассмат-
ривается как непрерывная, и для нее решаются уравнения
Навье—Стокса, а решение для дисперсной фазы находится
путем отслеживания большого количества частиц, пузырь-
ков или капелек в предварительно вычисленном поле тече-
173
Глава 7. Особенности расчёта многофазных течений
ния. При этом дисперсная фаза может обменивать импуль-
сом, массой и энергией с жидкой фазой.
Этот подход значительно упрощается, когда можно
пренебречь взаимодействиями частиц с частицами, и для
этого требуется, чтобы дисперсная фаза занимала незначи-
тельную объемную долю, даже если ее массовая доля дос-
таточно велика. Траектории частиц или капель определя-
ются индивидуально через определенные промежутки вре-
мени в процессе вычислении непрерывной фазы. Данная
модель оказывается весьма подходящей для моделирования
распылительных сушилок, сжигания угля и жидкого топли-
ва, а также и некоторых потоков с частицами, но она не го-
дится для моделирования жидкостно-жидкостных смесей,
псевдоожиженных слоев или любого другого случая, когда
объемной долей дисперсной фазы нельзя пренебрегать.
Рассмотрим отдельно эти два основных подхода.
7.1. Подход «Эйлер—Эйлер»
Предполагается, что для каждой фазы а справедливы
обычные уравнения механики жидкости и газа (МЖГ):
уравнения неразрывности, количества движения и энергии,
записанные в форме Эйлера:
Ра+^-(РЛа) = 0>	(7.1)
дт
) + V • (ЛЧЛг ) = -VPa + V •«« + F, ,	(7-2)
ot
аг " а' к а’	(7.з>
= V • (и„ • о„ -и„р„ -q„) + и  F + W„.
Здесь: ра — эффективная плотность фазы а определя-
ется как ра = кар, а, где ка — объемная доля фазы a, ps a —
физическая плотность фазы а, т.е. плотность самого мате-
риала, из которого состоит эта фаза; ио, ра, Еа — соот-
ветственно вектор скорости, давление, полная энергия фа-
зы а; оа, qa — соответственно тензор вязких напряжений
175
Глава 7. Особенности расчёта многофазных течений
и удельный тепловой поток фазы a; Fa — объемная внеш-
няя сила, в частности, сила межфазного взаимодействия,
т.е. сила, с которой остальные фазы воздействуют на фазу
a; Wa — объемное тепловыделение, в частности, то тепло,
которое фаза а получает от остальных фаз.
В уравнениях (7.1) - (7.3) мы для простоты пренебрегли
массообменом между фазами. Позже межфазному массо-
обмену будет посвящен специальный раздел.
Как уже говорилось, основную сложность описания
многофазных течений представляют дополнительные чле-
ны, появляющиеся в основных уравнениях. Они представ-
лены правыми частями в уравнениях (7.2), (7.3).
Для непрерывной среды используются обычные соот-
ношения, такие же, как для однофазной среды. А вот для
дискретной среды и для межфазного взаимодействия необ-
ходимо использовать дополнительные формулы.
Допущение: членами, содержащими давление ра, на-
пряжение трения оа и тепловой поток qa, в уравнениях
(7.2), (7.3) для дискретной фазы можно пренебречь. Такое
допущение вполне правомочно, особенно когда объемная
доля дискретной фазы не очень велика.
В этом случае задача замыкания системы сводится к
определению сил межфазного взаимодействия и теплооб-
мена между фазами.
176
7.1. Подход «Эйлер—Эйлер»
Сила межфазного взаимодействия складывается из не*
скольких составляющих, главной из которых является сила
лобового сопротивления FaD. Для этой силы используется
следующее выражение:
vD_ 3 Cp.ap(u-ua)|u-ug|pg _
Q
8 Г°	(7.4)
_ 3 CDi„p(u-u,)|u-n,|>r,
8
где p,u — соответственно плотность и вектор скорости
непрерывной среды; га — радиус дискретной частицы;
CD g — коэффициент сопротивления.
В зависимости от конкретной задачи для коэффициента
сопротивления используются различные выражения. Рас-
смотрим некоторые из них.
1)	Закон Стокса:
24
Cd“ = rT’ Re“^1’	(7.5)
где
(7.6)
р
где ц — коэффициент динамической вязкости непрерыв-
ной среды.
177
Глава 7. Особенности расчёта многофазных течений
Закон Стокса получен аналитически для очень малых
значений числа Рейнольдса и для большинства практиче-
ских задач требует корректировки.
2)	Модель Шиллера Наумана:
со,. =^-(l+0.15Re.ae’).
Re„	7
(7-7)
Модель Шиллера Наумана является наиболее универ-
сальной и чаще всего применяется для расчета дискретных
частиц, имеющих небольшую объемную долю.
3)	Модель Вэня Ю (Wen Yu Drag Model):
24 /	----0.687
max	1 + 0.15Re„
LRea '
,0.44
(7-8)
Rea = kc Ree
где кс — объемная доля непрерывной среды.
Модель применима для твердых частиц, занимающих в
течении достаточно большую объемную долю.
Существуют и другие модели для расчета коэффициен-
та сопротивления.
Для дискретной фазы источник в уравнении энергии
складывается из трех составляющих:
178
7.1. Подход «Эйлер—Эйлер»
Я. phase,a	<lconv,a	Я rad ,а9
(7.9)
где qphaseia — тепло фазового перехода; qcom а — тепловой
поток, сбрасываемый частицей за счет конвекции; qrad a —
лучистый тепловой поток, сбрасываемый частицей.
Рассмотрим более подробно конвективный тепловой
поток. Он определяется через критерий Нуссельта:
9^=Р.С,(Т.-Т),	(7.10)
Га Ps,a
где Я — коэффициент теплопроводности непрерывной
среды.
Для определения критерия Нуссельта используются
различные модели.
1) Модель Ранца—Маршалла:
Nua = 2 + 0.6 Re/2 * * * Рг,/3,	(7.11)
где Рг — число Прандтля непрерывной среды.
Эта модель является наиболее универсальной.
2) Модель Томиямы:
Nua = 2 + 0.15Re/8 Pr05.	(7.12)
Модель применяется к турбулентным пузырьковым те-
чениям с относительно низким числом Рейнольдса.
179
Глава 7. Особенности расчёта многофазных течений
3) Модель Кавано:
Nue0
NU“ “ l + (3.42Nue>0M/(RecPr)) ’	13)
где Nua 0 = 2 + 0.459 Rea°55 Рг033 , а число Маха равно:
1и—и I
M=J------£1,	(7.14)
а
где а — скорость звука в непрерывной среде.
Основное преимущество этой модели заключается в
том, что она применима и для очень разреженного газа.
Методы VOF (Volume of Fluid — объем жидкости)
Описанный выше метод «Эйлер—Эйлер» можно суще-
ственно упростить для случая, когда рассматриваются те-
чения с четко выраженной границей, когда течение разде-
лено на достаточно большие области, в которых объемная
доля одной из фаз равна единице, а всех других — нулю.
Примеры: течение реки в русле, каскады или небольшие
водопады, заправка топливного бака, движение объемных
пузырей в жидкости. При этом области с мелкодисперс-
ными каплями (пузырьками) должны отсутствовать.
Методы VOF используют уравнения вида, аналогично-
го уравнениям метода «Эйлер—Эйлер». Главная отличи-
тельная особенность методов VOF — это то, что в ходе
моделирования решается только один набор уравнений
180
7.1. Подход «Эйлер—Эйлер»
движения, энергии и турбулентности, а также специальное
уравнение для отслеживания положения межфазной гра-
ницы (аналог уравнения неразрывности). Значения пара-
метров, характеризующих свойства фазы (плотность, вяз-
кость, теплопроводность и т.п.) при переходе через меж-
фазную поверхность «переключаются», чтобы в каждой
точке пространства соответствовать той фазе, которая в
этой точке находится.
Методы VOF отличаются друг от друга в основном
способами моделирования межфазной границы.
Главным преимуществом описанных методов является
их простота по сравнению с методом «Эйлер—Эйлер»
(ведь в каждой точке пространства решается только один
набор уравнений), а главным недостатком — то обстоя-
тельство, что размер расчетной сетки должен соответство-
вать минимальному объему, который в рассматриваемом
течении может занимать каждая из фаз. Например, если в
течении присутствуют капли жидкости, то размер расчет-
ной сетки должен определяться размерами наименьшей из
возможных капель (при этом необходимо учитывать, что
капли при течении могут дробиться на более мелкие).
7.2. Подход «Эйлер—Лагранж»
Как уже говорилось, при использовании метода Ла-
гранжа для непрерывной фазы решаются уравнения На-
вье—Стокса, а решение для дисперсной фазы находится
путем отслеживания большого количества частиц, пузырь-
ков или капелек в предварительно вычисленном поле тече-
ния. При этом уравнения, решаемые для каждой дисперс-
ной частицы, записываются в т.н. постановке Лагранжа,
т.е. в системе координат самой частицы.
Дисперсная фаза может обмениваться импульсом,
массой и энергией с жидкой фазой. Траектории частиц
или капель вычисляются индивидуально через опреде-
ленные промежутки времени в процессе вычислении не-
прерывной фазы.
Рассмотрим случай, когда объемная доля частиц неве-
лика, и взаимодействием частиц можно пренебречь. При
этом частицы могут иметь размеры, находящиеся в неко-
тором диапазоне: ra)inin < га < га тт.
Полидисперсная смесь частиц представляется в виде
набора L групп частиц, каждая из которых характеризуется
182
7.2. Подход «Эйлер—Лагранж»
значениями радиуса га, плотности ра = па • та, компонент
скорости иа j и температуры Та. Здесь па — концентрация
частиц (число частиц в единице объема); та — масса од-
ной частицы. Предполагается, что частицы имеют сфери-
ческую форму, химически инертны по отношению к не-
прерывной фазе и не взаимодействуют между собой.
Для каждой группы частиц (а= 1,2, ...,£) уравнения,
описывающие движение частиц, включают в себя:
1) уравнение сохранения количества движения:
du .
Pa~f1L~^aj»	(7.15)
dr
2) уравнение сохранения энергии:
dT
PaCs-^ = Wa.	(7.16)
dr
Для подхода Лагранжа удобнее использовать не
векторную форму записи уравнений, а индексную.
Здесь —----полная производная по времени (берется
dr
вдоль траектории частицы); Fa i — сила воздействия газа
или жидкости на частицу (см. формулу (7.4); Wa — тепло,
поступающее в частицу (может включать в себя а —
тепловой поток, сбрасываемый частицей за счет
183
Глава 7. Особенности расчёта многофазных течений
конвекции, аги...г — тепло фазового перехода, а —
лучистый тепловой поток, сбрасываемый частицей); Cs —
теплоемкость частиц, которая предполагается постоянной.
Необходимо иметь в виду, что при рассмотрении взаи-
модействия частиц и газа/жидкости могут использоваться
две модели взаимодействия:
1) упрощенная (односторонняя), в которой учитывается
только воздействие газа или жидкости на частицы, а воз-
действием частиц на непрерывную среду пренебрегают;
2) полная (двухсторонняя), в которой учитывается и
воздействие частиц на газ или жидкость.
Во втором случае в уравнения движения непрерывной
среды добавляются дополнительные члены:
F> = -ffaj’ W =	-Fa^).	(7.17)
cr==l	a-1
7.3. Массообмен между частицами
и непрерывной средой
Изменение массы компонента А, которое происходит в
результате массообмена между непрерывной фазой и час-
тицей, удовлетворяет уравнению:
^ = -2^.pDSh(EC^-C,),	(7.18)
где тл — масса компонента А в частице; pD — коэффи-
циент динамической диффузии в непрерывной среде; Е —
равновесное отношение массовых долей (как правило,
принимается равным единице); Сл — массовая доля ком-
понента А в частице; Сл а — массовая доля компонента А в
непрерывной среде; Sh — критерий Шервуда, для которо-
го используется следующая формула:
z Х1/3
Sh = 2 + O.6Re05 — .	(7.19)
В случае, когда происходит испарение частицы, лучше
использовать модель испарения жидкости (Liquid Evaporation
Model), которая использует два соотношения для массопе-
185
Глава 7. Особенности расчёта многофазных течений
реноса в зависимости от того, находится ли капля выше
или ниже точки кипения.
Точка кипения определяется через уравнение Антуана
(соотношение между давлением паров и температурой для
чистых компонентов):
Pvap
= Рscaled А~
(7.20)
где Тр — температура частицы в градусах Цельсия; А, В, С—
эмпирические константы для каждого чистого компонента;
Аса/. — масштаб давления (см. ниже). В Табл. 7.1 приведе-
ны значения этих констант для некоторых жидкостей:
Табл. 7.1. Константы в уравнении Антуана
для некоторых жидкостей
Вещество	А	В	С	тш1п,°с	Т^’С
Вода	8.07131	1730.63	233.426	1	100
Вода	8.14019	1810.94	244.485	99	374
Этанол	8.20417	1642.89	230.300	-57	80
Этанол	7.68117	1332.04	199.200	77	243
Метан	6.34159	342.22	260.221	-181	-163
Метан	6.7021	394.48	264.609	-162	-83
Надо иметь в виду, что в формуле (7.20) используются
градусы Цельсия, и получена она для давления в милли-
метрах ртутного столба. Поэтому для того, чтобы получить
давление в Па (система СИ), необходимо использовать
следующий масштаб давления:
186
7.3. Массообмен между частицами и непрерывной средой
Рscale
101325 Па
760 ммрт.ст.
= 133.32.
(7.21)
Если давление выше давления окружающего газа,
частица закипает.
Расчет массообмена при испарении осуществляется
следующим образом.
Если температура частицы находится выше точки ки-
пения, массоперенос определяется по формуле:
^тр _ Qc+Qr	/7
dr 91	’	(	'
где 91 — скрытая теплота парообразования; QC,QR — со-
ответственно конвективный и радиационный тепловые по-
токи к частице.
Если температура частицы находится ниже точки кипе-
ния, массоперенос определяется по формуле:
dmp *	_ с ,
—- = -2nrppD Sh —- In
dr	MG
s__
'V
yap 7
(7.23)
где mp —масса частицы; rp —радиус частицы; Mc, MG —
молекулярные массы пара и газовой смеси в непрерывной
фазе; X* —равновесная мольная доля пара испаряющегося
компонента на поверхности капли; Х^ — мольная доля
испаряющегося компонента в газовой фазе.
7 А. Численный метод решения
уравнений переноса дискретной фазы
Два принципиальных подхода к решению уравнений
переноса дискретной фазы: метод Лагранжа (отслежива-
ются траектории каждой частицы) и метод Эйлера (реша-
ются уравнения подобные уравнениям для непрерывной
среды — газа), отличаются и численным методом реше-
ния. Во втором случае расчетные сетки частиц и газа сов-
падают, в методе Лагранжа траектории частиц естественно
отличаются от расчетной сетки газа.
Оба метода обладают как преимуществами, так и недос-
татками по сравнению друг с другом. Для ряда задач наи-
более важно то, что метод Лагранжа позволяет максималь-
но точно отследить эволюцию каждой частицы, в частно-
сти, фазовые переходы. При использовании метода Эйлера
с отслеживанием изменений частицы возникают проблемы —
частицы, в которых уже началась кристаллизация, на се-
точной модели непрерывной среды могут оказаться в об-
ласти, где кристаллизация еще невозможна с точки зрения
предыстории траекторий частиц; возможна и обратная си-
188
7.4. Решение уравнений переноса дискретной фазы
туация. Кроме того, метод Эйлера требует существенно
больше компьютерных ресурсов, так как для каждой груп-
пы частиц приходится решать фактически полную систему
уравнений типа Навье—Стокса. В методе Лагранжа реша-
ются обыкновенные дифференциальные уравнения — по-
этому она значительно экономичнее.
С другой стороны, при использовании метода Лагранжа
возникают проблемы с решением уравнения неразрывно-
сти (7.1) на сеточной модели частиц. Дело в том, что тра-
ектории частиц могут иметь весьма сложную геометрию и
неоднократно пересекаться между собой. По этой же при-
чине проблематичен пересчет сил воздействия частиц на
газ и соответствующих конвективных потоков с сеточной
модели частиц на сеточную модель газа/жидкости. Даже в
двумерном случае эта проблема решается довольно слож-
но. Для трехмерного течения задача еще более усложняет-
ся: частицы сбиваются в узкие области, и на очень близком
расстоянии могут оказаться частицы с совершенно разны-
ми скоростями и температурой, т.е. нарушается непре-
рывность решения. Интерполяция параметров для пересче-
та параметров с траекторий частиц на сеточную модель
газа и наоборот становится практически невозможной.
Общий вид уравнений (7.15) и (7.16) выглядит следую-
щим образом:
189
Глава 7. Особенности расчёта многофазных течений
% = C(g-ga) + a.
(7.24)
Уравнение (7.24) решается на каждом шаге по времени
bi = bxlu на участке от т = 0 до г = Аг. Считается, что
на этом участке С, g, а постоянны и определяются по па-
раметрам на n-ом слое. Аналитическое решение уравнения
(7.24) имеет вид:
gan^=g«(^) =
= g + (gan ~g)e-CAr+а[1 -e’CAr] /С.
(7-25)
Используя решение (7.24), находим конечные точки
траекторий:
V+1 =хап+ил.т + (иап - м)[1 - ехр(-QAг)] / Cf,
Уап+1=Уа” +vAr + (v/ -v)^l-exp(-C/Ar)J/C/,	(7.26)
zan+i = zan + wAt + (w/ - W)[l - exp(-CyAr^ / Cf.
7.5. Общие рекомендации
по численному моделированию
многофазных течений
В данном разделе были рассмотрены основные подходы
к моделированию гетерогенных (многофазных) течений.
Можно было заметить, что основной особенностью мо-
делирования таких течений вне зависимости от выбранно-
го подхода («Эйлер—Эйлер» или «Эйлер—Лагранж») яв-
ляется использование эмпирических (полуэмпирических)
зависимостей для учета межфазного взаимодействия.
Таким образом, при моделировании многофазных по-
токов необходимо очень внимательно следить за тем, что
применяемые модели межфазного взаимодействия исполь-
зуются в диапазонах значений основных параметров, для
которых эти модели были созданы. Данный перечень па-
раметров включает различные критерии подобия (числа
Рейнольдса, Прандтля, Вебера и т.д.), объемную долю фаз
в потоке и ряд других. В случае, когда какой-либо из пара-
метров выходит за пределы «рабочего» диапазона модели,
это может привести к физически неадекватным результа-
191
Глава 7. Особенности расчёта многофазных течений
там или существенно повлиять на сходимость численного
решения.
Исходя из этого, перед началом решения задачи о мно-
гофазном течении необходимо заранее представлять себе
его природу и приблизительные характеристики (в особен-
ности — относительно объемных долей каждой из фаз в
потоке и скоростные характеристики течения). Основываясь
на этих данных выбирают соответствующий метод моде-
лирования и модели межфазного взаимодействия. В про-
цессе решения значения определяющих критериев меж-
фазного взаимодействия тщательно отслеживают. В случае
их выхода за пределы «рабочего» диапазона какой-либо
модели межфазного обмена, рассматривается вопрос о
применении другой модели.
Заключение
При изучении данного пособия может возникнуть впе-
чатление о некоторой его незавершенности и недосказан-
ности. Однако, следует обратить внимание читателя, что
это всего лишь обзорный курс.
Проблемам вычислительной гидродинамики и тепло-
обмена посвящены обширные труды и монографии. Здесь
же была сделана попытка обзора самых важных идей и из-
ложения их как можно более простым языком.
Существует множество CFD пакетов, но правильно и в
полной мере ими может пользоваться лишь тот, кто обла-
дает хотя бы начальными знаниями в области теории.
В данном пособии как раз была сделана попытка дать этот
минимум теоретических знаний. Именно поэтому мы не
рассматривали вопросы лучистого теплообмена, примене-
ние метода конечных элементов к задачам механики жид-
кости и газа и т.д.
Список
использованной
литературы
1.	ДульневГ. Н. и др. Применение ЭВМ для решения задач теплооб-
мена, М.: Высшая школа, 1990. 207 с.
2.	Власова Е. А., Зарубин В, С., КувыркинГ. Н. Приближенные мето-
ды математической физики, М.: Изд-во МГТУ им. Н. Э. Баумана,
2001.700 с.
3.	Коренев Г. В. Тензорное исчисление: Учеб, пособие: Для вузов.
М.: Изд-во МФТИ, 2000. 240 с.
4.	FerzigerJ. Н., Peric М. Computational Methods for Fluid Dynamics 3, rev.
ed. Berlin; Heidelberg; New York; Barcelona; Hong Kong; London;
Milan; Paris; Tokyo: Springer, 2002.
5.	Патанкар С. Численные методы решения задач теплообмена и
динамики жидкости / Пер. с англ. М.: Энергоатомиздат, 1984.
152 с.
6.	Мышкис А. Д. Прикладная математика для инженеров. Специаль-
ные курсы. М.: Физматлит, 2007. 688 с.
7.	МышкисА. Д. Лекции по высшей математике. М.: Наука, 1973.
640 с.
8.	Андерсон Д., Таннехилл Дж., Плетчер Р. Вычислительная гидро-
механика и теплообмен: В 2-х т. Т. 1: Пер. с англ. М.: Мир, 1990.
384 с.
9.	Самарский А. А., Вабищевич П. Н. Вычислительная теплопередача.
w М.: Книжный дом «Либроком»/ЦК88, 2014. 784 с.
10.	Белов И. А., ИсаевС. А. Моделирование турбулентных течений:
Учебное пособие, Балт. гос. техн. ун-т. СПб., 2001. 108 с.
194
Список использованной литературы
И. Волков К. Н.9 Емельянов В. Н. Моделирование крупных вихрей в
расчетах турбулентных течений. М.: Физматлит, 2008. 368 с.
12.	Лойцянский Л. Г. Механика жидкости и газа. Изд. 5-е, перерабо-
танное, Главная редакция физико-математической литературы из-
дательства «Наука», М., 1978,736 с.
13.	Монин А. С.9 ЯгломА.М. Статистическая гидромеханика, Часть!.
Главная редакция физико-математической литературы издательст-
ва «Наука», М., 1965,640 с.
14.	Альбом течений жидкости и газа: Пер. с англ. / Сост. М. Ван-Дайк.
М.:Мир. 1986. 184 с.
15.	ДрейцерГ.А. Основы конвективного теплообмена в каналах:
Учебное пособие. М.: МАИ, 1989. 84 с.
16.	CiofaloM. Large-Eddy Simulation: A Critical Survey of Models and
Applications // Advances in Heat Transfer, Vol. 25. P. 321—419.
17.	Сергель О, С. Прикладная гидрогазодинамика: Учебник для авиаци-
онных вузов. М.: Машиностроение, 1981. 374 с.
18.	Ландау Л. Д., Лифшиц Е.М. Теоретическая физика: Учеб. Пособ.:
Для вузов. В 10 т. Т. VI. Гидродинамика. 5-е изд., стереот. М.: Физ-
матлит, 2006.736 с.
19.	ХайрерЭ., Ваннер Г. Решение обыкновенных дифференциальных
уравнений. Жесткие и дифференциально-алгебраические задачи.
М.: Мир, 1999.
20.	Молчанов А. М. Расчет струй с неравновесными химическими ре-
акциями И Современные проблемы теплообмена в авиационной
технике. М., 1983. С. 15-19.
МАТЕМАТИКО-ХУДОЖЕСТВЕННЫЙ МИРОВОЙ БЕСТСЕЛЛЕР
ядя Петрос
ПРОБЛЕМА
ГОЛЬДБАХА
И
«Архимед» будут помнить, кода Эсхила уме забудут, потому что языки умирают, а матема-
тические идеи - нет. Может быть, «бессмертие» - глупое слово, ио что бы оно ни значило,
вероятно,лучшие шансы его достигнуть- у математика»
Годфри Харди. Апология математика
Перед читателем - удивительная книга, переведенная на все основные языки
мира и ставшая мировым интеллектуальным бестселлером. Ома повествует об ис-
тории едкого математика, Петроса Папахрмстоса, оставшегося до конца верным
великой идее познания истины, всю свою жизнь посвятившего поиску решения
одной-единствениой научной проблемы, оказавшейся не по силам великим умам
прошлого.
Это не просто рассказ о юноше и его дяде, познакомившем племянника с простыми числами,
красота которых была столь велика,что предопределила путь всей его жизни; это трагедия мате-
матика, напомнившая об играх разума Джона Нэша, Алана Тьюринга, Курта Геделя и других
великих подвижников науки.
Автор книги Апостолос Дрксиадис - греческий писатель, яркий представитель современной
интеллектуальной прозы. Благодаря его таланту повествованиеостранномдядюшке,искавшем
решение одной из труднейших задач во всей математике, становится увлекательным пригла-
- шением в мир научных открытий и математических идей.
Книга адресована самому широкому кругу читателей - от математиков и философов, историков
и методологов науки до тех, кто только подходит к выбору дела своей жизни; она напомнит
им, что математика признаеттолько величайших, но хорошо помнит каждого, кто был ей предан.
ЗАГАДКИ СФИНКСА
и другие математические головоломки
Перед читателем - впервые переведенная
на русский язык работа выдающегося
американского математика
и популяризатора науки Мартина Гарднера.
Автор как всегда остается верен своему
уникальному стилю, который характеризуют
яркость, доходчивость, тонкий юмор, блеск
мысли, постоянное вовлечения читателя
в самостоятельное творчество.
В книге представлены занимательные математические задачи и головоломки,
публиковавшиеся автором в течение ряда лет на страницах «Журнала
научной фантастики Айзека Азимова». Эти задачи, составленные самим
М. Гарднером, увлеченными читателями его журнальной колонки, друзьями
и коллегами автора, в равной мере относятся как к миру собственно
математики, так и к миру логических парадоксов, многие из которых
изложены в виде фантастических историй на загадочных далеких планетах.
Материал книги построен таким образом, что к каждой из головоломок
дается ответ; который в большинстве случаев порождает новые
вопросы, связанные с соответствующей темой. Решение этих новых
вопросов приводится уже на следующем уровне ответов.
Всего в книге четыре раздела с ответами и решениями, так что читатель
может последовательно переходить с одного раздела на другой, постепенно
углубляя свое понимание.

БЕРИЯ И СОВЕТСКИЕ УЧЕНЫЕ
В АТОМНОМ ПРОЕКТЕ
Книга 1
Книга 2
Выдающиеся
УЧЕНЫЕ-ЯДЕРЩИКИ
Советского Союза
Судьба
Лаврентия
Берии
Судьба
Лаврентия
БЕРИИ |
Создание атомного оружия в Советском Союзе - одна из выдающихся страниц истории нашей страны.
В реализации этого невероятно трудного технологического прорыва принимало участие несколько сотен
тысяч человек, но особая ответственность легла на интеллектуальное «ядро» создателей атомного оружия.
В первой части монографии НА Кудряшова обсуждается участие в Советском атомном проекте
выдающихся ученых и организаторов производства: Б. Л. Ванникова, П. Л. Капицы, И. В. Курчатова,
Ю. Б. Харитона, К. И. Щёлкина, Я. Б. Зельдовича, И. Е. Тамма, А. Д Сахарова, А. П. Александрова, А. Л. Минца
и М. Г. Мещерякова. Связующим звеном сюжета книги являются отношения ученых и председателя
специального комитета при Совете Министров СССР Лаврентия Берии.
Книга содержит много любопытных эпизодов и фактов из истории создания атомного оружия
в Советском Союзе.
Вторая часть научно-исторического исследования посвящена описанию жизни
и деятельности руководителя специального комитета при Совете Министров СССР Лаврентия Берии.
Дается анализ его роли в Советском атомном проекте и его участия в разрешении идеологического спора
между физиками и философами, позволившего сохранить высокий уровень физической науки в СССР.
Обсуждается попытка проведения реформ в Советском Союзе, предпринятая Берией после смерти
Сталина. Описаны события, происходившие после ареста Лаврентия Берии: выступления на июльском
пленуме его «вчерашних друзей».
Рассматриваются материалы следствия по делу Берии и его расстрел. Даны детальные характеристики
его помощников и заместителей. Книга содержит много новых фактов о деятельности Берии
и его окружения.
ЭЛЕМЕНТЫ ПРИКЛАДНОЙ
МАТЕМАТИКИ--------
элементы й Ь. Зельдович, А. Д Ммиюк
ПРИКЛАДНОЙ математики
Я?	Книга,
которую следовало бь
казватд иначе...
Книга, которую
следовало бы

назвать иначе...
В задачах физики, техники
и практических вычислений
используются численные
и графические методы, ряды.
В книге содержатся полезные
приемы таких вычислений.
В наглядной форме даются основные сведения о комплексном
переменном, линейных дифференциальных уравнениях,
векторах и векторных полях и вариационном исчислении.
Формальные доказательства в большинстве случаев заменены
наводящими соображениями; за счет этого достигнуто
упрощение и облегчено применение математических понятий.
Подробно анализируются некоторые физические задачи,
в частности относящиеся к оптике, механике, теории
вероятностей.
iritfiirtWAiaarrt	и» if tortmir	Tnrr im  i	ом rui^ AairfthMMrui.irifei
Изобретая ИНСТРУМЕНТЫ
ИЗОБРЕТАЯ
ИНСТРУМЕНТЫ
НАУКИ БУДУЩЕГО
лТРИЗ
«ЖЛЦиЫЙ
ЕЛ&НА СЕРАЯ
НАУКИ БУДУЩЕГО
УСКОРЯЮЩАЯ
НАУКУ ТРИЗ
Физика ускорителей, лазеров, плазмы
Читателя ждет увлекательный, не имеммций себе аналогов
по содержательности обзор принципов работы ускорительной техники,
основанных на синтезе физики платы, лазеров и ускорителей.
Развитие инструментов ускорительной науки является ярчайшим примером научного
прогресса, а тот факт, что более трети Нобелевских лауреатов по физике напрямую
связаны с использованием ускорителей, весомое тому подтверждение.
Век назад Резерфорду для изучения структуры атома д остаточно было лишь образца
урана, а в качестве помощников хватило нескольких аспирантов; для задач
же современной физики требуются большие коллективы исследователей и гигантские
ускорители, что неуклонно требует изменений в самом подходе к процессу инноваций.
Отвечал вызовам науки будущего, авторы книги адаптируют для ученых
методику ТРИЗ, выделяя изобретательские принципы Ускоряющей науку
ТРИЗ (УН-ТРИЗ) на примере прогресса ускорительной науки. В книге
предлагаются задачи для размышления и примеры применения этого
нового подхода.
Книга будет интересна как физикам, работающим в областях физики плазмы, лазеров
и ускорителей,так и в целом ученым и изобретателям, желающим усовершенствовать
свой подход к исследовательскому процессу. Много полезного для себя найдут в книге
и студенты, а также широкий круг заинтересованных читателей.
Леонид Владимирович БЫКОВ
Кандидат технических наук, доцент, директор Института дополнительного
профессионального образования Московского авиационного института
(МАИ), доцент кафедры «Авиационно-космическая теплотехника», руково-
дитель группы моделирования газовой динамики и тепломассообмена.
Автор более 60 научных статей и учебных пособий. Основные научные
интересы—в области моделирования высокоскоростных течений.
Александр Михайлович МОЛЧАНОВ
Доктор технических наук. Доцент кафедры «Авиационно-космическая
теплотехника» Московского авиационного института (МАИ). Автор более
100 научных статей, учебных пособий и монографий в области вычисли-
тельной механики. Основные научные интересы: химически и термически
неравновесные высокоэнтальпийные течения, турбулентность, неравно-
весное излучение.
Дмитрий Сергеевич ЯНЫШЕВ
Кандидат технических наук, доцент кафедры «Авиационно-космическая
теплотехника» Московского авиационного института (МАИ). Автор более
30 работ и трех учебных пособий по прикладной вычислительной механи-
ке жидкости и газа. Основные научные интересы сосредоточены в области
турбулентных течений и нестационарных процессов.
Наше издательство предлагает следующие книги:
СПЛОШНЫХ СРЕД
в задачах
авидцмонной
И КОСМИЧЕСКОЙ
ТЕХНИКИ
8
АНАЛИТИЧЕСКИЕ МЕТОДЫ
ПОРКИ
ППЛОПР ОВ ОДИОСТИ
нее ।
ИП0ЖЕ1«
В. В. Нмыцкий
В. В. Степанов
U1ECIIEIIAI
EEZPII
KietEkJMMbl'II
УРЙВИЕВИЙ
URSSi
Каталог изданий
в Интернете:
http://URSS.ru
E-mail: URSS@URSS.ni
117335, Москва, Телефон/факс
Нахимовский	(многоканальный)
проспект, 56	+7 (499) 724 25 45
Отзывы о настоящем издании, а также обнаруженные
опечатки присылайте по адресу URSS@URSS.nj.
Ваши замечания и предложения будут учтены
и отражены на web-странице этой книги на сайте
httpVAJRSS-nj