ЖУРНАЛ ВЫЧИСЛИТЕЛЬНОЙ МАТЕМАТИКИ И МАТЕМАТИЧЕСКОЙ ФИЗИКИ, 2015, том 55, № 5, с. 886-894
УДК 519.634
МЕТОД ПРЯМОГО ЧИСЛЕННОГО МОДЕЛИРОВАНИЯ ТУРБУЛЕНТНОГО ТЕЧЕНИЯ ГАЗА В КРИВОЛИНЕЙНЫХ КООРДИНАТАХ
© 2015 г. И. М. Васенин, А. Ю. Крайнов, А. М. Липанов, Э. Р. Шрагер
(634050 Томск, пр-т Ленина, 36, ТГУ) e-mail: akrainov@ftf.tsu.ru Поступила в редакцию 18.03.2014 г.
Представлен подход к прямому численному моделированию турбулентного течения вязкого теплопроводного газа в криволинейных каналах. Приведены результаты прямого численного моделирования турбулентного течения вязкого теплопроводного газа в канале с препятствием в виде поперечной пластины заданной высоты. Библ. 4. Фиг. 6.
Ключевые слова: вязкий, теплопроводный газ, уравнения Навье—Стокса, турбулентность, прямое численное моделирование.
DOI: 10.7868/S004446691505018X
Большим достижением в теории турбулентности явились прямые расчеты турбулентных течений на основе интегрирования уравнений Навье—Стокса (см. [1]). Однако все они проводились для областей с прямолинейными стенками. Прямые расчеты турбулентных течений в каналах с криволинейными границами не известны. Для решения задач расчета течений в каналах с переменным контуром целесообразно применять системы координат, связанные с этим контуром. Координатная сетка должна быть достаточно гладкой, чтобы не вносить ошибок в аппроксимацию дифференциальных уравнений. Оптимальной является ситуация, при которой и границы, и координатная сетка описываются аналитическими функциями.
Пусть требуется решить задачу о турбулентном течении в канале с двумя криволинейными противолежащими стенками. Будем рассматривать канал, координаты стенок которого не изменяются в направлении, перпендикулярном некоторой плоскости. В то же время сечение канала упомянутой плоскостью является криволинейным. В таком канале введем декартову координату г в направлении, перпендикулярном плоскости, и не зависящие от г ортогональные криволинейные координаты q1 и q2 в плоскости г = const. И пусть при этом каждая часть границы криволинейной области совпадает с соответствующей координатной линией. Тогда при конформном отображении этой области в прямоугольник плоскости q1, q2 в силу сохранения углов при конформном отображении получим ортогональную прямоугольную систему координат. Введенная таким образом система координат q1, q2, г будет прямоугольной декартовой системой в трехмерном пространстве.
В качестве примера подходящего изображения рассмотрим отображение полуплоскости с выброшенным отрезком на полуплоскость. Если из полуплоскости Im^ > 0 исключить отрезок (a, a + h), то такое отображение дается формулой
где ю = q1 + iq2, 2 = * + iy.
При этом прямым линиям q2 = const в плоскости q1, q2 в плоскости x, y соответствуют линии
(1)
Фиг. 1. Сетка координат на плоскости x, y.
а прямым q1 = const в плоскости q1, q2 — линии
x = а + $ - а) 1 —-^-J.
V (#1 - а) + y
Эти линии в плоскости q1, q2 образуют прямоугольную сетку, а в плоскости x, y — криволинейную ортогональную систему координат (фиг. 1). При этом линии q2 = 0 соответствует линия y = 0 в плоскости x, y.
Производная
d&= (%- а) d^ fa- а)2 + к2
обращается в нуль при ^ = а и в бесконечность при ^ = а + ik. Если рассматривать любую часть области при q2 > 0, то отображение везде будет аналитическим.
Приняв за границы канала в области q1, q2 две прямые линии q11 = const > 0 и q12 = const > 0, в области x, y получим две криволинейные границы, между которыми отображение будет аналитической функцией. Если взять q12 > 0, то верхняя граница в плоскости x, y будет близка к прямой линии. Меняя высоту исключенного отрезка h и величину q11, мы можем изменять высоту и кривизну границы в области максимума, исследовать поведение решений уравнений Навье— Стокса в зависимости от параметров препятствия. Однако главное достоинство введенного отображения заключается в преобразовании криволинейного канала в прямоугольный канал в координатах q1, q2, г. Отметим, что рассматриваемое отображение описывает течение идеальной несжимаемой жидкости с линиями тока q2 = const, которые сгущаются в области препятствия. Последнее обстоятельство является немаловажным для повышения точности численного решения уравнений Навье—Стокса.
Аналогично конформное отображение можно использовать для канала прямоугольного сечения, который получается из ограничения рассмотренной выше области двумя плоскостями z1 = = const и z2 = const. Полученная область в результате описанного выше конформного отображения также будет прямоугольной областью в трехмерном пространстве.
Уравнения Навье—Стокса в криволинейных ортогональных координатах широко используются при решении задач гидродинамики. В наиболее общем виде они приведены в [1], где для случая вязкости, не зависящей от температуры, они записаны в дивергентной форме. В этой записи для учета преобразования координат применяются коэффициенты Ляме (см. [2]):
Запишем систему уравнений Навье—Стокса в криволинейной системе координат q1, z, в которой координаты q1 и q2 задаются преобразованием (1). Выделяя в (1) действительную и мнимую части и решая полученную систему уравнений относительно х и у, найдем:
х = а +
У =
- а)2 - д22 - к2] + - а)2 - я22 - к2] + 4( - а)2 _С - а) с!2_
(3)
{ - а)2 - С22 - й2] + ^[(с - а)2 - С22 - к2] + 4( - а)2 д2^
Уравнения (3) дают нам искомые связи х = /1 (с1,с2), у = f (с1,с2). С их помощью по формулам (2) мы можем вычислить коэффициенты Ляме. Однако последние легче получить, если заметить, что для конформного преобразования 2, = / (ю)
дх
+
,2 / \ 2 ду | _\ дх
кдси Удд1 а z = q3 не зависит от q1, q2. Поэтому
Н1 - н2 -
дС2
+
ду 12 _(дС2
дС2
дю
Н3 = 1.
(4)
Замечая, что равенства (3) не зависят от конкретного вида конформного отображения, в дальнейшем введем обозначения
Для отображения (1) имеем
Отсюда получаем
Н1 — Н2 — Н, Н3 — 1.
йк = ю- а йю л!(ю- а)2 - к2
(5)
ш-а
ш-а
С учетом формулы (2) имеем
йи
д/(ш-а)2 -к2 а)2 - к (1 - а)2 + С22_
^[((С1 - а)2 + С22)]2 - С[( - а)2 + ,22)]2 к2 + к4
(6)
Запишем уравнения Навье—Стокса в системе координат q1, q2, z. Запись произведем с учетом коэффициентов Ляме (4)—(6) и не подставляя для краткости конкретных выражений для Н. Такой подход позволяет получить уравнения, справедливые для произвольного конформного отображения 2, = f (ю).
Уравнение неразрывности имеет вид
|рН2 + ^р¥шН + ¥9Н + ¥9} Н2 = 0,
дг дд1 дд2 дд3
(7)
где р — плотность жидкости; Уд,Уд,Уд — проекции скорости жидкости на оси q1, q2, q3 соответственно, q3 = z.
Уравнения движения запишутся в форме
дН2рК, + — НрV? НрУУ2 + Н2рУУ = рV2дН-рУУ — +
дг н С1 5,1 С1 5с2 н С1 с 2 5с3 н С1 С3 н С2 дс1 С1 С2 5С2
+ Н2р^С| - Н— - 2Н—цШУ а + 2Н2 (У цф) ,
5с1 3 дс С1
2
2
ди 2. dt
д
д
т-2 , д
г— дн
дН
—HpVa +HpVja +HpV: + — HpVaya =pv; —-pVq,Vq—+
q2 3q, s q2 3q2 q2 dq3 F q3 q2 F s3q2 F s q— 3q2
+ H2pFq - H — - - H—| div a + 2H2 (Div |Ф) ,
2 dq2 3 dq2 q2
дH2pVa HpVq,Vq HpVq Vq H2pVq2 = H2pFq - H22 -
dt н q3 3q, qi q3 dq2 q2 q3 dq3 F S2 F S3 dq3
- 2 H2-^u div a + 2H2 (Div цф) .
3 8q3 У >q3
(9)
(10)
В этих уравнениях использованы обозначения: Р — давление; Еч,, Рд2, Е% — проекции массовых сил на оси q1, q2, q3 соответственно; ц — вязкость жидкости;
div a =
1 I д
д
д тГ2т
„ HVq, + — HVq + — HV I,
H [dq, qi dq2 q2 dq3 q3 1
входящие в уравнения (8)—(10) проекции вектора Div^ имеют вид
(Div =Тг
H2
(Div =i
H2
А
dq,
"_д_
dq,
dqi д
dq2 dq3
^ h \
dq3 д tt2
dq2 дqз
дН дН
H2 5q2 Н dq,'
^^^ q,q2 дН дН
Н2 дг, Н дq2
(Div И®),3 = р
а н^ф-3нн 2^ф»
где компоненты тензора скоростей деформации
ф = LdVaL + Vs^dH
qs Н dq, Н2 dq2
ф = idVqi + Vl dH
q2q2 H dq2 H2 dq,'
dVq
ф =-q3
S3S3
5q;
ф =
qs2
, dVq_ , 5Vq, Vq, 5H Vq„ 5H
H dq, H 3q2 H2 dq2 H2 dq, _
ф =
SS3
1 5V 5V
ф =
sr
, 5V 5V
, ^r S3 , 17' qi
ф = ф
q3q, q
\_H dq2 dq3 В заключение запишем уравнение энергии в виде
О Г,-2
ф = ф
S3S2 q2q3'
Н dq, dq3 _
ф = ф
q2q! qs2
—рН 2E + — Hp VE + — HpVqET —
df dq, ql 3q2 F q2 dq3
д
д
H pVq3 E = -f- HVqP -
dq,
(11)
- dtH lVqp+dt2цН (ф qlqlvql+фq2qlVq2+фqзqlVqз)
+dr 2^н (ф qlq2Vql + ф q2q2Vq2 + ф q3q2Vq3 )+dT 2цН 2 (ф ^^q3Vq, + ф q2q3Vq2 + ф q3q3Vq3) - (12)
dq2 dq3
-H^V divп - H^Vq divn - H2\iVq divП +
3 dq, s 3 dq2 ^ q2 3 dq3 * q3
+dr ^эфф dT+дгx ^ F+дтн ЭФФ dT+pH 2 (++VF)).
dq, dq, dq2 dq2 dq3 dq3
Здесь E = евн + 2 ¡2 — полная энергия единицы массы газа; X эф — эффективный коэффициент теплопроводности, записанный с учетом излучения; T — температура газа; проекции тензора напряжений на оси криволинейных координат следующие:
= (2иф s3S, ), FS2 = (2 Иф S^ 2 Иф s2sl,2иф S3S2 FS3 = (2 Иф SS^ Иф s2sз, 2 Иф ?3?3).
Система уравнений (7)—(12) замыкается выражением для внутренней энергии идеального газа
^вн _ CvT
и уравнением состояния идеального газа
Р = рЯГ.
В качестве начальных условий для решения системы (7)—(10), (12) в начальный момент времени I = 0 задаются скорости
Уъ (0, С1, С2, Сз) = УС1,0 (,1, С2, Сз), У?2 (0, С1, С2, Сз) = УС2,0 (,1, С2, Сз), У,3 (0, С1, С2, Сз) = УС3,0 (,1, С2, Сз),
а также давление Р (0,С1,С2,Сз) = Р0 (д1, д2, сз) и температура Т (0,С1, с2,сз) = Т0 (д1,д2,дз).
Рассматриваются два случая. В первом случае исследуется течение между верхней и нижней криволинейными поверхностями, неограниченными в направлении z, в качестве граничных условий на нижней и верхней границах используются условия прилипания
у\ = 0, у\ = 0, у\ = 0, у\ = 0, у\ = 0, у\ = 0.
Г! ' С2\Г! Сз\Г1 *?1|г2 ' Г С2\г 2 ' *?з1г2
Граничные условия для температуры на этих границах задаются из условия отсутствия тепловых потоков:
дТ дд
= 0, Т
= 0.
г2
дС2
На боковых границах задаются условия симметрии (см. [1]):
Уй1г=0 = УгJг=к, Ус\г=0 = УгJг=к, У?з1г=0 = У?з1 г=^ pz=0 = рг=к, рг=0 = рг=к' В другом случае рассматривается течение в канале прямоугольного сечения при наличии стенок, перпендикулярных оси z. На этих стенках используются условия прилипания
у\ = у\ = 0, у| = у\ = 0, У I = у\ = 0
*?1| г=0 г щ г=к ' ГЯ2\г=0 г С2\г=к ' г Сз\г=0 г Сз\г=к
и условия для температуры
дТ дд
= 0, дТ
Для дальнейшего прочтения статьи необходимо приобрести полный текст. Статьи высылаются в формате PDF на указанную при оплате почту. Время доставки составляет менее 10 минут. Стоимость одной статьи — 150 рублей.