книги / Метод крупных частиц в газовой динамике

3. Исследуем вязкостные эффекты полученной разностной схемы [205,214, 252]. Для этой цели используем аппарат дифференциальных приближений.

Выпишем для простоты в одномерном случае дифференциальное прибли­ жение уравнения импульса. При этом в качестве формул эйлерова этапа будем использовать выражения (12.47), на лангранжевом этапе — формулы для АМ п первого порядка точности с учетом направления потока, на заключи­ тельном — разностный аналог законов сохранения. Опуская промежуточные выкладки, получим

й Ш + дЛ Е + ^ 1 = (Л + 2 ) + р | и 1 ^ + рм2-£ --§ ■ и | Дх2-

- 1 р | Д х 2- 2 , и ( А + 2 ) | д г } | ^ + 0(Дг2, Дх2). (12.48)

Отсюда видно, что диффузионный коэффициент (выражение при д2и/дх2) мо­ жет в значительной степени зависеть от аппроксимационной вязкости. По­ этому в реальных расчетах эффективное число Рейнольдса зависит от молеку­ лярной н аппроксимационной вязкостей:

Величина схемной вязкости зависит, как известно, от параметров потока в данной точке течения. На рис. 12.3 приведена картина (схема) обтекания

334 РАЗВИТИЕ МЕТОДА КРУПНЫХ ЧАСТИЦ [ГЛ. XII

В обзоре методов для решения уравнений Навье — Стокса [414] при об­ суждении проблемы больших чисел Re А. А. Дородницын отмечает, что «при­ ближение к реальному решению может быть достигнуто, когда Re*/i<e<$l, иначе математическая вязкость превысит реальную и мы не сможем узнать, чему отвечает полученное численное решение».

4. В качестве примера приведем результаты расчета вязких течений по описанной выше методике [205, 214, 352] (см. гл. VII).

Был осуществлен следующий численный эксперимент. В невозмущенный

гидродинамических параметров изменятся и будут соответствовать течению газа, обладающего молекулярной вязкостью и теплопроводностью с опреде­ ленным числом Re. Проведенные расчеты показали, что полученное решение в основном ведет себя указанным образом [214].

Некоторые результаты данного численного эксперимента приведены на рис. 7.41, а, б, где показаны профили плотности р и числа Маха (М) вдоль оси симметрии AO(j—1) (рис. 12.3). Сплошные линии соответствуют течению без молекулярной вязкости, пунктирные — с ее учетом. Параметры вязкого газа были таковы: М „= 2 ,0; р=0,01; А = —0,67; £=2,15. В данном случае расчеты проводились при неизменном во всем поле течения коэффициенте молекуляр-

По приведенной разностной схеме можно проводить счет от произвольных начальных данных (как при решении уравнений Эйлера), в том числе и от первоначально невозмущенного плоскопараллельного потока. Однако при этом время счета будет значительно больше, чем в случае использования не­ вязкого модуля эйлерова этапа, так как введение навье-стоксовских членов накладывает более жесткие ограничения на условие устойчивости используе­

мой явной разностной схемы.

Разностные схемы такого типа предпочтительны при решении нестацио­ нарных аэродинамических задач при наличии разнохарактерных областей взаимодействия (течения со вдувом, струйные задачи, турбулентные движе­

дования различного вида гидродинамических неустойчивостей [320, 329 и др.]. Особый интерес представляет при этом не только начальная линейная фаза развития возмущений, но и развитая неустойчивость.

Исследование рэлей-тейлоровой неустойчивости (РТН) является весьма актуальной задачей. Помимо несомненного теоретического интереса, она имеет большое значение для ряда важных практических проблем, например, при изу­ чении устойчивости сжатия оболочек в задачах лазерного термоядерного син­ теза, при получении сверхсильных магнитных полей и др.

Большой вклад в развитие обсуждаемой проблемы внесли работы Э. Фер­ ми [321—323]. В них рассматривается развитие тейлоровской неустойчивости на границе «жидкость — вакуум» в линейном [321] и нелинейном [322] слу­ чаях, а также на границе двух жидкостей [323]. При этом в нелинейных слу­ чаях приводятся постановки, геометрически близкие к численным. В част­ ности, плавная поверхность раздела аппроксимируется кусочно-линейной ступенчатой функцией.

В развитии РТН прослеживаются следующие ярко выраженные стадии: линейная, промежуточная, регулярная асимптотическая и турбулентная [421, 422]. Линейная стадия характеризуется малостью амплитуды а, по сравнению с длиной волны возмущения L, и экспоненциальной скоростью роста. Когда величина амплитуды возмущения а достигает 0,4L, процесс вступает в стадию, промежуточную между линейной и регулярной асимптотической. На регулярной асимптотической стадии, наступающей при a«0,75L, окончательно форми­ руются пики тяжелой жидкости, проваливающиеся с постоянным ускорением, и пузыри легкой жидкости, всплывающие с постоянной скоростью. Эта стадия РТН является^ неустойчивой [421, 422] и сменяется турбулентной стадией, в ходе которой происходит интенсивное взаимодействие возмущений различ­ ных длин волн и перемешивание жидкости.

РТН наиболее исследована для случая плоской поверхности раздела и стремящегося к бесконечности отношения плотностей тяжелой и легкой жид­ костей. Линейная стадия подробно изучена в классических работах Рэлея, Тейлора и Льюиса [423—425], регулярная асимптотическая — в [426—428], в [429] развита феноменологическая теория турбулентной стадии, а в [422] высказаны некоторые соображения о механизме ее образования.

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

До сих пор во всех как аналитических, так и вычислительных работах исследовался только двумерный случай. Однако двумерная модель неадек­ ватна, вообще говоря, физике явления: в физическом эксперименте двумерные

структуры разрушаются из-за поперечной коротковолновой неустойчивости и превращаются в трехмерные.

Численные методы, примененные в [430, 431], в принципе могут быть рас­ пространены и на трехмерный случай. Однако это ведет к столь значительному увеличению времени счета и возрастанию объема необходимой машинной па­ мяти, что подобные расчеты на современных ЭВМ не реальны.

2. Алгоритм метода крупных частиц для расчета пространственно-трех­ мерных нестационарных течений весьма экономичен и позволяет рассматри­ вать обсуждаемую задачу [213].

Развитие двумерной РТН вплоть до больших амплитуд, когда процесс носит существенно нелинейный характер, впервые изучалось данным методом Ю. М. Давыдовым в [432]. В этой работе показано, что проведенные численные расчеты дают скорость развития неустойчивости на несколько десятков про центов (в зависимости от варианта) меньшую по сравнению с аналитическими оценками [322]. Согласно многочисленным исследованиям тяжелая жидкость должна проникать в легкую с постоянным ускорением. Экспериментальные данные Д. Льюиса [325] подтверждают это положение.

Методом крупных частиц здесь решается полная пространственно-трех­

ние; ускорение свободного падения g направлено по оси у .

Программно данный численный эксперимент был выполнен в рамках па­ кета прикладных программ КРУЧА (крупные частицы) [455], в который бы­ ли добавлены модули, отражающие специфику задачи (гравитацию, начальные данные). Отметим, что хотя в качестве уравнения замыкания системы (12.60) используется уравнение состояния идеального газа 3=р11(к—1)р] (х — по­ казатель адиабаты), сжимаемость мало влияет на характер процесса (значение числа Маха в расчетной области нигде не превышает 0,3).

Расчет производится сквозным образом без выделения контактной по­ верхности раздела тяжелой и легкой жидкостей. Двумерные расчеты показали, что, несмотря на размазывание контактной поверхности раздела (7 8 ячеек в районе пика и 4—5 ячеек в районе пузыря), ее с достаточной точностью можно отождествлять с поверхностью, на которой р=0,5 (pi+Рг)! где Рг плотность тяжелой жидкости, pi — легкой (близкие результаты дает определение этой поверхности разрыва другими эвристическими способами: по шах grad р и т. д.).

Для сравнения с данными других авторов по трехмерной программе были рассчитаны двумерные тестовые задачи. Начальные данные и параметры тестов были идентичны используемым в [431], за исключением того, что рассматри­ валась среда с уравнением состояния идеального газа. Из-за малости числа

мыми средами важно везде соблюдать условие div W ^O (за исключением по­ верхности раздела). Иначе возникнут возмущения плотности, которые могут исказить картину развития РТН.

Двумерные расчеты РТН методом крупных частиц показало хорошее сов­ падение с [431]. На рис. 12.5 приведены формы контактной поверхности при

Рис. 12.6. Динамика формы поверхности раздела при трехмерной рэлейтейлоровой неустойчивости в перспективной проекции.

ра/р1 = 10 в различные моменты времени: 0,25; 0,49; 0,73; 0,97; 1,21; сплошные линии соответствуют расчетам методом крупных частиц, штриховые ре­ зультатам работы [4311. График скорости движения пика также близок к при­ веденному в [431], среднее ускорение пика, как и в [431], составило примерно

к длине волны возмущения RI2L составило 0,40, что ближе к значению 0,39, приведенному в [430], и теоретическому значению 0,35 [426], чем 0,48 в [431].

w = A c o s ^ - s m ~ - [ 2 H ( y )~ 1]ехр ( ~ 2" ' у1- ) .

В остальном все начальные, граничные условия и парагоны задачи были такими же, как и в двумерном случае.

На рис. 12.6, а — г приведены поверхности раздела ь перспективной проекции, построенные в последовательные моменты времени ^=0,6; 0,9; 1,2; 1,5 соответственно. Отчетливо видно всплывание широких пузырей и проваливание достаточно узких пиков.

Заметим, что высоты (размеры по оси у) пузыря и пика на рис. 12.6 в дейст­

Более четкое представление о характере процесса дает рис. 12.7, на котором изображены сечения контактной поверхности в плоскостях х = 0 (z=0) и JC= Z, х = —г при р=0,5(р1+ра) в момент времени *=1,5.

Отношение радиуса кривизны пузыря к длине волны возмущения R/2L составляет 0,34 в сечениях х= 0 и z=0 и 0,31 в сечении х= г. Теоретическое и экспериментальное значения этой величины, полученные при исследовании подъема воздушных пузырьков в вертикальных трубах, заполненных жидко­ стью, составляют 0,35 [122]. Безразмерная скорость подъема пузыря uB(2Lg)-li* равнялась 0,31 ±0,02, в то время как согласно данным, приведенным в [122], для цилиндрических труб она составляет 0,32, а для труб прямоуголь­ ного сечения с отношением сторон прямоугольника 1 : 4 — 0,29. Ускорение пика было примерно на 20% больше, чем в двумерном случае, в то время как сам пик был несколько шире.