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

КОНФЕРЕНЦИИ, КНИГИ, ПОСОБИЯ, НАУЧНЫЕ ИЗДАНИЯ

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


Pages:     | 1 || 3 | 4 |   ...   | 6 |

«Т. Г. Елизарова КВАЗИГАЗОДИНАМИЧЕСКИЕ УРАВНЕНИЯ И МЕТОДЫ РАСЧЕТА ВЯЗКИХ ТЕЧЕНИЙ Москва Научный Мир 2007 ...»

-- [ Страница 2 ] --

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

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

2.8 Разностная аппроксимация уравнения Больцмана и кинетически-согласованные разностные схемы Численное решение уравнения Больцмана (2.1) представляет собой существенные сложности, связанные с большой размерностью зада 54 Гл.2. Элементы кинетической теории газов чи и проблемами аппроксимации и вычисления интеграла столкно вений. Развитые для этой задачи численные алгоритмы приведены, в частности, в [9, 124].

Левая часть уравнения (2.1) имеет вид уравнения переноса функции f со скоростью, и может быть аппроксимирована с помо щью разностной схемы первого порядка точности с направленными разностями схемы Куранта–Исааксона–Рисса. При согласованном выборе шагов по времени и по пространству h = ||t эта разност ная схема описывает перенос заданного возмущения функции f без искажений [87, 88].

Приведем вид соответствующей разностной схемы для одномер ного пространственного течения в отсутствие внешних сил на рав номерной пространственной сетке с шагом h = xi+1 xi и шагом по времени t = tj+1 tj fij+1 fij fi fi (2.34) + x = I(f, f )h,i, k 0, t h fij+1 fij fi+1 fi + x = I(f, f )h,i, k 0.

t h компонента скорости частицы, f = f j = f (xi,, tj ), Здесь k I(f, f )h,i разностный аналог интеграла столкновений.

Разностная схема (2.34) тождественным образом преобразуется к эквивалентному виду fij+1 fij fi+1 fi + k t 2h h fi+1 2fi + fi = I(f, f )h,i. (2.35) |k | h Пусть для разностного аналога интеграла столкновений I(f, f )h,i выполняются условия ортогональности сумматорным инвариантам h(). Разностную схему (2.35) можно осреднить в пространстве ско ростей, предполагая, что f = f (0). При этом все интегралы удается вычислить аналитически и сразу построить разностную схему для 2.8. Кинетически-согласованные разностные схемы газодинамических величин – плотности, скорости u и энергии E [30, 51, 105, 106, 114].

Полученная схема оказывается довольно громоздкой, поскольку она включает в себя интегралы ошибок, возникающие при осредне нии модулей скоростей k. Тем не менее эта схема оказалась весьма эффективной для решения уравнений Эйлера и была использована при проведении расчетов как одномерных, так и некоторых двумер ных сверхзвуковых газодинамических течений. Позднее близкие по структуре разностные уравнения были получены в работе [144].

Заменим в разностной схеме (2.35) коэффициент перед послед ним слагаемым, полагая || c, где c – скорость звука. Вводя в рассмотрение время = h/2c, характеризующее скорость пересече ния частицей разностной ячейки, получим hk h h2 (2.36) |k | = = k = k.

2 2|k | 2c Схема (2.35) примет вид fij+1 fij fi+1 fi1 2 fi+1 2fi + fi = I(f, f )h,i. (2.37) + k k h t 2h Выражение (2.37 ) в случае f = f (0) также допускает аналитиче ское осреднение с сумматорными инвариантами h() и приводит к построению более изящной по сравнению с предыдущим случа ем системы разностных уравнений для описания макроскопических величин в газе. Полученные таким образом схемы были названы кинетически-согласованными разностными схемами КСРС.

С использованием принятых в [92] обозначений для центрально разностных производных первого и второго порядка по простран ству, перепишем разностную схему (2.37) в виде fij+1 fij (2.38) + (f k ) (f k k )xx = I(f, f )hi.

x t Здесь f = (fi+1 fi1 )/h – центральная разностная производная x первого порядка, fx = (fi+1 fi )/h и fx = (fi fi1 )/h – левая и правая разностные производные, fxx = (fx fx )/h – центральная разностная производная второго порядка.

56 Гл.2. Элементы кинетической теории газов Применяя введенные обозначения, запишем КСРС для случая плоского одномерного течения в виде j+1 j + (u) = (u2 + p)xx, (2.39) i i x t j+1 j (u)i (u)i + (u2 + p) = (u3 + 3pu)xx, x t j+1 j Ei Ei p + (u(E + p)) = (u2 (E + 2p))xx + ( (E + p))xx.

x t Модифицированные позднее КСРС оказались очень эффективными при численном моделировании широкого круга газодинамических течений [6–8, 24, 25, 30, 31, 50–52, 105, 106].

Дифференциальные аналоги этих схем легли в основу первых вариантов КГД уравнений.

Глава Квазигазодинамические уравнения В этой главе приведено два варианта построения квазигазодинами ческой системы уравнений, которые позволяют получить конкрет ный вид выписанных ранее без вывода выражений для векторов плотности потока массы jm, тензора вязких напряжений и векто ра теплового потока q.

Первый вариант основан на использовании достаточно простой кинетической модели движения частиц. Эта модель представляет собой циклически-повторяющийся процесс бесстолкновительного разлета и последующих столкновений частиц с установлением максвелловского равновесия. Изложенный вывод КГД уравнений устанавливает связь параметра сглаживания с максвелловским временем релаксации. Построенные таким образом КГД уравнения представляются затем в виде законов сохранения. При этом выра жения для jm, и q записываются как соответствующие выражения для уравнений Навье–Стокса с добавками, пропорциональными.

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

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

Изложение этих результатов приведено в соответствии с работа ми [43, 49–53, 114, 115, 146, 158].

58 Гл.3. Квазигазодинамические уравнения 3.1 Регуляризованное кинетическое уравне ние Решение интегро-дифференциального уравнения Больцмана явля ется сложной задачей, поэтому к настоящему времени развиты упро щенные кинетических модели, которые позволяют находить прибли женные решения отдельных задач. Многие модели основаны на идее расщепления задачи по физическим процессам, которая иначе на зывается принцип суммарной аппроксимации. Именно такой подход лежит в основе DSMC метода, описанного в п. 2.7.2.

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

Эта модель представляет движение газа как цикличе ски повторяющийся процесс, состоящий из двух этапов: это бесстолкновительный разлет молекул газа и последующее мгновенное установление тер модинамического равновесия за счет столкновения частиц этап мгновенной максвелли зации. Схематическое изобра Рис. 3.1. Схема, поясняющая кинети- жение этой модели приведено на рис. 3.1. Детальный анализ ческую модель этой модели и вариантов полу ченных на ее основе модельных кинетических уравнений приведен в [51, 52, 105, 106] и здесь не обсуждается.

Пусть в некоторый момент времени t = t1 функция распределе ния имеет локально-максвелловский вид:

( u) f (0) (x,, t) = (3.1) exp.

(2RT )3/2 2RT Затем в течение времени t происходит бесстолкновительный раз 3.1. Регуляризованное кинетическое уравнение лет молекул, который описывается уравнением Больцмана для сво бодномолекулярного течения f + ( · )f = 0.

t Это уравнение представляет собой линейное уравнение переноса и имеет точное решение (3.2) f (x,, t) = f0 (x t, ), где f0 = f0 (x, ) известная функция распределения в момент вре мени t = 0.

Далее, в момент времени t = t2 = t1 + t функция распределе ния f мгновенно вновь становится локально-максвелловской (3.1), но уже с новыми значениями макроскопических параметров, u и T.

Мгновенная максвеллизация имитирует этап столкновения молекул, который в уравнении Больцмана описывается интегралом столкно вений I(f, f ). Затем оба этапа циклически повторяются.

Считая время бесстолкновительного разлета t достаточно ма лым, разложим функцию распределения (3.2) в ряд Тейлора, пола гая функцию f0 максвелловской f0 = f (0). Ограничиваясь членами второго порядка малости по t, получим:

t f (x,, t) = f (0) t( · )f (0) + ( · )( · )f (0). (3.3) Параметр разложения t · при больших значениях скоростей нельзя считать малым, и отбрасываемые при таком разложении члены могут оказаться существенными. Однако в данном случае их вкладом при || 2RT можно пренебречь, так как сама функция распределения f (0) и все ее производные экспоненциально убывают с ростом.

В равенстве (3.3) перенесем все слагаемые в левую часть и раз делим обе части уравнения на t, заменяя разностную производную по времени в первом слагаемом на дифференциальную. В результа те получим уравнение для функции распределения в виде f (0) f f t + ( · )f (0) ( · ) ( · )f (0) =.

t 2 60 Гл.3. Квазигазодинамические уравнения Здесь в правую часть полученного уравнения добавлен интеграл столкновений в БГК приближении, обеспечивающий релаксацию функции распределения на новом временном слое к максвелловско му распределению.

Отождествляя временной интервал t/2 с максвелловским вре менем релаксации, получим окончательный вид регуляризованно го кинетического уравнения f (0) f f + ( · )f (0) ( · ) ( · )f (0) = (3.4).

t Заметим, что это уравнение представляет собой дифференциаль ный аналог разностной схемы Куранта–Иссаксона–Рисса для БГК уравнения, выписанной в виде (2.38) для случая = const.

При осреднении уравнения (3.4) с сумматорными инвариантами, как будет показано далее, получается система КГД уравнений для = 5/3, P r = 1, Sc = 1.

Уравнение (3.4) формально можно выписать исходя из уравне ния Больцмана в БГК приближении:

f (0) f f (3.5) + ( · )f =, t заменяя функцию распределения в конвективном слагаемом урав нения (3.5) на ее приближенное значение вида f = f (0) ( · )f (0). (3.6) Заметим, что формальная замена функции распределения f в конвективном слагаемом (3.5) на функцию распределения Навье– Стокса (2.11) позволяет, после осреднения с сумматорными инвари антами, получить систему уравнений Навье–Стокса.

В [158] было показано, что соотношения (2.11) и (3.6) могут быть представлены в идентичном виде как f = f (0) (1 + P 3 (i )), где P 3 (i ) полином третьей степени. При этом коэффициенты по линомов для обеих функции распределения оказываются близкими.

3.2. Кинетический вывод КГД уравнений Регуляризованное кинетическое уравнение (3.4) построено здесь с использованием большого числа допущений, однако оно имеет це лый ряд свойств, сходных со свойствами БГК уравнения. В частно сти, для стационарных течений показано, что если функция распре деления f удовлетворяет БГК уравнению (3.5), то она удовлетворя ет и уравнению (3.4) с точностью O( 2 ). И наоборот, если функция распределения f удовлетворяет регуляризованному кинетическому уравнению, то с точностью O( 2 ) она удовлетворяет и стационарно му уравнению БГК. Таким образом, решения этих уравнений долж ны быть близки. Для уравнения (3.4) доказан аналог H-теоремы Больцмана [114].

На основе регуляризованного кинетического уравнения (3.4) уда ется построить КГД систему уравнений и родственные ей моментные уравнения для описания слабо-неравновесных газодинамических те чений (главы 8 и 9).

3.2 Кинетический вывод КГД уравнений Приведем способ построения квазигазодинамической системы урав нений на основе регуляризованного кинетического уравнения, вы писанного в предыдущем параграфе.

Последовательно интегрируя кинетическое уравнение (3.4) с ве сами 1,, 2 /2 и используя выражения f (0) d, = f d = u2 2 2 (0) + = f d = f d, c = u, 2 2 c2 c2 (0) (c c)f (0) d, (3.7) RT = f d = f d, Ip = 2 2 cf (0) d = 0, (c c c)f (0) d = 0, c2 (0) 5 p cc2 f (0) d = 0, (c c) f d = I, 2 62 Гл.3. Квазигазодинамические уравнения получим систему КГД уравнений в следующем виде:

(3.8) + div(u) = div div(u u) + p, t (u) + div(u u) + p = div { [div(u u u)+ t +( pu) + ( pu)T (3.9) + { [div(pu)]}, u2 u + + div u + + pu = t 2 u2 p = div div ++2 uu + 2 u2 p (3.10) + p ++.

2 Замыкание вида (3.7), основанное на локально-максвелловской функции распределения, автоматически приводит к построению си стемы моментных уравнений для газа твердых сфер, то есть для = 5/3. Обобщение на случай произвольных молекул проводится путем выделения слагаемых с удельной внутренней энергией в пред положении = 5/3, и затем обобщения этого выражения на случай произвольного значения.

Система (3.8)–(3.10), выведенная для случая произвольных кри волинейных эйлеровых координат, замыкается уравнениями состо яния идеального политропного газа p p = RT, =.

( 1) При записи уравнения энергии в последнем слагаемом имеются члены, не содержащие операций дивергенции от скорости. Из это го слагаемого можно сразу выделить слагаемые в форме теплового 3.3. КГД уравнения в форме законов сохранения потока Навье–Стокса qN S и ввести число Прандтля P r = E + i ui (E + p) = i (j (E + 2p)ui uj + i uk uk p) + t p 1 p (3.11) ( p) + P r ( p ).

1 1 Пример построения КГД уравнений для плоского одномерного течения газа детально выписан в Приложении A, где наглядно про слежены все использованные преобразования.

На основе осреднения регуляризованного кинетического уравне ния с учетом внешней силы F f + [( · x ) + (F · )]f (0) t f (0) f ( · ) + (F · ) ( · ) + (F · ) f (0) =. (3.12) Ю.В. Шеретов построил КГД уравнения для описания течений газа в поле внешних сил [113, 114].

Левые части системы (3.8) – (3.10) представляют собой левые ча сти уравнений Эйлера. Правые части имеют дивергентный вид и пропорциональны малому параметру. При этом построенные урав нения (3.8) – (3.10) не имеют вид законов сохранения, то есть в них в явном виде не выделены тензор вязких напряжений и векторы плотности потока массы и теплового потока. Диссипативные слага емые входящие в уравнение Навье–Стокса также явно не выделены.

Преобразование построенной выше КГД системы к форме законов сохранения проведена в следующем разделе.

3.3 КГД уравнения в форме законов сохранения Представим КГД систему в виде законов сохранения (1.54)–(1.56), которые для случая F = 0 имеют следующий вид:

(3.13) + div jm = 0, t 64 Гл.3. Квазигазодинамические уравнения (u) (3.14) + div(jm u) + p = div, t u2 u (3.15) + + div jm + = div A div q.

t 2 Уравнение неразрывности и вектор плотности потока массы Сравнивая уравнения переноса массы (3.8) и (3.13), найдем вектор плотности потока массы jm = u div(u u) + p.

Или, обозначив добавку к скорости через (3.16) w= div(u u) + p, представим jm в виде, предложенном ранее в первой главе, и удоб ном для дальнейшего анализа уравнений jm = (u w).

Уравнение импульса и тензор вязких напряжений Сопоставляя уравнения импульса (3.9) и (3.14), найдем вид тензора вязких напряжений. Для этого представим div(jm u) = div(u u) div( div(u u) u) div( p u).

Тогда (3.9) можно переписать в виде (u) + div(jm u) + p = div( div(u u) u) t div( p u) + div { [div(u u u)+ +( pu) + ( pu)T + div {I div(pu)}, (3.17) где использовано тождество (3.18) div {I div(pu)} = { div(pu)}.

3.3. КГД уравнения в форме законов сохранения Из (3.17) и (3.14) следует, что тензор вязких напряжений имеет вид = div(u u) u (p u) + [div(u u u)+ +( pu) + ( pu)T + I div(pu).

Представим этот тензор в виде суммы тензора вязких напряжений Навье–Стокса N S = µ[( u) + ( u)T (2/3) I div u] + I div u (3.19) и некоторой добавки. Для определения вида этой добавки восполь зуемся тождествами ( pu) = p( u) + (u p), ( pu)T = p( u)T + (p u), div(pu) = (u · )p + p div u, div(u u u) = div(u u) u + u (u · )u.

Тогда примет вид = p ( u) + ( u)T + u (u · )u + p + + I (u · )p + p div u.

Преобразуем полученное выражение, прибавляя и вычитая величи ны Ip div u и (2/3) pI div u = p ( u) + ( u)T (2/3)I div u + + u (u · )u + p + I (u · )p + p div u + + Ip div u Ip div u + (2/3) pI div u. (3.20) Группируя слагаемые, получаем:

66 Гл.3. Квазигазодинамические уравнения = p ( u) + ( u)T (2/3)I div u + + pI(5/3 ) div u + u (u · )u + p + + I (u · )p + p div u. (3.21) Сравнивая (3.21) с видом тензора Навье–Стокса (3.19), увидим, что (3.22) µ = p, = p.

Из первого соотношения сразу следует, что имеет смысл максвел ловского времени релаксации [134].

Из полученной формулы для коэффициента объемной (второй) вязкости следует, что коэффициент неотрицателен и связан с на личием внутренних степеней свободы молекулы, что соответствует теоретическим представлениям, изложенным в [18, 57, 72, 78]. Дей ствительно, для одноатомного газа = 5/3 и = 0, в противном случае, при наличии колебательных или вращательных степеней свободы молекулы, 5/3 и 0.

В результате этих преобразований получаем, что тензор вязких напряжений в КГД системе уравнений представляется в виде тензо ра вязких напряжений Навье–Стокса N S, (3.19) и добавки, которая пропорциональна и обращается в нуль при u = = N S + u (u · )u + p + I (u · )p + p div u. (3.23) Уравнение полной энергии и вектор теплового потока В заключение рассмотрим третью пару уравнений (3.10), (3.15), из которых определим вид векторов A и q. Для этого запишем (3.10) в 3.3. КГД уравнения в форме законов сохранения виде u2 u p + + div jm ++ = t 2 2 u2 p = div w ++ + 2 u2 p + div div ++2 uu + 2 u2 p + p ++, 2 а (3.15) перепишем в виде u2 u p + + div jm ++ = t 2 2 p = div jm + div A div q.

Сравнивая последние два выражения, найдем u p p A q = jm div(u u) + p ++ + 2 u2 u p p + div ++2 u u + p ++ + 2 2 u2 p + ++ p, 2 здесь для w использовано представление (3.16). Приводя подобные члены, получим u p p A q = jm ++ div(u u)+ 2 u2 p + ++2 div(u u)+ 2 u2 u p p + u(u · ) ++2 + p ++, 2 2 68 Гл.3. Квазигазодинамические уравнения где мы воспользовались тождеством u2 p div ++2 uu = 2 u2 u p p = ++2 div(u u) + u(u · ) ++2.

2 2 Приведем подобные слагаемые u2 p u p p A q = jm + ++2 div(u u)+ 2 2 u2 u p p + u(u · ) ++2 + p ++.

2 2 Применяя тождество div(u u) = u div(u) + (u · )u = = u div u + u(u · ) + (u · )u, (3.24) получим p p A q = jm + pu div u + u(u · ) + p(u · )u+ u2 + u(u · ) + u(u · ) + 2 pu(u · ) + 2 u2 p + 2 u(u · )p + p + p + p.

2 Сгруппируем слагаемые u p A q = jm + p + (u · )u + u + u (u · ) + (u · )p + u(u · )p+ p + pu div u + u(u · ) + pu(u · ) + 1 p + u (u · ) + p(u · ) + p +.

3.3. КГД уравнения в форме законов сохранения Применив тождества u ( u)u =, ( u)u = (u · )u, получим p A q = jm + p ( u) + ( u) 2/3I div u u+ + u u · (u · )u + p + u (u · )p + p div u + p + pu (5/3 ) div u + p + + + u (u · ) + p(u · ) + p + u(u · ) up 2 (u · ).

Из полученного выражения выделим тензор в виде (3.21) p p A q = u jm + p + +. (3.25) + u (u · ) + p(u · ) В теории Навье–Стокса [72] AN S = N S u pu представляет собой работу поверхностных сил давления и внутрен него вязкого трения в единицу времени. По аналогии положим jm A = u p.

70 Гл.3. Квазигазодинамические уравнения Остальные слагаемые в (3.25) представим как тепловой поток и обо значим через p q = p + u (u · ) + p(u · ).

Выделим в тепловом потоке q часть, связанную с потоком Навье– Стокса qN S = T.

Примем в качестве уравнений состояния уравнения состояния идеального политропного газа RT p = RT, =, представим формулу для теплового потока в виде R q = p T u (u · ) + p(u · ).

1 Сопоставляя формулу для q и qN S, получим, что коэффициент теп лопроводности в КГД модели определяется как R = p, и выражение для теплового потока q в КГД уравнениях имеет вид (3.26) q = qN S u (u · ) + p(u · ).

Таким образом, полученная на основе кинетического уравнения (3.4) КГД система вида (3.8)–(3.10) представлена в виде законов сохранения (3.13)–(3.15). При этом вектор плотности потока мас сы, тензор вязких напряжений и вектор теплового потока представ ляются в виде суммы соответствующих величин в форме Навье– Стокса и малых добавок существенно нелинейного вида, пропорци ональных коэффициенту (3.16, 3.23, 3.26).

3.4. Коэффициенты диссипации 3.4 Коэффициенты диссипации 3.4.1 Формулы для диссипативных коэффициентов и их обобщения На основании кинетического вывода КГД уравнений и представле ния их в виде законов сохранения удается получить вид диссипатив ных коэффициентов µ, и. Эти коэффициенты связаны между собой через параметр и выражаются как 5 R (3.27) µ = p, = p, = p.

3 Именно в таком виде коэффициенты вязкости µ и теплопровод ности получаются при выводе системы уравнений Навье–Стокса методами Чепмена–Энскога из уравнения БГК. Коэффициент объ емной вязкости, совпадающий с формулой (3.27), был получен в [193] на основе БГК приближения для немоноатомного газа, обла дающего вращательными степенями свободы. Действительно, при веденная в [193] формула для этого коэффициента имеет вид 2i =µ, 3(i + 3) где i число вращательных степеней свободы. Выражая общее чис ло степеней свободы n = i + 3 через как n = i + 3 = 2/( 1), получим коэффициент объемной вязкости вида (3.27).

Вводя число Прандтля P r = 1, запишем коэффициент теплопро водности в виде R =µ.

1 Pr В п. 1.5.3 с точностью до численного коэффициента была найдена связь релаксационного параметра и коэффициента динамической вязкости µ в виде µ =, pSc где Sc число Шмидта. Для газа его величина близка к единице.

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

72 Гл.3. Квазигазодинамические уравнения 3.4.2 Коэффициент объемной вязкости Диссипативные эффекты, связанные с релаксацией внутренних сте пеней свободы в газе, могут оказывать заметное влияние на течения с ударными волнами и на быстропеременные по времени процессы.

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

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

Способы вычисления коэффициентов второй вязкости и получа ющиеся при этом выражения в общем случае достаточно сложны (см., например, [18, 78, 102]).

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

Пример такого выражения для приведен в предыдущем разделе.

Более точное выражение для коэффициента второй вязкости, за висящей от вращательных степеней свободы, предложено в [57] в виде:

pR (3.28) = rot rot, cV здесь rot доля внутренней энергии, содержащейся во вращатель ных степенях свободы.

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

i rot =, n 3.4. Коэффициенты диссипации где i число вращательных степеней свободы, n общее число степеней свободы. Если принять во внимание, что n = i + 3, где 3 число поступательных степеней свободы, то формулу для rot можно переписать в виде n rot =.

n Учитывая связь числа степеней свободы с показателем адиабаты n= получим, что (3.29) rot = ( ).

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

rot = Zrot c, где c среднее время свободного пробега молекул;

Zrot коэф фициент обмена энергий между вращательными и поступательны ми степенями свободы, который равен среднему числу межмолеку лярных столкновений, необходимых для релаксации к равновесно му состоянию вращательной моды. Выражения для Zrot приведены, например, в [57, 134]. Эти выражения для Zrot использовались для расчетов неравновесных течений в струях CO2 [172] и N2 [159], а также для численного моделирования структуры ударной волны в азоте [55]. Например, для азота Zrot = Z /[1 + ( 3/2 /2)(T /T )1/2 + ( + 2 /4)(T /T )], где Z = 23, T = 91.5K, если температура газа в невозмущенном потоке составляет 273K [57]. В этом случае значение Zrot изменяется в пределах 5 16 в области температур от 300 до 6000 K (рис. 3.2).

Выразим c через динамическую вязкость газа µ. Cредняя длина свободного пробега частиц в газе определяется средним временем свободного пробега частиц и средней тепловой скоростью c как (3.30) = c c.

74 Гл.3. Квазигазодинамические уравнения Величина средней тепловой скорости частиц в равновесном газе составляет (п. 2.5) c = 8RT /. Таким образом, (3.31) c =.

8RT / Средняя длина свободного пробега может быть связана с вязкостью газа как µ (3.32) =A RT, p где A= формула Чепмена [3] или 2(7 2)(5 2) A= 15 формула Берда [134]. Здесь показатель степени в законе вяз кости µ T.

Соотношения (3.31) и (3.32) позволяют выразить c через динами ческую вязкость газа µ в виде µ (3.33) c = A.

p Подставив выражения для cv, rot (3.29) и c (3.33) в (3.28), запишем исходную формулу для второй вязкости (3.28) в виде 5 3 (3.34) = µ( )( 1)A Zrot = µ( )B.

3 2 8 Последняя формула (3.34) с точностью до численного множителя 3 B = ( 1)A Zrot 2 совпадает с формулой для коэффициента второй вязкости (3.27), полученной ранее при анализе КГД уравнений. На рис. 3.2 приве дены значения B (пунктир) и Zrot (сплошная линия) в зависимости 3.5. Система Навье–Стокса как асимптотика КГД системы Рис. 3.2. Коэффициенты B (пунктир) и Zrot (сплошная линия) для азота от температуры для азота. Видно, что коэффициент B изменяется в пределах от 1 до 4.

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

3.5 Система Навье–Стокса как асимптотика КГД системы Диссипативные слагаемые N S и qN S, входящие в КГД уравнения, при условии достаточной гладкости соответствующих производных, имеют порядок малости O(µ) O( ), или, в безразмерном виде, O(Kn). Добавки к тензору вязких напряжений и векторам плотно сти потока массы и теплового потока, выписанные в п. 3.3, имеют такой же порядок малости. Покажем, что для стационарных течений эти КГД добавки имеют уже не первый, а второй порядок малости 76 Гл.3. Квазигазодинамические уравнения по, то есть O( 2 ).

В п. 3.3 величины jm,, A и q в КГД системе были выписа ны в виде суммы соответствующих слагаемых типа Навье–Стокса и КГД добавок. Для удобства изложения приведем отдельно вид этих добавок:

(3.35) jQGD = div(u u) + p, (3.36) QGD = u (u · )u + p + I (u · )p + p div u, p (3.37) AQGD = QGD u jQGD, (3.38) qQGD = u (u · ) + p(u · ).

единичный тензор.

Здесь I Покажем, что в стационарном случае выписанные величины формально имеют второй порядок малости O( 2 ).

Для доказательства используем подход, изложенный в [114]. За пишем систему (3.8)–(3.10) в стационарном случае:

(3.39) div(u) = M, (3.40) div(u u) + p = I, u (3.41) div u + + pu = E.

Здесь правые части соответсвующих уравнений обозначены как M = div div(u u) + p, I = div div(u u u) + ( pu) + ( pu)T + + { [div(pu)]}, u2 p E = div div ++2 uu + 2 u p + p ++.

2 3.5. Система Навье–Стокса как асимптотика КГД системы Сразу заметим, что (3.42) M O( ), I O( ), E O( ).

КГД добавка к вектору потока массы Из (3.35), (3.40) и (3.42) следует, что jQGD = I O( 2 ).

КГД добавка к тензору вязких напряжений Применяя тождество (3.24), преобразуем (3.40) (3.43) u div(u) + (u · )u + p = I.

Тогда первое слагаемое в QGD запишется в виде (3.44) (u · )u + p = I u div(u) = I uM.

p Раскроем дивергенцию в (3.41), учитывая = ( 1) u2 u div u + + pu = div(u) + u(u · )u+ 2 + (u · )p + p div u.

Откуда следует, что u (u · )p + p div u = ( 1) E M u(u · )u.

Далее прибавим и вычтем (u · )p (u · )p + p div u = ( 1) u M u(u · )u (u · )p. (3.45) E Преобразуем выражение (u·)p. Для этого умножим скалярно (3.43) на u. В результате получим 78 Гл.3. Квазигазодинамические уравнения (u · )p = (u · I) u2 M u(u · )u.

Итак, учитывая (3.5), запишем (3.45) следующим образом u (u · )p + p div u = ( 1) E M (u · I).

Преобразуя (3.44) с учетом (3.45), окончательно получим u QGD = u (I uM ) + I( 1) E M (u · I).

Из последнего равенства следует, что QGD O( 2 ).

КГД добавка к вектору работы сил давления и вязкого тре ния Из (3.35), (3.36) и (3.37) вытекает, что AQGD O( 2 ).

КГД добавка к вектору теплового потока Преобразуем левую часть выражения (3.41) u E= M + u(u · )u + M + (u · ) + (u · )p + p div u и выразим u (u · ) + p div u = E M + u(u · )u (u · )p.

Из (3.39) выразим div u как M (u · ) M (3.46) div u = = + (u · ).

Учитывая (3.5) и (3.46), получим u2 p (u · I). (3.47) (u · ) + p(u · ) =EM 2 3.6. Модель для течений с внешними источниками Из (3.47) и (3.38) следует, что u2 p qQGD = u E M (u · I).

2 Из последнего соотношения заключаем, что qQGD O( 2 ).

Результаты, полученные в этом разделе, означают, что для стационарных течений КГД уравнения аппроксимируют систему Навье–Стокса с точностью O( 2 ), или, в безразмерном виде, с точностью O(Kn2 ) при достаточной гладкости соответствующих пространственных производных. Таким образом, можно полагать, что при стремлении величины к нулю, решение КГД уравнений с соответствующим порядком сходится к решению системы урав нений Навье–Стокса. Заметим, что для нестационарных течений аппроксимацию второго порядка доказать не удается.

Другой важной асимптотикой уравнений гидродинамики явля ется приближение пограничного слоя [72, 79]. Было показано (см., например, [51, 105, 114]) что приближением пограничного слоя для КГД уравнений, так же как и для уравнений Навье–Стокса, слу жат уравнения Прандтля. Вид уравнений Прандтля для двумерных течений приведен в приложении B.

3.6 КГД уравнения для течений газа при наличии внешних сил и тепловых источников В этом разделе приведен другой способ построения КГД системы. С его помощью построены уравнения для описания течений в присут ствии внешних сил и источников тепла. Полученные системы можно использовать для численного моделирования целого ряда актуаль ных задач, например, задачи горения, исследования возможностей управления потоками с помощью энерговложения, расчетов актив ных сред в резонаторах газовых лазеров, течений излучающего или поглощающего газа [105].

Для построения КГД системы следуя методике, изложенной в [114], выпишем систему уравнений Эйлера для идеального полит ропного газа 80 Гл.3. Квазигазодинамические уравнения (3.48) + ui = 0, t xi (3.49) ui + ui uj + p = Fi, t xj xi u2 u p (3.50) + + ui ++ = ui Fi + Q t 2 xi 2 с уравнениями состояния R (3.51) p = RT, = T, и следующие из нее дифференциальные тождества 1 1 (3.52) + ui ui = 0, t xi xi (3.53) ui + uj ui + p Fi = 0, t xj xi p Q (3.54) + ui + ui = 0, t xi xi (3.55) p + ui p + p ui ( 1)Q = 0.

t xi xi Здесь и далее использованы обычные обозначения: удельная плотность, ui компонента скорости, p давление, Fi ком понента внешней силы, удельная внутренняя энергия, Q мощность тепловых источников. В формулах подразумевается сум мирование по повторяющимся индексам. Справедливость тождеств (3.52) – (3.55) можно проверить непосредственной подстановкой в них выражений для производных по времени из уравнений Эйлера (3.48) – (3.50) с последующим приведением подобных слагаемых.

Запишем интегральные законы сохранения для малого непо движного объема V в конечно-разностном виде, заменяя производ ные по времени разностным отношением для моментов t и t+t, где 3.6. Модель для течений с внешними источниками t малый промежуток времени. Тогда законы сохранения массы, импульса и полной энергии можно представить, соответственно, в виде:

u di = 0, (3.56) d x+ i t V ui ui u u dj + p di = d x+ ij t V Fi d3 x + ij dj, (3.57) = NS V (2 /2 + ) (u2 /2 + ) u d x+ t V (u )2 p u + + di + qNS i di = i 2 u Fi d3 x + ij u dj + Q d3 x. (3.58) = i i NS V V Здесь вектор теплового потока и тензор вязких напряжений вы числяются как:

(3.59) qNS i = T, xi (3.60) NS ij = µ uj + ui ij uk, xi xj 3 xk µи коэффициенты вязкости и теплопроводности. Величинами со звездочками в (3.56) – (3.58) обозначены значения газодинамиче ских параметров в момент времени t, где t t t+t. Обозначим t/2 =, и определим параметры газа в средней точке t = t, огра ничиваясь первым порядком малости по :

= +, (3.61) t ui u = ui +, i t p p = p +.

t 82 Гл.3. Квазигазодинамические уравнения Подставим эти выражения в формулы (3.56) – (3.58). Величины Fi и Q считаем мало меняющимися за время, поэтому звездочки возле них опустим. В получающихся формулах пренебрежем также сла гаемыми порядка O( 2 ), O( · µ) и O( · ). Возвращаясь вновь к дифференциальному виду производной по времени в интегралах по объему, получим d3 x + (3.62) ui + ui di = 0, t t V ui d3 x + ui + ui uj dj + t t V + ui uj dj + p+ p di = t t Fi d3 x + NS ij dj, (3.63) = + t V u + d3 x+ t V u p + ui + ui · ++ di + t 2 1 + ui uj uj + +p + p di + t t t t ui Fi d3 x+ + qNS i di = ui + t V Qd3 x. (3.64) + NS ij ui dj + V Считаем, что в нулевом по приближении для рассматриваемого газа справедливы уравнения Эйлера (3.48) – (3.50). Используем эти уравнения и следующие из них тождества (3.52) – (3.55) для того, чтобы исключить производные по времени в слагаемых, линейных по. Приведем подобные члены, сгруппируем однотипные слагае мые и используем произвольность объёма V, чтобы перейти от ин тегральной формы уравнений (3.62) – (3.64) к дифференциальной.

3.6. Модель для течений с внешними источниками В результате получим систему квазигазодинамических уравне ний в виде законов сохранения с учетом влияния внешних сил и тепловых источников:

(3.65) + jmi = 0, t xi (3.66) ui + jmj ui + p = Fi + ji, t xj xi xj u2 u p + + jmi ++ + qi = t 2 xi 2 xi ij uj + Q, (3.67) = jmi Fi + xi где введены обозначения:

(3.68) jmi = (ui wi ), (3.69) wi = ui uj + p Fi, xj xi (3.70) = uk, xk ij = NS ij + ui uk uj + p Fj + xk xj uk ( 1)Q, (3.71) + ij uk p + p xk xk 1 Q (3.72) qi = qNS i ui uj + puj.

xj xj 84 Гл.3. Квазигазодинамические уравнения 3.7 Уравнение баланса энтропии В работе [114] было показано, что в случае Q = 0 производство эн тропии для системы (3.65) – (3.67) является неотрицательным. Обоб щим этот результат на случай Q = 0.

Обозначим полную производную по времени оператором D· = · + (ui wi ) · = ·+ jmi · t xi t xi Будем считать, что параметры рассматриваемого нами газа удо влетворяют тождеству Гиббса T Ds = D + pD.

Используя систему (3.65) – (3.67), получим из этого равенства урав нение для энтропии s. Перенесем в уравнении (3.66) все слагаемые в левую часть и домножим его на ui :

0 = ui ui + jmj ui + p Fi ji = t xj xi xj u2 = D + ui p Fi ui ui ji.

2 xi xj Отсюда следует соотношение u2 D = ui p + Fi ui + ui ji.

2 xi xj Уравнение (3.67) можно записать в виде u2 p D + + jmi + qi = jmi Fi + ij uj + Q, 2 xi xi xi откуда, используя предыдущее равенство, получаем 3.7. Уравнение баланса энтропии p D = jmi Fi + ji ui qi jmi xj xi xi Fi ui ui ji + ui p = (jmi ui ) Fi + xj xi p + ji ui qi jmi + ui p + Q.

xj xi xi xi Заметим, что 1 D = (ui wi ), xi поэтому 1 p Ds = D + D = T T 1 p = (jmi ui ) Fi + ji ui qi jmi + T xj xi xi p qi + ui p+Q + (ui wi ) = + qi + xi T xi xi T xi T 1 + (jmi ui ) Fi + ji ui T xj p jmi + ui p+Q+p (ui wi ).

xi xi xi Последнее соотношение представляет собой уравнение баланса эн тропии вида qi Ds = + X, xi T где производство энтропии X представлено выражением 1 1 X = qi + (jmi ui ) Fi + ji ui xi T T xj p (ui wi ). (3.73) jmi + ui p+Q+p xi xi xi 86 Гл.3. Квазигазодинамические уравнения Критерием физической корректности гидродинамической моде ли является неотрицательность производства энтропии X 0 при 0. Покажем, что указанное свойство для КГД модели выпол Q няется. Для этого преобразуем выражение (3.73), учитывая (3.68) – (3.72). Рассмотрим отдельно каждое слагаемое, выделяя, где это воз можно, неотрицательные комбинации величин.

Для первого слагаемого получим:

1 qi = T ui uj + xi T xi xj T 1 Q + puj = + xj xi T T 1 Q + ui uj + puj T= T 2 xi xj xj ui uj T = + T T xi xj ui uj p ui Q.

T xi xj T xi Преобразуем второе слагаемое 1 1 ((jmi ui )Fi ) = jmi ( uj )ui Fi = T T xj Fi = wi + ui uj = T xj Fi = ui uj ui uj p + Fi = T xj xj xi Fi Fi = Fi Fi uj ui p.

T T xj T xi Заметив, что NS ij NS ij µ ui ui 2 uk um ui uj = +, 2µT T xj xj 3 xk xm xj xi 3.7. Уравнение баланса энтропии преобразуем следующий фрагмент в формуле для X:

1 1 2 ji ui = µ ui + uj ij uk + T xj T xj xi 3 xk + uj uk ui + p Fi + xk xi + ij uk p + p uk ( 1)Q ui = xk xk xj NS ij NS ij uj 1 = + uk ui + p Fi ui + 2µT T xk xi xj + uk p + p uk ( 1)Q ui.

T xk xk xi Наконец, рассмотрим оставшиеся слагаемые:

1 p 1 p jmi + ui p+ (ui wi ) = T xi T xi T xi = ui uj + p Fi p.

T xj xi xi Собирая все фрагменты исходной формулы вместе, получаем:

T 2 NS ij NS ij X= + + uk ui + p Fi + T 2µT T xk xi ui uj pui uj + T xi xj T xi xj ui Q + uk p + p uk ( 1)Q ui + T xi T xk xk xi Q + ui p uj + ui uj p +.

T xi xj T xi xj T Для краткости здесь и далее использовано обозначение (ai ) ij ai aj.

Используя соотношение p = ( 1), 88 Гл.3. Квазигазодинамические уравнения можно записать предыдущее равенство в виде 2 NS ij NS ij T X= + + uk ui + p Fi + T 2µT T xk xi 2 p Q p Q, (3.74) + ui + ui +2 ui + T xi xi T xi T где 1 Q T = T · 1 + ui + ( 1) ui.

xi xi Два слагаемых равенства (3.74), содержащих величину Q, можно переписать в виде p 1Q Q ( 1) Q ui + ui + 1.

T xi xi 2 T Как видно из полученного выражения, при Q 0 (при Q 0 эн тропия может убывать), производство энтропии является неотрица тельным, если величина является достаточно малой ( 1)Q 1.

4p При = 0 производство энтропии (3.74) совпадает с соответствую щей величиной, вычисленной для уравнений Навье–Стокса.

Это преобразование было выполнено А.А.Злотником Глава КГД уравнения и системы координат В этой главе квазигазодинамические уравнения выписаны в тензор но-индексном представлении в произвольных ортогональных коор динатах, что позволяет использовать КГД систему в удобной для решения задачи системе координат. Выписанный вид уравнений рас ширяет класс решаемых задач на сложные, в том числе трехмерные, расчетные области, в которых координатная сетка может задаваться как аналитически, так и численно. Это позволяет строить однород ные разностные схемы на квазиортогональных пространственных сетках в преобразованном пространстве координат [100]. В двух по следних параграфах КГД уравнения выписаны в декартовой и ци линдрической системах координат [45].

4.1 КГД уравнения в произвольной системе координат В соответствии с [29,66,81] приведем формулы тензорного анализа, которые понадобятся для дальнейших выкладок.

Рассмотрим в евклидовом пространстве некоторый базис ej, лю бой вектор в котором задается координатами (x1, x2, x3 ). Пусть ri новый базис, любой вектор в котором задается координатами (x1, x2, x3 ).

Определим переход от базиса ej к новому базису ri с помощью следующего преобразования (здесь и далее подразумевается сумми рование по повторяющимся индексам):

xj ri = bj ej, bj = (4.1), i i xi где bj матрица перехода от базиса ej к базису ri.

i 90 Гл.4. КГД уравнения и системы координат В дифференциальной геометрии и тензорном анализе вводится понятие метрического тензора, который отражает дифференциаль ные свойства внутренней геометрии области. Компоненты метриче ского тензора определяются как:

xk xk (4.2) gij =, xi xj g det(gij ) определитель метрического тензора.

Рассмотрим некоторый вектор U = U i ei в базисе ei, где U i координаты вектора в данном базисе. При переходе к новому нор мированному базису (ei = ri / |gii |) координаты вектора U преоб разуются по следующей формуле:

Uj U j = U i j Ui = bi, (4.3) bi |gjj |, j |gjj | где = b1 матрица обратного перехода от базиса ri к базису ej.

b Координаты вектора U i в нормированном базисе ej связаны с координатами вектора U i в базисе ri следующим соотношением:

Ui = |gii |U i. (4.4) Аналогичная формула для компонент тензора ik имеет вид:

ik = |gii | |gkk |ik. (4.5) В формулах (4.4) и (4.5) суммирование по повторяющимся ин дексам не производится.

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

Скалярное произведение векторов:

(A · B) = gkl Ak B l. (4.6) Ковариантная и контравариантная производная скалярной функции:

A i A = gij j A. (4.7) i A =, xi 4.1. КГД уравнения в произвольной системе координат Ковариантная производная тензора первого и второго ранга:

U i j U i = + U m i, (4.8) mj xj Aij + i Alj + j Aim, i Aij = (4.9) li mi xi где glj gij 1 gil k = gkl (4.10) + ij i j xl 2 x x символы Кристоффеля второго рода.

Дивергенция вектора (частный случай (4.8) при i = j):

1 ( |g|Ai ) i Ai = (4.11).

xi |g| В дальнейшем мы рассматриваем ортогональные системы ко ординат, в которых элементы метрических тензоров подчиняются i i условию gik = gii k, где k символ Кронекера, для которого 1, если i = k i k =.

0, если i = k В системе ортогональных координат для вычисления символов Кристоффеля k справедливы следующие соотношения:

ij k = 0 (4.12) (i = j, j = k, k = i), ij 1 gii k = (4.13) (i = k), ii 2gkk xk 1 gkk 1 (ln gkk ) k = k = (4.14) =.

ik ki i 2 xi 2gkk x Выпишем систему КГД уравнений в тензорно-индексном пред ставлении.

Уравнение баланса массы имеет вид:

1 i (4.15) + ( |g|jm ) = 0, |g| xi t 92 Гл.4. КГД уравнения и системы координат где вектор плотности потока массы определяется как Ui i wi jm =, |gii | Uj Ui Ul Ui + j wi = + lj xj |gjj | |gii | |gll | |gii | Um Uj p i + gij j.

mj x |gjj | |gmm | Уравнение баланса импульса записывается следующим образом Ui Ui Ui + j j l + jm jm + lj xj t |gii | |gii | |gii | Um ji p + j li + i jm, (4.16) i j + gij jm = mj mj lj xj xj |gmm | а компоненты тензора вязких напряжений определяются по форму лам (U i / |gii |) Uj Uk Um p ji = i ik + mk + g + xk xk |gjj | |gkk | |gmm | Uk Um p 1 + ji S, gji + p |g| N k |g| xm |gkk | x |gmm | (U i / |gii |) Um ji S =µ gjk i + + mk N xk |gmm | (U j / |gjj |) Um j µ gik + |gmm | mk xk Um 1 2/3gji |g|.

|g| xm |gmm | 4.2. Декартова система координат Уравнение баланса энергии имеет вид:

E 1 E+p 1 j |g|q j = + |g|jm + |g| xj |g| xj t Ul 1 |g|gkl jk ), (4.17) ( |g| xj |gll | где энергия и тепловой поток записываются как Uk Ul gkl 1p E = +, 2 |gkk | |gll | µ p qj = gjk k Pr 1 x Uj Uk Uk 1 p +p.

|gkk | xk |gkk | xk |gjj | 1 Тензорно-индексный вид (4.15) – (4.17) позволяет адаптиро вать КГД систему уравнений к произвольным ортогональным системам координат: декартовым, цилиндрическим, сферическим, эллипсоидальным, а также к координатным системам, построен ным численно.

4.2 Декартова система координат В декартовой системе координат положим x1 = x, x2 = y, x3 = z.

Метрический тензор имеет диагональный вид:

gij = 0 1 0, (4.18) g det gij = 1.

Формулы (4.6) – (4.11) принимают более простой вид:

94 Гл.4. КГД уравнения и системы координат Скалярное произведение векторов:

(A · B) = Ai Bi. (4.19) Прямое тензорное произведение векторов:

(A B) = Ai B j. (4.20) Производная скалярной функции (градиент):

p grad p = i p = i p = (4.21).

xi Производная тензора первого и второго ранга:

Ai j Ai = (4.22), xj ij i ij = (4.23).

xi Дивергенция:

Ax Ay Az div A = ( · A) = i Aj = (4.24) + +, x y z div = ( · ) = i ij. (4.25) Дивергенция от прямого тензорного произведения:

div(A B) = i (Ai B j ). (4.26) Дивергенция от скалярного произведения тензора и вектора:

div( · A) = ( · ( · A)) = i (ij Aj ). (4.27) Система КГД уравнений (4.15) – (4.17) запишется в виде:

i (4.28) + i jm = 0, t ( uj ) + i jm uj + j p = i ij, i (4.29) t u2 u p i + i q i = i (ij uj ).

+ + i jm + + t 2 2 (4.30) 4.2. Декартова система координат Выражения (4.4) для искомых компонент векторов jm и q i будут i вычисляться как:

1 2 jmx = jm, jmy = jm, jmz = jm, qx = q 1, qy = q 2, qz = q 3.

Система уравнений (4.28) – (4.30) в переменных x, y, z записыва ется как jmx jmy jmz (4.31) + + + = 0, t x y z (ux ) (jmx ux ) (jmy ux ) (jmz ux ) p + + + + = t x y z x xx yx zx, (4.32) = + + x y z (uy ) (jmx uy ) (jmy uy ) (jmz uy ) p + + + + = t x y z y xy yy zy, (4.33) = + + x y z (uz ) (jmx uz ) (jmy uz ) (jmz uz ) p + + + + = t x y z z yz xz zz, (4.34) = + + x y z E (jmx H) (jmy H) (jmz H) qx qy qz + + + + + + = t x y z x y z = (xx ux + xy uy + xz uz ) + (yx ux + yy uy + yz uz ) + x y (zx ux + zy uy + zz uz ). (4.35) + z Здесь компоненты вектора плотности потока массы jm имеют вид:

(u2 ) (ux uy ) (ux uz ) p x jmx = jm = ux + + +, x y z x (ux uy ) (u2 ) (uy uz ) p y, (4.36) jmy = jm = uy + + + x y z y 96 Гл.4. КГД уравнения и системы координат (ux uz ) (uy uz ) (u2 ) p z jmz = jm = uz + + +.

x y z z В формулах (4.36) первые слагаемые соответствуют вектору плотности потока массы Навье–Стокса, а последующие (с мно жителем в качестве коэффициента) представляют собой КГД добавки.

Элементы тензора вязких напряжений имеют вид:

xx = 11, xy = 12, xz = 13, yx = 21, yy = 22, yz = 23, (4.37) 31 32 zx =, zy =, zz =, ux xx = µ 2 div u + p div u+ x 3µ ux ux ux p p p + ux ux + uy + uz + 2ux + uy + uz, x y z x y z uy uy uy uy ux 1 p xy = µ + + ux ux + uy + uz +, x y x y z y uz ux uz uz uz 1 p xz = µ + + ux ux + uy + uz +, x z x y z z ux uy ux ux ux 1 p yx = µ + + uy ux + uy + uz +, y x x y z x uy yy = µ 2 div u + p div u+ y 3µ uy uy uy p p p + uy ux + uy + uz + ux + 2uy + uz, x y z x y z 4.2. Декартова система координат uz uy uz uz uz 1 p yz = µ + + uy ux + uy + uz +, y z x y z z ux uz ux ux ux 1 p zx = µ + + uz ux + uy + uz +, z x x y z x uy uy uy uy uz 1 p zy = µ + + uz ux + uy + uz +, z y x y z y uz zz = µ 2 div u + p div u+ z 3µ uz uz uz p p p + uz ux + uy + uz + ux + uy + 2uz.

x y z x y z В выражениях для компонент тензора вязких напряжений (4.37) сла гаемые с множителем µ соответствуют тензору Навье–Стокса, а слагаемые с множителем представляют собой КГД добавки.

Приведем далее более удобный и экономичный в численной ре ализации вариант записи компонент вектора плотности потока и тензора вязких напряжений для случая = 0:

jmx = (ux wx ), jmy = (uy wy ), jmz = (uz wz ), где (u2 ) (ux uy ) (ux uz ) p x (4.38) wx = + + +, x y z x (ux uy ) (u2 ) (uy uz ) p y (4.39) wy = + + +, x y z y 98 Гл.4. КГД уравнения и системы координат (ux uz ) (uy uz ) (u2 ) p z (4.40) wz = + + +.

x y z z Компоненты тензора вязких напряжений выписаны с помо щью вспомогательных величин, обозначенных индексом ( ) в виде:

ux xx = N S + ux · wx + R, N S = 2µ · µ · div u, xx xx x uy ux xy = N S + ux · wy, N S = N S = µ · +, xy xy yx x y uz ux xz = N S + ux · wz, N S = N S = µ · +, xz xz zx x z yx = N S + uy · wx, yx uy yy = N S + uy · wy + R, N S = 2µ · µ · div u, yy yy y uy uz yz = N S + uy · wz, N S = N S = µ · +, yz yz zy y z zx = N S + uz · wx, zx zy = N S + uz · wy, zy uz zz = N S + uz · wz + R, N S = 2µ · µ · div u, zz zz z а выражения для вспомогательных величин wx, wy, wz и R зада ются формулами:

ux ux ux p wx = ux + uy + uz +, x y z x uy uy uy p wy = ux + uy + uz +, x y z y uz uz uz p wz = ux + uy + uz +, x y z z p p p R = ux (4.41) + uy + uz + · p · div u.

x y z 4.3. Цилиндрическая система координат Компоненты вектора теплового потока q вычисляются как:

qx = q 1 = qx S ux · Rq, N qy = q 2 = qy S uy · Rq, N (4.42) qz S N q qz = q = uz · R, где 1 p p p Rq = (4.43) ux + uy + uz 1 x y z 1 1 + p ux + uy + uz, x y z Слагаемые Навье–Стокса qx S, qy S и qz S определяются по форму N N N лам:

T T T qx S = N qy S = N qz S = N,,.

x y z 4.3 Цилиндрическая система координат В цилиндрической системе координат положим x1 = r, x2 =, x3 = z.

Метрический тензор имеет диагональный вид:

gij = 0 r 2 0, (4.44) а его определитель равен g det gij = r 2.

Воспользуемся приведенными выше формулами тензорного ана лиза (4.4)–(4.11) и выпишем r,, z – компоненты вектора плотности потока массы jm, тензора вязких напряжений и вектора теплового потока q.

Компоненты вектора плотности потока массы jm :

100 Гл.4. КГД уравнения и системы координат 1 (ru2 ) 1 (ur u ) (ur uz ) p r jmr = jm = ur + + +, r r r z r 1 (u2 ) u 1 (ru ur ) jm = rjm = +2 + r r r r 1 (u uz ) 1 p, (4.45) + r z r 1 (ruz ur ) 1 (uz u ) (u2 ) p z jmz = jm = uz + + +.

r r r z z В каждой из вышеприведенных формул первое слагаемое пра вой части соответствует вектору потока массы Навье–Стокса jm S, а N последующие (с множителем в качестве коэффициента) представ QGD ляют собой КГД добавки jm.


Элементы тензора вязких напряжений имеют вид:

ur rr = 11 = µ 2 div u + r 3µ u ur u ur ur 1 p + ur ur + + uz + + r r r z r p u p p + ur +3 + uz + p div u, r r z u 1 u 1 ur r = r12 = µ 2+ + r r r r ur u u u 1 p r 2 ur u + u + rur + ruz +, r r z ur uz rz = 13 = µ + + z r uz u uz uz 1 p + ur ur + + uz +, r r z z 4.3. Цилиндрическая система координат u 1 u 1 ur r = r21 = µ 2+ + r r r r u ur u ur ur 1 p + u ur + + uz +, r r r z r r 1 u µ ur = r 2 22 = +2 div u + r2 r r 3 µ u u u u 1 p + rur + u + ur u + ruz + + r r z p u p p + ur +3 + uz + p div u, r r z 1 u 1 uz z = r23 = µ +2 + r z r u uz uz uz 1 p + u ur + + uz +, r r z z ur uz zr = 31 = µ + + z r u u ur ur ur 1 p + uz ur + + uz +, r r r z r 1 u 1 uz z = r32 = µ +2 + r z r uz u u u 1 p + rur + u + ur u + ruz +, r r z 102 Гл.4. КГД уравнения и системы координат uz zz = 33 = µ 2 div u + z 3µ u uz uz uz 1 p + uz ur + + uz + r r z z p u p p + ur + + uz + p div u, r r z В выражениях для компонент тензора вязких напряжений сла гаемые с множителем µ соответствуют тензору Навье–Стокса, а сла гаемые с множителем являются КГД добавками.

И, наконец, выпишем компоненты вектора теплого потока q:

qr = q 1 = qr S ur · Rq, N (4.46) = qz S N q qz = q uz · R, q 3 = q S N u · Rq, q = где u 1 p p p Rq = (4.47) ur + uz + 1 r z r 1 1 u + p ur + uz +, r z r Слагаемые Навье-Стокса qr S, qz S и q S определяются по форму N N N лам:

T T T qr S = N qz S = N q S = N,,.

r r z Глава Численные алгоритмы решения задач газовой динамики В этой главе приведены численные алгоритмы решения задач газо вой динамики, основанные на квазигазодинамических уравнениях.

Разностные аппроксимации строятся в потоковой форме непо средственно для векторов плотности потока массы jm, вектора теп лового потока q и тензора вязких напряжений. Это соответствует записи уравнений газовой динамики в виде законов сохранения и делает алгоритм достаточно компактным и экономичным. Устойчи вость разностного алгоритма обеспечивается присутствием искус ственной диссипации, вид которой определяется КГД добавками.

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

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

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

Последняя задача описана в приложении D.

В основу этого раздела положены материалы [44, 46, 53, 114, 116].

5.1 Система уравнений для плоских двумерных течений Приведем систему квазигазодинамических уравнений в декартовых координатах для двумерных течений:

jmx jmy (5.1) + + = 0, t x y 104 Гл.5. Алгоритмы решения задач газовой динамики (ux ) (jmx ux ) (jmy ux ) p xx yx (5.2) + + + = +, t x y x x y (uy ) (jmx uy ) (jmy uy ) p xy yy (5.3) + + + = +, t x y y x y E (jmx H) (jmy H) qx qy + + + + = t x y x y (yx ux + yy uy ), (5.4) = (xx ux + xy uy ) + x y где E полная энергия единицы объема и H полная удельная энтальпия:

u2 + u2 p (E + p) x y (5.5) E= +, H=, p = RT.

2 1 Компоненты вектора плотности потока массы jm вычисляются по формулам:

(5.6) jmx = (ux wx ), jmy = (uy wy ), где (u2 ) (ux uy ) p x wx = + +, x y x (ux uy ) (u2 ) p y (5.7) wy = + +.

x y y Компоненты тензора вязких напряжений определяются с по мощью выражений, экономичных с точки зрения программной реа лизации:

ux xx = N S + ux wx + R, N S = 2µ µ div u, xx xx x uy ux xy = N S + ux wy, N S = N S = µ +, xy xy yx x y yx = N S + uy wx, (5.8) yx uy yy = N S + uy wy + R, N S = 2µ µ div u, yy yy y 5.1. Система уравнений для плоских двумерных течений здесь N S, N S, N S, N S компоненты навье-стоксовского тен xx xy yx yy зора вязких напряжений.

Величины wx, wy, R вычисляются по формулам:

ux ux p wx = ux + uy +, x y x uy uy p (5.9) wy = ux + uy +, x y y p p R = ux + uy + pdivu, x y а дивергенция вектора скорости div u равна:

ux uy (5.10) div u = +.

x y Компоненты вектора теплового потока q имеют вид:

qx = qx S ux Rq, N qy = qy S uy Rq, N (5.11) uy ux p p 1 Rq = + + pux + puy, 1 x 1 y x y где навье-стоксовские слагаемые qx S и qy S вычисляются как:

N N T T qx S = N qy S = N (5.12),.

x y Зависимость коэффициента динамической вязкости µ от темпе ратуры выберем в следующем виде T (5.13) µ = µ, T где µ = µ(T ) известное значение µ при температуре T.

Коэффициент теплопроводности и релаксационный параметр связаны с коэффициентом динамической вязкости µ соотношени ями:

R (5.14) = µ, = µ, ( 1)P r p Sc где P r число Прандтля, Sc число Шмидта.

106 Гл.5. Алгоритмы решения задач газовой динамики 5.2 Система уравнений в цилиндрической геометрии В цилиндрических координатах в двумерном случае система КГД уравнений имеет вид:

1 (rjmr ) jmz (5.15) + + = 0, t r r z (ur ) 1 (rjmr ur ) (jmz ur ) p 1 (rrr ) zr + + + = +, t r r z r r r z r (5.16) (uz ) 1 (rjmr uz ) (jmz uz ) p 1 (rrz ) zz, (5.17) + + + = + t r r z z r r z E 1 (rjmr H) (jmz H) 1 (rqr ) qz + + + + = t r r z r r z 1 (zr ur + zz uz ), (5.18) = [r (rr ur + rz uz )] + r r z где E полная энергия единицы объема и H полная удельная энтальпия:

u2 + u2 p (E + p) r z (5.19) E= +, H=, p = RT.

2 1 Здесь ur и uz проекции вектора скорости u на оси r и z. Влияние внешних сил не учитывается.

Компоненты вектора плотности потока массы jm вычисляются по формулам:

(5.20) jmr = (ur wr ), jmz = (uz wz ), где 1 p (ru2 ) + wr = (ur uz ) +, r r r z r 1 p (u2 ) + wz = (rur uz ) +.

z r r z z 5.2. Система уравнений в цилиндрической геометрии Компоненты тензора вязких напряжений выпишем в виде, удобном для программной реализации:

ur rr = N S + ur wr + R, N S = 2µ µ div u, rr rr r ur uz rz = N S + ur wz, N S = N S = µ +, rz rz zr z r zr = N S + uz wr, (5.21) zr ur = N S + R, N S = 2µ µ div u, r uz NS NS zz = zz + uz wz + R, zz = 2µ µ div u, z N S, N S, N S, N S, N S здесь rr компоненты навье-стоксовского rz zr zz тензора вязких напряжений.

Величины wr, wz, R вычисляются по формулам:

ur ur p wr = ur + uz +, r z r uz uz p wz = ur + uz +, r z z p p R = ur + uz + p div u, r z а дивергенция вектора скорости div u равна:

1 (rur ) uz div u =+.

r r z Компоненты вектора теплового потока q находятся по форму лам:

qr = qr S ur Rq, N qz = qz S uz Rq, N (5.22) ur p uz p 1 Rq = + + pur + puz, 1 r 1 z r z где навье-стоксовские слагаемые qr S и qz S задаются формулами:

N N T T qr S = N qz S = N (5.23),.

r z 108 Гл.5. Алгоритмы решения задач газовой динамики 5.3 Граничные условия Приведенные системы уравнений следует дополнить начальными и граничными условиями. Постановка этих условий определяется кон кретной решаемой задачей. В качестве примера приведем вариант постановки граничных условий для задачи о сверхзвуковом течении газа в окрестности цилиндра в осесимметричной геометрии. Схема расчетной области для этого течения приведена на рис. 5.7.

Профиль течения на входной границе, расположенной в плоско сти z = L (вертикальная граница), зададим в виде (5.24) =, uz = uz, ur = 0, p = p.

На оси симметрии, совпадающей с осью z, поставим условия сим метрии:

uz p (5.25) = 0, = 0, ur = 0, = 0.

r r r На свободных границах, где предполагается, что газ вытекает из рассматриваемой области, поставим так называемые "мягкие" гра ничные условия, или условия "сноса". В этих условиях предполага ется равенство нулю нормальных производных плотности, давления и компонент скорости. Например, если граница находится в плоско сти z = 0, то такие условия запишутся как:

uz ur p (5.26) = 0, = 0, = 0, = 0.

z z z z Пусть одной из границ является непроницаемая твердая стенка.

На твердой стенке можно поставить условия непротекания un = и прилипания ut = 0 или скольжения ut /n = 0 для скорости, где un и ut соответственно нормальная и тангенциальная компоненты скорости течения. Эти условия следует дополнить граничным усло вием для температуры, которое определяется физическими услови ями на твердой стенке.

В КГД системе уравнений вектор плотности потока массы вы числяется как QGD (5.27) jm = u w = u div(u u) + p.

5.3. Граничные условия Для того, чтобы нормальная составляющая потока массы, проте кающая через твердую границу и вычисляемая по формуле (5.27), была равна нулю, условия непротекания для скорости необходимо дополнить условием для давления вида p/n = 0. Это соотноше ние представляет собой дополнительное условие, необходимое для замыкания КГД уравнений и отличающее постановку задачи в рам ках КГД системы от постановки этой же задачи в рамках модели Навье–Стокса.

Для определенности на боковой поверхности r = Rc цилиндра (рис. 5.7) в качестве граничных условий поставим условия прили пания для скорости, зададим постоянную температуру и будем по лагать, что поверхность непроницаема. Тогда граничные условия запишутся следующим образом p (5.28) uz = 0, ur = 0, T = T0, = 0.

r Выпишем выражения для теплового потока и силы трения на границе с условием непротекания un = 0. Пусть граница располо жена в плоскости z = Lc (торец цилиндра). Сила трения на этой границе определяется компонентой тензора вязких напряжений zr, а тепловой поток qz равен:

qz = qz S uz Rq, N ur ur p zr = N S + uz uz (5.29) + ur +.

zr z r r С учетом условия непротекания на стенке (uz = 0), сразу полу чаем выражение для теплового потока и коэффициента трения на стенке в виде:

T qz = qz S = N, z ur zr = N S (5.30) =µ.

zr z То есть для КГД уравнений выражения для теплового потока и силы трения на твердой стенке (5.30) совпадают с соответствующи ми величинами для уравнений Навье–Стокса.


110 Гл.5. Алгоритмы решения задач газовой динамики 5.4 Безразмерный вид уравнений Для численного решения уравнений газовой динамики бывает удобно представить их в безразмерном виде. Это позволяет, во первых, оперировать при расчетах с величинами порядка единицы, во-вторых, позволяет выделить в уравнениях безразмерные коэф фициенты, от которых зависит решение задачи так называемые параметры подобия.

Выберем в качестве основных размерных параметров характер ный линейные размер Rc (например, радиус цилиндра), плотность и скорость звука c в набегающем невозмущенном потоке.

Запишем соотношения между размерными и безразмерными ве личинами, обозначая безразмерные величины знаком "тильда " :

p = p c2, =, u = u c, c = c c, Rc (5.31) x = x Rc, t=t, c p c2 p p 1 c c = T.

T= = = R R R R Введем числа Маха и Рейнольдса:

u u Rc (5.32) M a =, Re =.

c µ Уравнение полной энергии не изменяет своего вида:

u2 u p p (5.33) E= +, E= +.

2 ( 1) 2 ( 1) При обезразмеривании скорость звука преобразуется как (5.34) c= RT, c= T.

Уравнение состояния:

c T p c2 = T (5.35) p = RT, p=.

R 5.5. Разностная аппроксимация Таким образом, уравнения связи (5.34) и (5.35) после обезразмери вания изменили свой вид.

Подставляя соотношения (5.31) в уравнения КГД системы, мож но убедиться, что приведение к безразмерному виду не меняет вида исходных уравнений. При этом безразмерные коэффициенты вязко сти, теплопроводности и параметр вычисляются как M a (5.36) µ= T, Re 1 M a (5.37) = T.

( 1)P r Re Ma T (5.38) =.

ReSc p Далее, используя безразмерные КГД уравнения и уравнения свя зи, знак "тильда" будем опускать.

5.5 Разностная аппроксимация Рассмотрим, для определенности, аппроксимацию уравнений в ци линдрической геометрии. Обозначим символом h множество узлов (i, j) пространственной сетки. Для простоты рассмотрим равномер ную сетку, в которой обозначим шаги по пространству как hz и hr.

Все газодинамические величины (плотность, компоненты ско ростей uz и ur, давление p) будем относить к узлам сетки h. Зна чение газодинамических величин в полуцелых узлах (5.39) (zi±0.5, rj ), (zi, rj±0.5 ), (zi±0.5, rj±0.5 ) будем вычислять как среднее арифметическое их значений в прилегающих узлах. То есть значение функции из множества {, ur, uz, p} в полуцелых узлах вычисляется как (5.40) i±0.5,j = 0.5(i±1,j + i,j ), i,j±0.5 = 0.5(i,j±1 + i,j ), i±0.5,j±0.5 = 0.25(i±1,j±1 + i,j±1 + i±1,j + i,j ).

112 Гл.5. Алгоритмы решения задач газовой динамики Таким образом, для вычисления пространственных производных ис пользуется девятиточечный шаблон (рис. 5.1).

Рис. 5.1. Шаблон для аппроксимации разностных производных Для функций газодинамических параметров f = f (, uz, ur, p) положим (5.41) fi,j = f (i,j, (ur )i,j, (uz )i,j, pi,j ), fi±0.5,j = f (i±0.5,j, (ur )i±0.5,j, (uz )i±0.5,j, pi±0.5,j ), fi,j±0.5 = f (i,j±0.5, (ur )i,j±0.5, (uz )i,j±0.5, pi,j±0.5 ), и т.д.

Для численного решения системы (5.15)–(5.18) используем яв ную по времени разностную схему. Производные по времени аппрок симируются разностями вперед с первым порядком точности. Про странственные производные аппроксимируем центральными разно стями со вторым порядком точности.

Аппроксимируем дифференциальное уравнение (5.15) разност ным:

ij ij + [(rjmr )i,j+0.5 (rjmr )i,j0.5 ] + t rj hr + [(jmz )i+0.5,j (jmz )i0.5,j ] = 0, (5.42) hz 5.5. Разностная аппроксимация где t шаг по времени. Величина, помеченная верхним индексом i,j, вычисляется на следующем временном слое.

Остальные уравнения системы (5.15)–(5.18) аппроксимируются аналогично:

(ur )ij (ur )ij + t + [(rjmr )i,j+0.5 (ur )i,j+0.5 (rjmr )i,j0.5 (ur )i,j0.5 ] + rj hr 1 + [(jmz ur )i+0.5,j (jmz ur )i0.5,j ] + [pi,j+0.5 pi,j0.5 ] = hz hr = [(rrr )i,j+0.5 (rrr )i,j0.5 ] + rj hr ( )ij, (5.43) + [(zr )i+0.5,j (zr )i0.5,j ] hz rj (uz )ij (uz )ij + t + (rjmr )i,j+0.5 (uz )i,j+0.5 (rjmr )i,j0.5 (uz )i,j0.5 + rj hr 1 + (jmz uz )i+0.5,j (jmz uz )i0.5,j + pi+0.5,j pi0.5,j = hz hz = rj+0.5 (rz )i,j+0.5 rj0.5 (rz )i,j0.5 + rj hr (zz )i+0.5,j (zz )i0.5,j, (5.44) + hz Eij Eij + (rjmr )i,j+0.5 Hi,j+0.5 (rjmr )i,j0.5 Hi,j0.5 + t rj hr + (jmz H)i+0.5,j (jmz H)i0.5,j + hz 114 Гл.5. Алгоритмы решения задач газовой динамики + rj+0.5 (qr )i,j+0.5 rj0.5 (qr )i,j0.5 + rj hr + (qz )i+0.5,j (qz )i0.5,j = hz = (rrr )i,j+0.5 (ur )i,j+0.5 (rrr )i,j0.5 (ur )i,j0.5 + rj hr + rj+0.5 (rz uz )i,j+0.5 rj0.5 (rz uz )i,j0.5 + rj hr + (zr ur )i+0.5,j (zr ur )i0.5,j + hz (zi, rj ) h. (5.45) + (zz uz )i+0.5,j (zz uz )i0.5,j, hz Далее строятся разностные аппроксимации непосредственно для компонент векторов плотности потока массы jm, теплового потока q и для компонент тензора вязких напряжений, записанных в виде (5.20)–(5.23). Эти величины аппроксимируются в полуцелых точках. Например, величины, входящие в уравнение (5.42), аппрок симируются следующим образом:

(rjmr )i,j±0.5 = i,j±0.5 [rj±0.5 (ur )i,j±0.5 (rwr )i,j±0.5 ], (jmz )i±0.5,j = i±0.5,j [(uz )i±0.5,j (wz )i±0.5,j ].

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

Численное решение системы разностных уравнений (5.42)–(5.45) на каждом временном слое осуществляется в следующей последова тельности. По известным с предыдущего временного слоя газодина мическим параметрам, uz, ur и p вычисляются значения величин в полуцелых точках (5.40). Затем вычисляются значения потоков jm, q и в соответствующих полуцелых точках сетки, после чего эти значения подставляются в разностные выражения (5.42)–(5.45).

К системе разностных уравнений (5.42)–(5.45) добавляются на чальные и граничные условия.

Для единообразного вычисления газодинамических величин во всех внутренних точках расчетной области, включая приграничные 5.5. Разностная аппроксимация узлы, вводится система фиктивных ячеек, примыкающих к каждой из границ. При этом граница области располагается в полуцелых точках (см. рис. 5.8). Значения плотности, компонент скорости и давления в фиктивных ячейках задаются таким образом, чтобы ап проксимировать соответствующее граничное условие.

Например, пусть на границе, расположенной в точке i = 1/2, задано значение функции fw. Точка i = 0 является фиктивной, а точка i = 1 является ближайшей прилегающей к границе внутрен ней точкой, значение в которой вычисляется на каждом временном слое. Тогда f0 в фиктивной точке выбирается из условия f0 + f fw = f0 = 2fw f1.

Если на границе задано условие на производную вида f /n = 0, то f0 в фиктивной точке равно f0 f = 0 f0 = f1.

h Таким образом, алгоритм нахождения плотности, компонент скорости и давления на каждом временном слое состоит из двух этапов. Сначала, на основе вычисленных значений заполняются фиктивные ячейки по указанному выше правилу, затем вычисля ются значения ij, (ur )ij, (uz )ij и Eij во всех внутренних точках на следующем временном слое.

Стационарное решение находится методом установления и счи тается достигнутым при выполнении следующего критерия:

1 ij ij (5.46), Nr Nz ij t (zi,rj )h в котором невязка может варьироваться в зависимости от требова ний к точности решения. Здесь Nr и Nz число узлов сетки по r и z. В приведенных далее расчетах значение выбиралось в диапазоне от 103 до 107.

Возможны и другие критерии установления стационарного ре шения, например:

ij ij max.

t (zi,rj )h 116 Гл.5. Алгоритмы решения задач газовой динамики 5.6 Введение искусственной диссипации При численном моделировании газодинамических течений возника ют различные особенности решения, например, разрывы решения ударные волны и контактные разрывы, пограничные слои и другие зоны с большими градиентами параметров. В этих случаях оказы вается невозможным провести расчет непосредственно по разност ным схемам, полученным центрально-разностными аппроксимация ми исходных уравнений.

Известным способом расчета ударной волны без явного выделе ния на сетке ее фронта является метод "размывания"фронта за счет введения в систему разностных уравнений некоторых диссипатив ных членов, называемых псевдовязкостью или искусственной вязко стью. Такие подходы к решению задач газовой динамики детально разработаны и изложены, например, в монографиях [5,70,87, 88,91].

Вводимые добавки моделируют действие реальной вязкости, то есть приводят к диссипации кинетической энергии. Например, хо рошо известны добавки к давлению вида p p +, где величина интерпретируется как искусственная вязкость. Наи более часто рассматривается линейная вязкость u =, x или квадратичная вязкость u =.

x Здесь – некоторый коэффициент, пропорциональный шагу про странственной сетки h. При решении задач газовой динамики для политропного газа, теплопроводность которого равна нулю, исполь зуют вязкость Неймана–Рихтмайера (см., например, [87]), где u = h.

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

В алгоритмах, основанных на КГД уравнениях, в качестве регу ляризатора выступают слагаемые, пропорциональные коэффициен ту.

Для обеспечения устойчивости КГД алгоритма в формулу (5.38) для вычисления релаксационного параметра вводится добавка, пропорциональная шагу пространственной сетки. Эта добавка опре деляет дополнительную искусственную диссипацию. Безразмерный коэффициент вычисляется как:

Ma T h (5.47) = +, Re p Sc c где численный коэффициент порядка единицы, который опре деляется подбором. Для выписанной выше разностной схемы h = = h2 + h2. Коэффициенты динамической вязкости µ и теплопро z r водности вычисляются через релаксационный параметр как:

pSc (5.48) µ = pSc, =.

P r( 1) Таким образом, стабилизирующая добавка h/c включается в µ и, а следовательно, и в выражения для теплового потока и силы трения на границе (5.30).

Исходя из формулы (5.47) сеточную вязкость формально можно считать малой, когда M a T +1/ (5.49) h.

Re p 118 Гл.5. Алгоритмы решения задач газовой динамики Введение искусственной диссипации, пропорциональной шагу сетки h, делает результирующую разностную схему схемой первого порядка точности по пространству.

Таким образом, разностные схемы на основе КГД уравнений ти па (5.42)–(5.45) аппроксимируют начально-краевую задачу с пер вым порядком точности по времени и по пространству. Схема явля ется явной, однородной и консервативной. Дополнительные слагае мые, пропорциональные h/c, интерпретируются как искусственные регуляризаторы.

5.7 Задача о распаде сильного разрыва В качестве примера использования построенного выше алгоритма рассмотрим одномерную задачу о распаде сильного разрыва в при ближении невязкого нетеплопроводного газа. То есть будем решать задачу в рамках уравнений Эйлера. Эта задача является известным тестом для оценки устойчивости и точности численных алгоритмов, и различные подходы к рассмотрению этого течения приведены, в частности, в [23, 87, 98].

В этом случае параметр в (5.47) определяет величину искус ственной диссипации и вычисляется как h (5.50) =.

c КГД уравнения для одномерного плоского течения имеют вид jm (5.51) + = 0, t x (u) (jm u) p xx (5.52) + + =, t x x x E (jm H) q (xx u) (5.53) + + =.

t x x x Здесь E и H полная энергия единицы объема и полная удельная энтальпия, которые вычисляются по формулам: E = u2 /2 + p/( 1) и H = (E + p)/. Вектор плотности потока массы вычисляется 5.7. Задача о распаде сильного разрыва как jm = (u w), где (u2 + p).

w= x Компонента тензора вязких напряжений, входящая в систему уравнений (5.51)–(5.53), определяется как 4 u u p p u xx = µ + u u + + u + p.

3 x x x x x Вектор теплового потока q равен T u p q = u + pu.

x 1 x x Обезразмеривание не изменяет вида уравнений. Релаксационный параметр и коэффициенты вязкости и теплопроводности в безраз мерном виде вычисляются как h · p · Sc =, µ = · p · Sc, =.

c P r( 1) Введем равномерную сетку по координате x с шагом h, коорди натами узлов xi и значениями индекса i = 0... Nx 1, а также сетку по времени с шагом t. Значения всех газодинамических величин скорости, плотности, давления определяются в узлах сетки.

Значения потоков определяются в полуцелых узлах. Для решения задачи (5.51)–(5.53) используем явную по времени схему следующе го вида:

t i = i j jm,i1/2, h m,i+1/ t i ui = i ui + xx,i+1/2 xx,i1/2 pi+1/2 pi1/ h jm,i+1/2 ui+1/2 jm,i1/2 ui1/2, 120 Гл.5. Алгоритмы решения задач газовой динамики t Ei = Ei + xx,i+1/2 ui+1/2 xx,i1/2 ui1/ h jm,i+1/ qi+1/2 qi1/2 Ei+1/2 + pi+1/ i+1/ i u jm,i1/2 i Ei1/2 + pi1/2, pi = ( 1) Ei.

i1/2 Дискретный аналог потока массы jm имеет вид jm,i+1/2 = i+1/2 (ui+1/2 wi+1/2 ), i+1/2 (i+1 u2 + pi+1 i u2 pi ).

wi+1/2 = i+1 i i+1/2 h Дискретные выражения для xx и q выписываются аналогично.

Задача решается на отрезке 0 200. Начальные условия x представляют собой разрыв в точке x = 100. Значения газодинами ческих параметров справа и слева от разрыва составляют 8, x 99 480, x (x, 0) =, p(x, 0) =, u(x, 0) = 0.

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

Расчеты проведены для = 5/3, P r = 2/3, Sc = 1. Параметры проведенных расчетов систематизированы в табл. 5.1, где указан номер варианта расчета, величина пространственного шага h, число точек сетки Nx, шаг по времени t, число шагов по времени Nt, необходимое для достижения безразмерного времени t0 = 4. Этот момент времени соответствует целому ряду тестов, приведенных, например, в [23, 52, 59].

5.7. Задача о распаде сильного разрыва N варианта h Nx t Nt t 2 · 1 1 200 2 · 103 2 · 2 0.5 400 2 · 103 2 · 3 0.25 800 2 · 103 2 · 4 0.125 1600 2 · 104 2 · 5 0.0625 3200 2 · 104 2 · 6 0.03125 6400 2 · 104 Таблица 5.1. Варианты расчета задачи о распаде разрыва, = 0. Результаты расчетов на момент времени t0 = 4 представлены на рис. 5.2–5.6. На рис. 5.2–5.4 показаны распределения плотности, дав ления и скорости, проведенных на сгущающихся пространственных сетках. Там же сплошной линией приведено автомодельное реше ние этой задачи согласно [87]. Приведены распределения газодина мических величин во всей области расчета, а также представлены фрагменты решения в зоне ударной волны. Пунктирные линии соот ветствуют вариантам с 1 по 6, соответственно таблице. Варианту соответствует пунктирная кривая, которая больше всего сглаживает решение. Из приведенных рисунков наглядно видна сходимость чис ленного решения к автомодельному при сгущении пространственной сетки.

На рис. 5.5 5.6 показано влияние выбора коэффициента, опре деляющего величину параметра регуляризации, на точность чис ленного решения для h = 0.03125. Видно, что с уменьшением точность численного решения возрастает. Однако при слишком ма леньких значениях = 0.02 и 0.1 в области фронта ударной волны появляются осцилляции, величина которых неограниченно возрас тает при уменьшении. В ряде случаев подавление этой неустойчи вости возможно при уменьшении шага по времени t.

Для уточнения влияния коэффициента на рис. 5.6 изображен фрагмент графика плотности в зависимости от величины при = 0.5, 0.4, 0.3 и 0.2 для h = 0.03125 (h/c = 0.003125). Видна быстрая сходимость решения с уменьшением коэффициента.

При = 0.2 ударная волна размывается примерно на 4 шага сетки, что соответствует алгоритмам сквозного счета повышенного 122 Гл.5. Алгоритмы решения задач газовой динамики порядка точности для таких задач. Однако для = 0.1 уменьше ние шага по времени даже в 100 раз приводило к осциллирующему решению. Оптимальное значение параметра регуляризации, обеспе чивающее максимальный шаг по времени и наилучшую точность в этой задаче, достигается при = 0.2 0.3 (см. параграф 5.11).

0 50 100 150 x 130 135 140 145 150 155 x Рис. 5.2. Распределение плотности вдоль оси x (внизу в увеличенном масштабе) 5.7. Задача о распаде сильного разрыва p 0 50 100 150 x p 130 135 140 145 150 155 x Рис. 5.3. Распределение давления вдоль оси х (внизу в увеличенном масштабе) 124 Гл.5. Алгоритмы решения задач газовой динамики u 0 50 100 150 x u 130 135 140 145 150 155 x Рис. 5.4. Распределение скорости вдоль оси х (внизу в увеличенном масштабе) 5.7. Задача о распаде сильного разрыва 143 143.5 144 144.5 x p 143 143.5 144 144.5 x Рис. 5.5. Графики плотности (вверху) и давления (внизу) при = 0.02 – тонкая линия, 0.1 – штрих-пунктир, 0.5 – пунктир, 1 – сплошная линия 126 Гл.5. Алгоритмы решения задач газовой динамики Рис. 5.6. = 0.5 соответствует линия 1, = 0.4 – линия 2, = 0.3 – линия 3, = 0.2 – линия 4.

5.8 Задача о течении в окрестности цилиндра Постановка задачи и численный алгоритм В качестве примера реализации двумерного разностного алгоритма приведем результаты расчета задачи о течении в окрестности ци линдрического торца радиуса Rc, который помещен в однородный сверхзвуковой поток вязкого сжимаемого теплопроводного газа, па раллельный оси цилиндра. Расчетная область для этой задачи при ведена на рис. 5.7.

Введем цилиндрическую систему координат, направив ось z вдоль оси цилиндра, а ось r вдоль поверхности торца (рис. 5.7).

Размеры области расчета L H, размеры цилиндра Lc Rc, где Rc радиус цилиндра.

Рассматриваем течение одноатомного газа твердых сфер со сле дующими параметрами: = 5/3, P r = 2/3, Re = 1000, Sc = 0.77, = 0.5. В качестве начальных условий для системы (5.15)–(5.18) 5.8. Задача о течении в окрестности цилиндра Рис. 5.7. Схема расчетной области в задаче об обтекании цилиндрического торца используем параметры набегающего потока:

(5.54) |t=0 =, uz |t=0 = U, ur |t=0 = 0, p|t=0 = p.

После введения безразмерных величин, эти параметры равны = = 1, U = M a и p = 1/.

На правой входной границе (z = L, 0 r H) задаются пара метры набегающего потока (5.54):

(5.55) = 1, uz = M a, ur = 0, p=, на оси симметрии (Lc z L, r = 0) "условия симметрии":

uz p (5.56) = 0, = 0, ur = 0, = 0, r r r на торце цилиндра (z = Lc, 0 r Rc ) "условия скольжения" для скорости:

ur p (5.57) = 0, uz = 0, = 0, = 0, z z z 128 Гл.5. Алгоритмы решения задач газовой динамики на боковой поверхности цилиндра (0 z Lc, r = Rc ) "условия прилипания" для скорости:

p (5.58) = 0, uz = 0, ur = 0, = 0.

r r Поверхность цилиндра предполагается адиабатической.

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

Выходная граница (z = 0, Rc r H):

uz ur p (5.59) = 0, = 0, = 0, = 0, z z z z свободная граница (0 z L, r = H):

uz ur p (5.60) = 0, = 0, = 0, = 0.

r r r r При аппроксимации граничных условий с помощью фиктив ных узлов ячейка, соответствующая угловой точке с координатами Lc, Rc, оказывается неоднозначной (рис. 5.7). А именно, при вычис лении значения функций в фиктивном узле (i 1, j 1) в одном случае требуется использовать значение функции в узле (i 1, j), а в другом – значение в расчетном узле (i, j 1) на предыдущем временном слое (рис. 5.8).



Pages:     | 1 || 3 | 4 |   ...   | 6 |
 





 
© 2013 www.libed.ru - «Бесплатная библиотека научно-практических конференций»

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