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

Том XXXIX

УЧЕНЫЕ ЗАПИСКИ ЦАГИ 2 0 08

№4

УДК 532.529

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

А. Г. ЗДОР

Приведены физико-математическая модель и результаты численных исследований эволюции капель водного конденсата при обтекании плоской пластины сверхзвуковым двухфазным полидисперсным потоком. В учет приняты скоростная и температурная неравновес-ность фаз, а также межфракционное взаимодействие частиц. Из-за малости массовой доли конденсата в набегающем потоке его обратное влияние на газ не рассматривалось. Течение газа описывалось уравнениями Рейнольдса с использованием двухпараметрической модели турбулентности. Для разных режимов обтекания изучены характерные особенности полей основных параметров дисперсной фазы: числовых концентраций, радиусов, температур и ос-редненных скоростей микрокапель. Исследована зависимость динамики дисперсной фазы от числа Маха набегающего потока, показано влияние процессов коагуляции/дробления капель на трансформацию их массового спектра.

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

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

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

Газодинамический «фон». Обтекание пластины вязким теплопроводным совершенным газом (при у» 1.4) рассмотрено в условиях, соответствующих полету в атмосфере Земли на высоте

приблизительно 10 км под углами атаки 10 и 15° для чисел Маха набегающего потока 2 и 4. Для описания турбулентности применялась двухпараметрическая q- со модель [5]. Численное реше-

ние полных нестационарных уравнений Рейнольдса проводилось методом «крупных частиц» [6]. В зонах косых скачков уплотнения и вееров волн разрежения результаты расчетов находятся в хорошем качественном соответствии с известными классическими решениями задач обтекания тонких тел идеальным совершенным газом [7].

Дисперсная фаза. Для описания динамики конденсата использовалась модель полидис-персного аэрозоля [2—4], учитывающая упомянутые выше процессы: дробление и коагуляцию частиц, влияние турбулентных пульсаций на частоту столкновений капель, обмен импульсом и тепломассообмен (как межфракционный, так и с несущим газом). В рамках гидродинамической модели двухфазная система представляется многокомпонентной сплошной средой (при расчетах рассматривались три фракции). Считается, что все частицы данной фракции являются сферическими и в фиксированной точке пространства-времени имеют одинаковые скорости, размеры и температуры, причем для разных фракций эти величины могут отличаться. Таким образом, каждая фракция в данной точке пространства и в данный момент времени может быть охарактеризована вектором состояния ^ [ns, , \s, Ts, as). Здесь ns, ps — числовая и массовая концентрация частиц 5-й фракции; v5 — вектор их осредненной скорости; Ts, as — среднеобъемная температура капель фракции и их радиус.

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

капли обусловлены диффузией и теплопроводностью, полные тепловой поток qc и скорость из-

• с

менения массы капли т в предположении степенной зависимости динамическои вязкости газа

от температуры \im (r) = constr“ подробно описаны в работе [1]. В выражениях для тс и qc

используются параметры газа, соответствующие условиям адиабатического восстановления. Коэффициент восстановления при расчетах принимался равным единице по нескольким причинам. Во-первых, он пропорционален квадратному корню из числа Прандтля (Рг = 0.74), которое само

близко к единице. Затем, как показали вычисления, числа Маха обтекания капель всех фракций не превосходят значения, приближенно равного 0.8. Поэтому вносимые сделанным предположением ошибки являются достаточно малыми и находятся в рамках погрешностей расчета. Аналогия процессов тепло- и массопереноса в несжимаемых потоках дает основания распространить ее на случай чисел Маха, близких к единице. Этим объясняется структура использованных выра-

, Г'

жений для т при произвольных числах Маха [ 1 ].

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

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

F = |rc«2pCD|u-v|(u-v).

В континуальном случае использована трехчленная аппроксимация стандартной кривой коэффициента сопротивления с поправочным множителем, описывающим его зависимость от числа Маха (Хернер, 1958 г.). При наличии разности температур капель и несущего газа также требуется вводить поправки к коэффициенту сопротивления. Использованное при расчетах выраже-

ние для CD подробно описано в работах [1, 8]. Коэффициент сопротивления частиц конденсата

П

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

Большой интерес представляет интерполяция выражений для потоков массы, тепла и коэффициента аэродинамического сопротивления, полученных для предельных случаев континуального и свободномолекулярного режимов обтекания, на случай промежуточных значений числа Кнудсена. Даже при малом количестве фракций конденсата вычислительные затраты становятся очень существенными (особенно при использовании метода установления для плоских и пространственных задач). Поэтому представляется целесообразным отдать предпочтение простым вариантам интерполяции. Этим требованиям простоты и адекватности удовлетворяют, например, формулы типа «параллельных проводимостей», проверенные при дозвуковых режимах обтекания капель [1]:

ш^шС

'р = ^Г^с- Ч'Н™.9.С’0}-

Дробление капель под действием аэродинамических сил. Существует сложная классификация режимов дробления капель в потоках [4]. При расчетах в качестве критерия использовалось число Вебера: Ше = 2ра|у -и|2/о№ (о№. —коэффициент поверхностного натяжения воды).

Дробление моделировалось по схеме удвоения: если локальное значение числа Вебера данной фракции превосходило критическое, считалось, что капля распадается на две одинаковые части. Затем находились новые величины Ч^е и снова сравнивались с WeKp. В случае We> Weкp =16

процесс повторялся.

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

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

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

При моделировании процессов коагуляции важную роль играют | — средние относи-

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

21 I

К^п^п5, где =%{а]+ I I — ядро интеграла столкновений. Величины Дзависят как от осредненных скоростей фракций Yj,vs) так и от их стохастических составляющих

АхТ™Л, которые обусловлены вовлечением капель в турбулентное движение несущего газа и их участием в броуновском движении (играет важную роль для мелких частиц) [2]:

ДуВг

При континуальном обтекании частиц на основании данных работы [2] принимается: = 2kвTfЗnajas\lm (Т) (£в — константа Больцмана). Исследование взаимодействия капель

с турбулентными пульсациями потока потребовало отдельного рассмотрения. Это взаимодействие является одной из причин взаимных столкновений частиц, а также вызывает их диффузионное движение. В работе использовались результаты исследования параметров относительного движения капель на основании статистического подхода [10—12] и следующие из них простые аппроксимирующие зависимости, удобные для применений в численных алгоритмах. Трехмерный спектр турбулентности моделировался путем сращивания аппроксимаций Кармана и Пао [13]. С его учетом получены строгие результаты для функций, описывающих вовлечение частиц в пульсационное движение. Для функции вовлечения частиц справедливо [11, 12]:

ШГ

где к% — средняя удельная энергия пульсационного движения частиц 5-й фракции; £Шг ’ єіиг — средняя удельная энергия пульсационного движения несущей среды и скорость ее диссипации;

= — параметр инерционности (х5 — время динамической релаксации капель фракции

по скорости, ТЕ ~ кш /є,иг — временной интегральный масштаб турбулентности, характеризующий время жизни энергоемких вихрей); уя = |у5 - и |/м1иг — параметр скольжения (мШг — характерная скорость турбулентных пульсаций газа). Для среднего значения модуля относительной скорости двух частиц использована зависимость [11, 12]:

Ду1ш

16 к

Зп

-1/2

Коэффициенты диффузии в рамках модели оцениваются по формуле:

Основные уравнения. Ниже приводится использованная при расчетах система уравнений динамики полидисперсного конденсата [2—4]. Учитывается обмен массой, импульсом, теплом между частицами конденсата и газом, а также между каплями разных фракций:

N

а-7 + с1Іу(уЛ) = сііу(^У^)-^м«2-я5 £ *тК„„(1-Ф,*,),

/Я=5+1

дt

Эр5

ді

+ £1іу(у5рі) <1іу(/)іУрі) + «5 15Ф 15р; р5 ^ К5тФ5т п1

5-І

7=1

N

£

ЯІ=5+1

Эр,у“

Э*

+

сііу(у,р,у8а) = с1іу(/)і5у8аурі) + «д/;)а +иХ{'”ЛИ* +

7=1

+ л*ХА’^Ру(У5"(1"ФлЮ“Р' Е 0СЄ{1’2}’

/я=і+1

+ <1ІУ ( у, р, ея) = СІІУ [р,е;V р,) + п,д, + п3 (Ж„ у, ) + пя (е, + ){т5 }РЬ +

Э*

7=1

N

I

т=я+1

До

-|^ + (11у(ир1,) = с11у(1>ур5;Уа1,)-^{т5}Р11«5, ру=а„р2, рг=р + ру.

Ш 5=1

Здесь е3 = смТ3 + ^-----удельная энергия частиц фракции; ру — плотность водяного пара

в несущем потоке; {/й5 }РЬ — скорость изменения массы частицы сорта 5 при фазовом переходе.

Результаты расчетов. Приводятся результаты расчетов для тонкой плоской пластины (длина в направлении вниз по потоку — 10 м, толщина — 6 см), обтекаемой сверхзвуковым (М^ = 2

и 4) двухфазным потоком под углом атаки 10°. Расчетная область взята прямоугольной, охватывающей треть длины пластины вверх по потоку и приблизительно две длины вниз по потоку. Для уменьшения влияния неточностей в моделировании краевых условий на границах области ее размеры подбирались так, чтобы все косые скачки уплотнения и веера волн разряжения сносились на ту ее часть, через которую поток вытекает. Область покрывалась сеткой из 533 х 322 прямоугольных ячеек со сгущениями вблизи поверхности тела. Обратное влияние конденсата на газ не рассматривалось, поэтому сначала проводились расчеты полей только газодинамических величин течения. На пластине использовались стандартные граничные условия прилипания, непро-текания и адиабатичности. Вычисления совершались в два этапа: сначала методом установления получалась картина обтекания в невязкой постановке с условиями непротекания на теле, затем эти результаты принимались за начальные для уравнений д - ю модели, которые также решались методом установления. Для удовлетворительного разрешения пограничного слоя предварительно проделаны оценки его толщины (по условиям в набегающем потоке) в рамках задачи о ламинарном обтекании бесконечно протяженной вниз по потоку плоской пластины [7]. Мелкость сетки вблизи поверхности тела подбиралась так, чтобы в окрестности его задней кромки на полученную в оценке толщину пограничного слоя приходилось порядка пяти ячеек. Метод крупных частиц позволил проводить сквозной счет, без предварительного выделения особенностей течения. С целью установления влияния различных входных данных задачи проведена серия параметрических расчетов. Здесь целесообразно провести сравнение численных результатов для разных физических ситуаций.

Сначала рассмотрены характерные особенности поведения капель разных фракций конденсата в рамках фиксированного режима обтекания. На рис. 1 представлены изолинии числовых концентраций двух фракций дисперсной фазы для режима Мтс =2. В носовой части пластины изменения концентраций частиц обоих сортов по сравнению с невозмущенными значениями приблизительно одинаковы (увеличиваются в несколько раз за фронтом отошедшей ударной волны). При этом области, где такие изменения выражены, для крупных капель имеют геометрически более компактный характер. Это можно объяснить большей инерционностью последних, вследствие которой их траектории искривляются лишь в зонах более интенсивных изменений параметров газодинамического фона.

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

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

Рис. 1. Изолинии поля числовой концентрации (и-10 6 м 3Режим М00=2, угол атаки а = 10° (стрелка указывает направление набегающего потока):

а — мелкая фракция

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

Также представляется важной роль коагуляции частиц этой фракции при взаимных столкновениях, которая усиливается в результате их вовлечения в развитое турбулентное движение газа, что сопровождается ростом их радиуса (рис. 3,6). Скорость изменения температуры частицы обратно пропорциональна положительной степени ее радиуса, поэтому частицы мелкой фракции в пристеночной области хвостовой части пластины практически полностью испаряются в горячем пограничном слое. Как следствие, образуется свободная от них зона, распространяющаяся вниз по потоку приблизительно на десятую часть длины обтекаемого тела. В принятой при расчетах модели коагуляции рост капель данного сорта возможен при соударениях с более мелкими. Частота столкновений частиц прямо пропорциональна скорости пульсаций. В рассматриваемой зоне последняя увеличивается на несколько порядков по сравнению с фоновым значением. В результате на рис. 3 наблюдается увеличение радиуса частиц крупной фракции на десятки процентов на подветренной стороне хвостовой части. Графики рис. 4, а позволяют сравнить степени релаксации по скорости к несущему газу для частиц двух фракций при разных режимах обтекания. В целом можно отметить, что чем частицы мельче, тем они лучше «подстраиваются» к потоку.

Рис. 2. Траектории частиц мелкой (слева) и крупной (справа) фракций дисперсной фазы (входные данные — см. рис. 1): а — носовая часть пластины; б — хвостовая часть пластины

Х.н

♦

Рис. 3. Изолинии полей радиусов капель

223К 223К

Рис. 4. Зависимость горизонтальных составляющих скоростей (а) и температур (б) дисперсной фазы и газа от вертикальной координаты в носовой части пластины:

-------фракция с радиусом «на бесконечности» 1 мкм; - --10 мкм;----несущий газ

Меньшие отличия горизонтальных составляющих скоростей капель от скорости газа для режима М^ = 4 по сравнению с режимом М„ = 2 объясняются предварительным дроблением частиц в головной ударной волне. Однако полной релаксации может и не быть, что видно для режима Мм = 2, где разница составляет приблизительно 30% для крупной фракции. На рис. 4, б сравниваются температуры фракций и газа в том же сечении расчетной области. Заметно, что релаксация по температуре выражена меньше, чем по скорости, что обусловлено большим характерным временем процессов теплообмена. Разность температур частиц крупной фракции и газа может быть более 50%.

Интересно сравнить результаты расчетов обтекания пластины двухфазным потоком при различных данных «на бесконечности». На рис. 5 и 6 приводятся изолинии полей числовой концентрации и размеров капель двух фракций конденсата для другого режима: Мте = 4. В целом можно отметить те же характерные особенности, что и в предыдущем случае, хотя следует указать некоторую специфику. Зона полного испарения мелкой фракции утолщается вблизи поверхности тела, начинаясь непосредственно за фронтом головной ударной волны, и оканчивается приблизительно на расстоянии в одну пятую часть длины пластины вниз по потоку от донного среза. Здесь сказывается более сильный, чем в предыдущем случае, нагрев несущего газа в неоднородностях течения и в пограничном слое, что ведет к увеличению тепловых потоков к каплям, их более интенсивному разогреву и быстрому испарению. Из трех имеющихся фракций этот эффект наиболее выражен для самой мелкой.

2 и з и п ' т хы 35

Рис. 5. Изолинии поля числовой концентрации |и -10 6 м 3|. Режим М„ = 4, угол атаки а = 10° (стрелка указывает направление набегающего потока):

а— мелкая фракция

Рис. 6. Изолинии полей радиусов капель

—....................................................................Ъ...............................................................................' .............................................................................................................................^-......................................................................................................................—

Рис. 7. Изолинии поля радиуса крупной фракции (размер «на бесконечности» равен 10 мкм) дисперсной фазы с учетом коагуляции/дробления (слева) и без учета этих процессов (справа).

Режим Мос=4, угол атаки а = 15°

В носовой части пластины в отличие от предыдущего режима наблюдаются зоны с характерным размером порядка нескольких толщин тела, в которых существенно увеличивается числовая концентрация обеих фракций (приблизительно в десятки раз), а размеры уменьшаются на порядок. Этот эффект объясняется влиянием дробления капель при прохождении через отошедшую ударную волну, что не происходило при меньшем числе Маха набегающего потока. Рис. 7 иллюстрирует результаты расчетов радиуса капель крупной фракции при учете процессов коагуляции/дробления и в пренебрежении этими эффектами. Можно отметить различные пространственные картины линий уровня и более компактный характер областей выраженных изменений радиуса в первом случае. Из-за коагуляционного роста крупных капель уменьшение их радиуса в приповерхностных зонах выражено слабее. На основании результатов, приведенных на рис. 7, можно сделать вывод о важности учета процессов коагуляции капель при изучении трансформации их массового спектра в отекающем ЛА многофазном потоке.

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

2. Исследовано влияние коагуляции/дробления капель на трансформацию их массового спектра и показана важность учета этих эффектов.

Автор работы благодарит Стасенко А. Л. за постановку задачи и плодотворное обсуждение результатов, а также Моллесон Г. В. за помощь в проведении расчетов.

Работа выполнена при финансовой поддержке РФФИ, гранты 08-01-00540-а, 07-01-00678-а и 07-08-00732.

ЛИТЕРАТУРА

1.ГилинскийМ. М., Стасенко A. JL Сверхзвуковые газодисперсные струи. —

М.: Машиностроение, 1990.

2. С т а с е н к о А. Л. Физическая механика многофазных потоков. — М.: Изд. МФТИ,

2004.

3. Васенин И. М., Архипов В. А., Бутов В. Г., Глазунов А. А., Трофимов В. Ф. Газовая динамика двухфазных течений в соплах. — Томск: Изд. Томского ун-та,

1986.

4. С т е р н и н JI. Е., Ш р ай б е р А. А. Многофазные течения газа с частицами. — М.: Машиностроение, 1994.

5. С oak ley Т. J. Turbulence modeling methods for the compressible Navier — Stokes equations//AIAA-83-1693.

6. БелоцерковскийО. М., Давыдов Ю. М. Метод крупных частиц в газовой динамике.—М.: Наука, 1982.

7. ЛойцянскийЛ. Г. Механика жидкости и газа. — М.: Наука, 1987.

8. К о г а н М. Н. Динамика разреженного газа. — М.: Наука, 1968.

9. М и л л е р А. Б., М о л л е с о н Г. В., С т а с е н к о А. Л. Механика и оптика сверхзвукового мелкодисперсного потока около освещаемой сферы // Ученые записки ЦАГИ. 2007. Т. XXXVIII, №3—4.

10. Д е р е в и ч И. В. Столкновения частиц в турбулентном потоке // Изв. РАН. МЖГ. 1996. №2.

11. Derevich I. V. Statistical modeling of particles relative motion in a turbulent gas flow // International J. of Heat and Mass Transfer, N 49 (2006).

12. D e r e v i с h I. V. Coagulation kernel of particles in a turbulent gas flow // International J. of Heat and Mass Transfer, N 50 (2007).

13. Pao Y. H. Structure of turbulent velocity and scalar fields at large wave numbers // Phys. Fluids. 1965. V. 8.

Рукопись поступила 27/1112008 г