A method of direct numerical simulation of turbulent flows of viscous heat-conducting gas in curved channels.pdf Большим достижением в теории турбулентности явились прямые расчеты турбулентных течений на основе интегрирования уравнений Навье - Стокса [1]. Однако все они проводились для течений, ограниченных прямыми стенками. Из экспериментов и теоретических расчетов известно о возникновении вихрей и турбулентности за плохо обтекаемыми телами. Прямые расчеты турбулентных течений в каналах с криволинейными границами авторам не известны. Для решения задач расчета течений в каналах с переменным контуром целесообразно применять системы координат, связанные с этим контуром. Координатная сетка должна быть достаточно гладкой, чтобы не вносить ошибок в аппроксимацию дифференциальных уравнений. Оптимальной является ситуация, при которой и границы, и координатная сетка описываются аналитическими функциями. Пусть требуется решить задачу о турбулентном течении в канале с двумя криволинейными противолежащими стенками. Будем рассматривать канал, координаты стенок которого не изменяются в направлении, перпендикулярном некоторой плоскости. В то же время сечение канала упомянутой плоскостью является криволинейным. В таком канале удобно ввести декартову координату z в направлении, перпендикулярном плоскости, и не зависящие от z ортогональные криволинейные координаты qx и q2 в сечении канала плоскостью z = const. Если при этом каждая часть границы криволинейной области совпадает с соответствующей координатной линией, то при конформном отображении этой области в прямоугольник плоскости qj, q2 в силу сохранения углов при конформном отображении получим ортогональную прямоугольную систему координат. Введенная таким образом система координат qx, q2, z будет прямоугольной декартовой системой в трехмерном пространстве. В качестве примера подходящего изображения рассмотрим отображение полуплоскости с выброшенным отрезком на полуплоскость. Если из полуплоскости > 0 исключить отрезок (a, a + h), то такое отображение дается формулой ro = V(§-a )2 +h2 +a, (1) где ю = q1 + iq2 , § = х + iy . При этом прямым линиям q2 = const в плоскости q1, q2 в плоскости х, y соответствуют линии h2 y = qiA1+- (х - a)2 + q2 а прямым q1 = const в плоскости q1, q2 - линии h2 • = a + (q1 - a2 1 -v (q1- a )2 + y2 Эти линии в плоскости q1, q2 образуют прямоугольную сетку, а в плоскости х, y - криволинейную ортогональную систему координат (рис. 1). При этом линии q2 = 0 соответствует линия y = 0 в плоскости х, y . (§-a 2 т-т d ю Производная — = - 0.3 -0.2 -0.1 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 х Рис. 1. Сетка координат на плоскости х, y обращается в нуль при § = a и в бесконечd§ V(§-a 22 +h2 ность при § = a + ih . Если рассматривать любую часть области при q2 > 0 , то отображение везде будет аналитическим. Приняв за границы канала в области q1, q2 две прямых линии q11 = const > 0 и q12 = const > q11, в области х, y получим две криволинейные границы, между которыми отображение будет аналитической функцией. Если взять q12 >> 0, то верхняя граница в плоскости х, y будет близка к прямой линии. Меняя высоту исключенного отрезка h и величину q11, можно изменять высоту и кривизну границы в области максимума, исследовать поведение решений уравнений Навье - Стокса в зависимости от параметров препятствия. Однако главное достоинство введенного отображения заключается в преобразовании криволинейного канала в прямоугольный канал в координатах q1, q2, z. Отметим, что рассматриваемое отображение описывает течение идеальной несжимаемой жидкости с линиями тока q2 = const , которые сгущаются в области препятствия. Последнее обстоятельство является немаловажным для повышения точности численного решения уравнений Навье - Стокса. Аналогично конформное отображение можно использовать для канала прямоугольного сечения, который получается из ограничения рассмотренной выше области двумя плоскостями z1 = const и z2 = const. Полученная область в результате описанного выше конформного отображения также будет прямоугольной областью в трехмерном пространстве. Уравнения Навье - Стокса в криволинейных ортогональных координатах широко используются при решении задач гидродинамики. В наиболее общем виде они приведены в [1], где для случая вязкости, не зависящей от температуры, они записаны в дивергентной форме. В этой записи для учета преобразования координат применяются коэффициенты Ляме [2]. + f— dx dy dy_ dq2 dz dq2. cX 5q2 H1 = H 2 = + Г dz l5q3 dq3 J f5q3 Запишем систему уравнений Навье - Стокса в криволинейной системе координат ql, q2, z, в которой координаты q1 и q2 задаются преобразованием (1). Выделяя в (1) действительную и мнимую части и решая полученную систему уравнений относительно x и y , найдем H3 = (2) £ {[( - a )2 - q22 - h2 ]+>/[(Г- a )2 - q22 - h2 ]2 + 4 ( - a) 2 2 q2 x = a + (q - a )q2 (3) y = 2{[( - a)2 - q2 - h2 ] + ^[(q - a)2 - q22 - h2 ]2 + 4(q - a)2 q22 | Уравнения (3) дают нам искомые связи x = f (q1, q2), y = f (q1, q2). С их помощью по формулам (2) можно вычислить коэффициенты Ляме. Однако последние легче получить, если заметить, что для конформного преобразования 1 = f и ^l2 =Г ^12+Г 2=f jl Uq2 дю М1J Wq J \dq2 а z = q3 не зависит от q1, q2. Поэтому dL H3 = 1. (4) Замечая, что равенства (4) не зависят от конкретного вида конформного отображения, в дальнейшем введем обозначения H1 = H2 = H , H3 = 1. (5) H1 = H 2 = d ю Для отображения (1) = . ю a . Откуда получаем d(D ^(ю-a 22-h2 ю-a ю-a 2 d § d ю ^(ю-a22 -h2 V(ra-a22-h йдем (q1 - a 22 + q2 Раскрывая обозначения ю и ю , найдем 2 / \2 . „2 (6) d § d ю q1- a22 + q22 2] - q [( - a22 + ql 2] •h2 + h4 Запишем уравнения Навье - Стокса в системе координат q1, q2, z. Запись произведем, учитывая коэффициенты Ляме (5) и не подставляя для краткости конкретных выражений для H . Такой подход позволяет получить уравнения, справедливые для произвольного конформного отображения § = f (ю2. Уравнение неразрывности имеет вид 1HpV>1HpV»2HvH2 (7) где p - плотность жидкости; Vq , Vq , Vq - проекции скорости жидкости на оси q1, q2, q3 соответственно, q3 = z . Уравнения движения запишутся в форме |h ^ +a h K h pV,V, 2 + A. h ^ = = f -p".". f-+H И -H Ц -f H ± .div a+2H 2 (Div ц.ф2„ , W iH 'p'* +£ H p"v +i|rH < H 2f = P>i f-^'',2 + Н^,2 -H f-"J H ^ " + 2H2 (Div ^Ф2,2, (9) H p"*"* H p'«2 '«3 HJ p'2 = dp 2 я = H2pFq -H2---H2-цdivQ + 2H2 (Divц-Ф2 . (10) 3 dq3 3 dq3 q3 В этих уравнениях использованы обозначения: P - давление; F^, Fq , Fq^ -проекции массовых сил на оси q1, q2, q3 соответственно; ц - вязкость жидкости; div о = -L fA HVn + —HVn + —H V H2 ydq1 q1 dq2 q2 dq3 q3 входящие в уравнения (8) - (10) проекции вектора Div ц • Ф имеют вид = fA H 2цФ, dq1 qq +dq2 H ЦФ 9291 dq: =J_ 91 " H' 1 d + d H цФ dq1 9192 dq2 +A H 2цФ dq3 д h 9292 q3q2 d H (DiV = HpT dH 1 н2 d92 Цфq1q2 dH Цфq,q, dH dH ~ q2q2 q1q1 d 2 H цФ99 +-H цФ99 +-H 2цФ [dq1 q1q3 dq2 q2q3 dq; H2 dq1 H2 dq2 q3q3 где компоненты тензора скоростей деформации Ф Ф 9393 H dq2 H dq1 dq3 Ф + ^2._ q1q1 H dq1 H2 dq2 1 dFq F dH , 1 dF Fq, dH +_q1_ q2 q2 U p,„ tt2 1 dFq,. + .dFq1 1 dF* + _ 1 dF* ^ dH Fq21 dH ф = 1 9192 2 Ф =q1q3 2 (11) H dq1 H dq2 H2 dq2 H2 dq1 H dq1 dq3 Ф =ф ф =ф ф =ф q3q1 9193' q3q2 q2 93' q2q1 q\q2 1 dFq dFq q3. + . q2 Ф q2q3 H dq2 dq3 В заключение запишем уравнение энергии —pH2E+—HpF E+—HpF E+—H2pF E= dt dq1 1 dq2 2 dq3 3 =-—HFq P -—HFq P -—H 2Fq P+—2цЯ (Ф qqFq +Ф9 q Fq +Ф9 q Fq )+ dq1 q1 dq2 q2 dq3 q3 dq3 V 9191 q1 q2q1 q2 q3q1 93 7 +—2цЯ(qq Fq +Ф9 q Fq +Ф9 q Fq )+—2цИ2 (qq Fq +Ф9 q Fq +Ф9 q Fq )dq V q1q2 q1 q2q2 q2 q3q2 q3 / dq ^ q1q3 q1 q2q3 q2 q3q3 q3 / -H^Fqi divQ-—— H^Fq2 dr^——H2^Fq3 divQ+ 3dq2 3d93 1 ^ d -2 эфф ^ dq1 эфф dq1 dq2 3dq1 ' q1 dT d d „ +—X, dq2 dqH21'эфф d|+pH2 ( РЪ +F* ^ +Fq3 ^ ). (12) Наряду с обозначениями, использованными выше, в (12) применены обозначе ния: E = евн + |Q| j2 - полная энергия единицы массы газа; 1эфф - эффективный коэффициент теплопроводности, записанный с учетом излучения; Т - температура газа; проекции тензора напряжений на оси криволинейных координат: ^ =(2ЦФ^ 2ЦФ929:,2^Ф9391 ) , Fq2 = (2ЦФ 9192 , 2ЦФ929!,2ЦФ9392 ) , F93 = ( 2цф 9193 , 2цф92 93 ,2цф9393 ) . Система уравнений (7) - (12) замыкается выражением для внутренней энергии идеального газа евн = cvT , (13) и уравнением состояния идеального газа P = pRT . (14) В качестве начальных условий для решения системы (7) - (14) в начальный момент времени t = 0 задаются скорости Vq1 (0 q^ q^ q3 2 = Vq1,0(ql, q^ q32, v«2 (0, ql, ql, q3 2 = Vq2,0(ql, ql, q32, (15) Vq3 (0, q^ ql, q3 2 = Vq3,0 (ql, ql, q32 > а также давление P (0, q1, q2, q3 2 = P0 (q1, q2, q3 2 и температура T (0, q1, q2, q3 2 = T0 (q1, q2, q3 2. Рассматриваются два случая. В первом случае исследуется течение между верхней и нижней криволинейными поверхностями, неограниченными в направлении z , в качестве граничных условий на нижней и верхней границах используются условия прилипания: = 0, у \ = 0, v \ = 0, v\ = 0, v \ = 0, v \ = 0. (16) V ql Г q2 ir, q3 ir2 ' «О it ' q iT ' «i |r q iT ' «I iT v 7 q2 q3 2 Граничные условия для температуры на этих границах задается из условия адиабатичности: = 0, ^ Г ^ dT =0. (17) г 2 На боковых границах задаются условия симметрии [1]: v\ = v\ , v\ =v\ ,v\ =v\ , q1 lz=0 q1 lz=h q2 lz=0 q2 Iz =h q3 lz=0 q3 lz=h Я=0 = Я=h , plz=0 = plz=h . В другом случае рассматривается течение в канале прямоугольного сечения при наличии стенок, перпендикулярных оси z . На этих стенках используются условия прилипания: Vq| = Vq| = 0, v\ = v\ = 0, vl = vl = 0, (19) q1 lz=0 q1 lz=h q2 lz=0 q2 lz=h q3 lz=0 q3 lz=h v ' и условия адиабатичности для температуры: (18) dT 5q3 dT 5q3 =0. =0 z=h z=0 На входе в канал задаются компоненты вектора скорости и температура. Полагаем V\ = 0, V\ = 0 и согласно рекомендации [1] q2 q1 =0 q3 q1 =0 Vq1 |qi=0 = f (q 2, q3 2 + C ( - P)fl (q 2, q3 2, \ =0 = f, (q q, 2 , (20) где f (q 2, q3 2 - заданное распределение скорости Vql во входном сечении, f3 (q 2, q3 2 - заданное распределение температуры. Слагаемое C (P0 - P 2 в выражении для Vql | =0 добавляется согласно [3] с целью вывода за пределы области решения возмущений, достигших левой границы. с, - константа порядка единицы, определяемая равенством C = (21) 1 к • M где к - показатель адиабаты, M - число Маха на входе. Выражение (P0 - P) есть разность между текущим значением давления как функции координат q 2, q3 на входе в канал и средним значением этого давления при q1 = 0. Поскольку Ц P0dq2dq3 = P, в среднем P0 - P равно нулю. Однако в каждой точке (q 2, q3) разность P0 - P отлична от нуля и способствует выводу возмущений за пределы области решения [1]. На выходе из канала в дозвуковом потоке задаем одно условие PL = ^ар + C2 (h ср - hL ) , где PL - давление на выходе из канала, P - наружное давление, даср - средний в момент времени t расход на длине канала, hL = Ц pFq dq2 dq3 - расход на вы- 5 ( L) ходе канала. Эти расходы соответствуют одному и тому же моменту времени t. Для расчета коэффициента C2 используется выражение [3] к • ML c = ■ • l C2 _ 2D + hi cp • Ml в котором D - половина расстояния в направлении координаты q2 на выходе из канала, ML - число Маха на выходе. Все остальные параметры на выходе из канала должны определяться путем экстраполяции значений этих параметров во внутренней части области решения, примыкающей к выходной границе канала [1]. При численном решении сформулированной выше задачи использовались подходы, предложенные в [1] для интегрирования уравнений Навье - Стокса в случае турбулентных течений. В основе этих подходов лежат высокоточные аппроксимации производных от параметров потока по пространственным переменным. При этом для интегрирования по времени применяются явные схемы вида /я \п 1 ( я2 Лп , (д3 Лп фИ+! =фп + Г3ф) +1 ] (д,)2 +1 Л (At)3 +..., (22) l5t J 2! ^ dt2 )К J 3ltst3 J1 ' (дк фЛ v wi / в которых в правые части подставляются точные выражения для производных вычисленные на n-м слое по времени с помощью дифференциальных v dtk / уравнений Навье - Стокса. В настоящей работе производные от параметров потока по пространству вычислялись с 8-м порядком аппроксимации. Производные по времени от параметров потока вычислялись по соотношению (22) с точностью до второго порядка. При этом вместо использования чрезвычайно громоздких аналитических выражеприменялась центральная разность f£V dt2 нии для вычисления аг 2 At что позволило существенно упростить программу ЭВМ-расчетов и сократить время вычислений. Приведем два примера расчетов на основе применения разработанной методики. В первом примере рассматривается течение между двумя поверхностями, координаты точек которых не зависят от направления z . На нижней поверхности была расположена перегородка (турбулизатор), высота которой составляла 0,1b, где b - расстояние между поверхностями. Число Рейнольдса описываемого течения было выбрано равным 50. Число Маха M = 0,5, число Прандтля 0,7 . На рис. 2 показано распределение давления в плоскости (x, у), перпендикулярной оси z. На рисунке видно, что за исключением резкого падения за преградой Рис. 2 у- 0.4 -0.3 -0.2 -0.1 -1. 8 1.9 2 2.1 2.2 2.3 2.4 2.5 2.6 2.7 x Рис. 3 изменение давления в канале носит плавный характер. На рис. 3 показаны векторы скорости в той же плоскости. Видно, что за преградой возникает вихрь. Этот вихрь является стационарным. Так как, согласно расчетам, скорости не имеют проекций на ось z , то этот вихрь представляет собой вихревую трубку с параметрами, не зависящими от z . В другом примере рассчитывалось течение при числе Re = 104 в канале прямоугольного сечения. Высота турбулизатора составляла 0,2b, где b - расстояние между стенками канала. На рис. 4 приведена картина мгновенного распределения давления, где наблюдаются нерегулярные колебания его величины, увеличивающиеся за турбулизатором. На рис. 5 показаны мгновенные значения проекций векторов скорости на плоскость (x,у), а на рис. 6, а - проекции этих векторов на плоскость, параллельную турбулизатору и расположенную на расстоянии 0,5D позади него. Видно, что течение за турбулизатором носит нерегулярный вихревой характер. На рис. 6, б приведены проекции скорости на плоскость сечения канала, находящуюся на расстоянии D за турбулизатором. Видно, что уносимые потоком вихри распадаются. На рисунке наблюдаются вихревая пелена у дна канала и два небольших вихря у плоскости симметрии. Этот факт подтверждает гипотезу Колмогорова о распаде крупных вихрей на более мелкие [4]. 1.6 1.8 2 2.2 2.4 2.6 2.8 3 3.2 3.4 x Рис. 5 У 1.4 1.2 1 б 0.2 0.4 0.6 0.8 1.2 1.4 Рис. 6 1 z У 1.4 1.2 1 0.8 0.6 0.4 0.2 а ■f t *. 0.2 0.4 0.6 0.8 1 1.2 1.4 z * '< ^ «'У * 0.8 0.6 0.4 0.2 Приведенные результаты рассчитывались на разностной сетке размером 64 х 64 х 160. Такая сетка не позволяет рассчитывать еще более мелкие вихри, которые по теории Колмогорова обеспечивают диссипацию турбулентной энергии. Для того чтобы при Re = 104 проводить расчеты таких вихрей, необходима разностная сетка, содержащая порядка 109 ячеек. Сетки таких размеров можно создавать на кластере Томского государственного университета. Поэтому в качестве дальнейшего развития созданной методики в настоящее время производится ее распараллеливание на MPI. Таким образом, разработан подход к прямому численному моделированию турбулентного течения вязкого теплопроводного газа в криволинейных каналах. Представлены результаты прямого численного моделирования турбулентного течения вязкого теплопроводного газа в канале с препятствием в виде поперечной пластины заданной высоты.

Скачать электронную версию публикации