МЕХАНИКА ЖИДКОСТИ И ГАЗА <5 • 2008
УДК 533.6.011.8
© 2008 г. О. И. РОВЕНСКАЯ
ПРЯМОЕ ЧИСЛЕННОЕ МОДЕЛИРОВАНИЕ ЭВОЛЮЦИИ ДВУМЕРНОЙ ВИХРЕВОЙ СИСТЕМЫ В РАЗРЕЖЕННОМ ГАЗЕ
В рамках кинетического подхода решается двумерная задача с начальными условиями Тейлора-Грина и с периодическими граничными условиями в вязком сжимаемом слаборазреженном газе. Цель данной работы - моделирование эволюции заданной вихревой системы на основе прямого численного решения уравнения Больцмана. Для этого использовался метод дискретных ординат, причем интеграл столкновений вычисляется консервативным проекционным методом Ф.Г. Черемисина, сохраняющим плотность, импульс и энергию. Полученные данные позволяют проследить эволюцию вихревой системы, заданную начальными условиями, и определить спектральные свойства течения. Представлены распределения параметров потока в последовательные моменты времени.
Ключевые слова: вязкий сжимаемый газ, уравнение Больцмана, вихревой каскад, метод дискретных ординат.
Уравнение Больцмана [1] предлагает иную модель по сравнению с уравнениями На-вье-Стокса и дает новые возможности для описания течений и явлений. Решение задач с помощью кинетического подхода позволяет описать механизмы, недоступные уравнениям Навье - Стокса. Примеры использования данного подхода можно найти в [2-5].
Рассматривается эволюция вихревого каскада - дробление крупных вихрей на мелкие, что является основной характеристикой турбулентности. Изучение каскадных вихревых явлений в разреженном газе - актуальная задача, так как размер мелких вихрей в турбулентных течениях может превосходить длину свободного пробега в газе на 1-2 порядка величины [6]. Моделирование вихревых течений при малых конечных значениях числа Кнудсена позволяет исследовать режимы, близкие к сплошной среде, но с учетом более общего диссипативного механизма [3-5].
В настоящей работе изучается влияние интенсивности начальных условий и числа Кнудсена (степени разреженности газа) на эволюцию заданной вихревой системы. Следует отметить, что на вихревую динамику оказывают влияние различные факторы, в первую очередь небаротропность и вязкость. Влияние этих факторов на спектральные свойства течения может быть оценено при использовании различных моделей среды. Для этой цели помимо уравнения Больцмана используется система уравнений Эйлера.
Кроме того, большое внимание уделяется выявлению спектрального закона, поскольку основной показатель развития вихревого каскада - распределение спектральной плотности кинетической энергии в пространстве волновых чисел. Спектральные законы в турбулентных течениях представляют большой интерес, они изучаются теоретически, экспериментально и численно [2-4, 7, 8].
1. Постановка задачи. В работе рассматривается двумерная задача о течении разреженного газа с начальными условиями Тейлора - Грина и с периодическими граничными условиями по х и у [2, 4, 9].
В момент времени ? = 0 задаются распределения безразмерных плотности, компонент скорости и температуры
р0(x, y) = 1 + Csin(2пx) sin(2ny)
и0(x, y) = A sin(2nx) cos(2ny)
(1.1)
v0 (x, y) = B cos (2nx) sin (2ny)
T0 (x, y) = 1+ D cos (2nx) cos (2ny)
Расчетная область представляет собой квадрат, на границах которого задаются периодические граничные условия: F(x + L, y + L) = F(x, y), 0 < x < L, 0 < y < L. В начальном поле задавалась система из четырех вихрей. Время моделирования достигало t = 0.6.
Начальное число Рейнольдса определяется по формуле Re = Vyn/2 Mmax/Kn¿, где Mmax - максимальное значение числа Маха в начальном поле, Кпг = Kn/p - локальное число Кнудсена для невозмущенного газа, у = 5/3 - показатель адиабаты для одноатомного газа.
Динамику заданной вихревой системы можно описать следующим образом: при развитии начального поля течения во времени происходит возникновение вихрей меньшего размера, с передачей им кинетической энергии и последующей ее диссипацией в тепло.
2. Численный метод. Для решения задачи используется метод дискретных скоростей. Для аппроксимации интеграла столкновений используется метод, разработанный Ф.Г. Черемисиным [10].
В качестве масштабов для обезразмеривания уравнения Больцмана использовались:
v = JRT = Vt/42 - масштаб скорости, где Vt = J2R T0 - тепловая скорость; L - масштаб длины, некоторый характерный размер течения; т = L/v - масштаб времени. В результате уравнение Больцмана для одноатомного газа примет вид
1+< I = T¿nJ (f•f >
где f=f(t, x, £) - функция распределения, £ - вектор скорости частицы, J(f, f) - интеграл столкновений, Kn = X/L - число Кнудсена, X = 1/( J2 пп0 bf) - длина свободного пробега, bf - эффективный радиус, bef = bm для случая твердых сфер.
Численный метод основан на схеме расщепления, в рамках которой раздельно вычисляются изменения плотности распределения за счет столкновений и переноса. Для численного решения уравнения Больцмана производится дискретизация пространства скоростей с помощью соответствующей расчетной сетки. Интеграл столкновений вычисляется как сумма вкладов от различных пар точек скоростного пространства. Однако в общем случае скорости частиц после столкновения не попадают в узлы скоростной сетки, поэтому необходима регуляризация распределения вкладов в узлы сетки при столкновениях. Такой подход обеспечивает консервативность алгоритма, отсутствие искажений плотности распределений, а также равенство нулю интеграла столкновений в случае максвелловской функции распределения.
Для аппроксимации переноса применялась явная консервативная по потокам схема второго порядка точности в рамках метода конечного объема (TVD - подход) с ошибкой O(At, Ax2) при ограничениях TVD метода [11]. Шаг по времени ограничивался условием Куранта: At = CFL min(Ax/Vmax, Ay/Vmax), где CFL - число Куранта, Ax, Ay - размер ячейки.
Как показали расчеты, удобнее использовать симметричную схему, включающую в себя три шага: перенос в течении At/2, релаксацию на промежутке времени At и еще один перенос в течении At/2. В этом случае условие Куранта - CFL < 1 для схемы второго порядка точности. Кроме того, этот метод позволяет использовать только два массива переменных на временных слоях t и t + At.
Если задаваемая область в пространстве скоростей достаточно большая и шаг сетки, заданной в физическом пространстве, больше или равен длине свободного пробега и условие Куранта выполнено, тогда гарантируется выполнение условия At < т0, где т0 = A/v - среднее время между столкновениями.
Для улучшения эффективности алгоритма в работе проводится распараллеливание в физическом пространстве. Программа на языке C++ была написана с применением библиотеки MPI (Message Passing Interface) [12, 13]. Расчеты проводились на многопроцессорной системе MBC - 1000/16. Компьютерное время составило от 96 до 1216 ч в зависимости от начальных условий.
Вычисления в работе проводились при следующих параметрах: -5 < bjv < 5 - область в скоростном пространстве, шаг по пространству скоростей менялся A^/v = 0.33-0.4. Число узлов в физическом пространстве nx = ny = 112 для Kn = 0.01; nx = ny = 208 для Kn = 0.005; nx = ny = 448 для Kn = 0.0025. Шаг по времени ограничивался условием CFL = 0.8.
При расчетах в основном применялась симметричная схема, так как она сокращает число массивов данных, необходимых для вычислений, и требования по памяти. Для сравнения были проведены расчеты, дополнительно включающие схему Рунге-Кутта второго порядка точности, которые показали, что на небольших временах счета ошибка использования только симметричной схемы незначительна - порядка 1-2%. Уменьшение шага дискретизации в пространстве скоростей не оказывало существенного влияния на картину течения.
3. Построение спектра. Спектральная плотность кинетической энергии изотропного турбулентного течения не зависит от направления волнового вектора и в двумерном случае определяется в виде
E( k ) = E( k ) link
где E(k) - многомерная спектральная плотность кинетической энергии, k - волновой вектор. Построение спектра кинетической энергии E(k) по волновым числам k = |k| для конкретного течения является сложной задачей, так как диапазон масштабов вихрей всегда ограничен. Это приводит к необходимости выбора ряда направлений и осреднения спектров.
Для вычисления спектральных амплитуд коэффициенты Фурье находятся по методу наименьших квадратов. Вычислив коэффициенты Фурье разложения одномерного пространственного распределения кинетической энергии по заданному набору функций, можно найти амплитуды Ek, соответствующие волновым числам k
E = ak sin ( kx) + bk cos ( kx)], Ek = Ja2k + b2k
где диапазон изменения волнового числа к = 1, ..., N, где N - размерность расчетной сетки.
Спектральная плотность кинетической энергии в задаче определяется с помощью
осреднениЯ по формуле £сред = (Ек, х = 0 + Ек, х = Ь + Ек, у = 0 + Ек, у = ь)/4, где Ек, х = 0 - спектральная плотность кинетической энергии по направлению х = 0, Ек х = Ь - по направлению х = Ь, Ек, у = 0 - по направлению у = 0, Еку = Ь - по направлению у = Ь.
4. Результаты численного моделирования. В работе приводятся картины распределения плотности в физическом пространстве в последовательные моменты времени. Представлены распределения спектральной плотности кинетической энергии в пространстве волновых чисел. В силу анизотропии течений графики спектральной плотности кинетической энергии получены с помощью дискретного преобразования Фурье и выделения характерных направлений. Начальное распределение плотности, заданное условиями (1.1), представлено на фиг. 1.
k
Фиг. 1. Изолинии поля плотности в момент времени t = 0
Для сравнения были проведены расчеты в рамках уравнений Эйлера на сетках 100 х 100, 200 х 200, 400 х 400 ячеек с помощью конечно-разностного метода типа Годунова второго порядка точности по пространственным переменным и времени (TVD подход и метод Рунге-Кутта). На фиг. 2 представлено сопоставление полученных результатов для начальных условий: А = 0.5, B = - 0.4, C = D = 0.1, Kn = 0.01, Mmax = 0.387 и Re = 63. Такой случай соответствует пределу больших чисел Рейнольдса, если схемная диссипация метода энтропийно согласована.
Сопоставление с данными,
Для дальнейшего прочтения статьи необходимо приобрести полный текст. Статьи высылаются в формате PDF на указанную при оплате почту. Время доставки составляет менее 10 минут. Стоимость одной статьи — 150 рублей.