Численное решение некоторых задач для модельного кинетического уравнения Больцмана

Как известно, методы расчета течения, основанные на представлении о газе как о сплошной среде (континууме), не подходят для описания движений разреженного газа, в котором средняя длина свободного пробега молекул между двумя соударениями становится сравнимой с характерным размером рассматриваемой области течения. Для анализа таких течений необходимо использовать уравнения кинетической теории газов. В круг задач кинетической теории газов входят, например, задачи об обтекании летательных аппаратов, движущихся на больших высотах, о движении газов в вакуумных аппаратах, ультразвуковых колебаниях в газах, неравновесных течениях и т. д. Как правило, аналитические методы решения применимы только в специальных (упрощенных или предельных) ситуациях. В общем случае необходимо обращаться к численным методам решения.

В последние десятилетия выделились два подхода к численному моделированию течений разреженного газа: статистический метод Монте-Карло и численное решение основного уравнения механики разреженного газа — уравнения Больцмана, и его различных аппроксимаций.

Вероятно, наиболее популярным способом моделирования течений разреженного газа является метод Монте-Карло, который направлен на то, чтобы обойти прямое численное решения уравнения Больцмана и моделировать непосредственно движение фиксированного числа частиц в ячейках численной сетки. Популярность метода Монте-Карло также в значительной степени связана с его физической наглядностью и ростом быстродействия современных вычислительных машин. Получили распространение модификации метода предложенные Бердом [10] и Белоцерковским и Яницким [8, 9]. В нашей стране расчетами по методу Монте-Карло занималась группа Когана М. Н. в ЦАГИ им. Н. Е. Жуковского, в которую также входили Ерофеев Л. И., Перепухов В. А., Горелов С. Л. и другие.

Второй подход к моделированию течений разреженного газа основан на численном решении уравнения Больцмана с точным или приближенным интегралом столкновений. В нашей стране этим занимается группа Вычислительного Центра Академии Наук РФ, в которую входят (или входили) Аристов В. В., Жук В. И., Ларина И. Н., Лимар Е. Ф., Попов С. П., Рыков В. А., Шахов Е. М, Черемисин Ф. Г. и другие. Основные результаты опубликованы в сборниках [55, 56, 57, 58, 59] и множестве статей.

Для одноатомного газа уравнение Больцмана для функции распределения молекул по скоростям f{t, r,?) в отсутствие внешних сил в общепринятых обозначениях имеет следующий вид [14, 25, 67]: + = J (ffi — ffi)9d*d& = J (/, fl, (B.l).

Здесь t — время, r = (a?i, = y, z) — координата физического пространства, (?x, ~ вектор молекулярной скорости, J (f,?) больцмановский интеграл столкновений, da = bdbde, д = — При известной функции распределения макропараметры газа выражаются в виде интегралов по молекулярной скорости. Для числовой и массовой плотности п, р, вектора средней скорости u = (ui, и2, щ), составляющих тензора напряжений Р^, давления р, температуры Т, потока энергии Е{ и потока тепла имеем п = j nui = J Mij = m J €i?jfd4,.

1 p.

Pij = Mij — рщщ, p = -Paa, P = mn, RT = -,.

3 P (B.2).

Ei = m jqi = ~m J- «,-)" — и)2 / d?, ф d? = d€xd?yd€z, г = 1,2,3 где га — масса молекулы.

Кинетическое уравнение Больцмана представляет собой сложное нелинейное интегро-дифференциальное уравнение, аналитическое решение которого возможно лишь в специальных случаях. Трудность численного решения состоит в высокой размерности уравнения и запутанной структуре интеграла столкновений. Заметим, что в общем случае размерность интеграла столкновений равна пяти даже для пространственно — одномерных задач.

Получили развитие два основных подхода к решению кинетического уравнения Больцмана: 1) прямое численное решение непосредственно уравнения Больцмана с точным интегралом столкновений- 2) построение и численное решение аппроксимирующих модельных уравнений. Следует также упомянуть метод моментов [25, 67], который практически не используется в настоящее время.

Один из первых методов прямого численного решения точного уравнения Больцмана был развит в работах [47, 48, 49]. Решения одномерных и двухмерных стационарных задач строились методом итераций по неявной разностной схемедля вычисления интеграла столкновений применялся метод Монте-Карло.

В [6, 7] введено понятие консервативных численных методов для решения кинетического уравнения. Консервативные методы определяются как такие методы, из которых при численном интегрировании (суммировании на дискретной скоростной сетке) конечно-разностной схемы по молекулярной скорости с весами 1, ?2 следуют разностные законы сохранения, т. е. отсутствуют нефизические источники массы, импульса и энергии, вызванные ошибками численного решения. Там же разработан метод расщепления кинетического уравнения по физическим процессам на оператор свободного разлета и оператор релаксации, и на этой основе построен консервативный на каждом шаге алгоритм решения нестационарных задач.

В работе [50] предложен консервативный способ вычисления интеграла столкновений, который обеспечивает выполнение законов сохранения без использования искусственной процедуры коррекции функции распределения. В дальнейшем данный способ вычисления интеграла столкновений был включен в алгоритм расщепления по физическим процессам, причем для оператора свободного разлета применялась монотонная схема второго порядка точности. Полученный численный метод показал свою эффективность [38, 39, 40, 51].

Метод решения уравнения Больцмана, основанный на расщеплении со вторым порядком точности на этап свободно-молекулярного движения газа и этап релаксации газа, предложен в [31, 32]. Новизна метода состоит в линеаризации интеграла столкновений около локально-максвелловской функции распределения на этапе релаксации. На этапе свободно-молекулярного движения используется монотонная схема второго порядка аппроксимации. Для обеспечения консервативности применяется процедура коррекции функции распределения.

Дискретные модели для уравнения Больцмана были предложены в работах [13, 14, 15], для которых доказаны ряд свойств, таких как консервативность и отсутствие нефизических инвариантов.

Альтернативой прямому численному решению уравнения Больцмана является метод построения последовательности модельных кинетических уравнений, которые аппроксимируют уравнение Больцмана [60, 67, 71]. Заметим, что при расчетах течений разреженного газа основной интерес представляет не функция распределения, а макроскопические параметры газа, такие как плотность, скорость, давление, то есть представляет интерес некоторое число первых моментов функции распределения. Поэтому для практических расчетов можно заменить точное уравнение Больцмана некоторым приближенным уравнением, таким чтобы макропараметры газа, вычисленные с помощью этого приближенного уравнения, совпадали или были бы близки с макропараметрами, вычисленными с использованием точного уравнения.

Изложим коротко основную идею построения аппроксимирующей последовательности модельных уравнений, предложенную в [60, 67].

Рассмотрим наряду с точным уравнением Больцмана (В.1) некоторое другое кинетическое уравнение где Q (f, a (r, t)) — приближенный (модельный) оператор столкновений, вид которого будет выбран ниже, а (г, t) — некоторая совокупность макропараметров. Говорят, что уравнение (В.З) аппроксимирует уравнение Больцмана в смысле моментов, если некоторое число первых мо-ментных уравнений от аппроксимирующего и точного кинетических уравнений совпадают. Так как дифференциальные части уравнений (В.1) и (В.З) одинаковы, необходимое условие для совпадения момент-ных уравнений состоит в равенстве соответствующих моментов от приближенного и точного интегралов столкновений:

J 0(Sq (/, I «с®-, *М= J ^(Й Af, 0.

Уравнение (В.З) называется аппроксимирующим уравнением п-го приближения, если условия (В.4) выполнены для всех многочленов п-го порядка от компонент молекулярной скорости. Если в аппроксимации используется только часть многочленов п-го порядка, то такая аппроксимация называется аппроксимацией неполного п-го приближения.

При заданном виде приближенного оператора столкновений из равенств (В.4) определяется совокупность параметров a (r, t). Однако вид самого приближенного оператора столкновений остается в значительной степени произвольным. Для выбора конкретного вида Q (/, а (г, ?)) перепишем уравнение Больцмана в следующем виде: at + идха nJ л = J figdadtu /+ = A J+ = J f’f[gdcTd4i.

B.5).

Здесь и — частота столкновений, J+ - интеграл обратных столкновений. По аналогии с (В.5) запишем приближенный оператор столкновений в релаксационной форме.

Q = *,(/+ /) (В.6) где теперь v — приближенная частота столкновений и /+ - приближенный интеграл обратных столкновений, деленный на частоту столкновений.

В дальнейшем широкое распространение получила модель неполного третьего приближения, получившая название модели Шахова, или S-модели [60, 67] :

0/, Л 0/.

1 + ^(1 — Pr)5acQ (c2 — 5/2) 5 = I.

Здесь R — газовая постоянная, Si — приведенный вектор теплового потока, fm — локально-максвелловская функция распределения, Рг число Прандтля. Коэффициент вязкости /л связан со средней длиной свободного пробега Л зависимостью вида.

5 ц = —тп V2tТКГ, А (В.8).

При Pr = 1 модель Шахова переходит в модель второго приближения, известную как БГК модель, или модель Крука [71]. Заметим, что при переходе к режиму сплошной среды Кп —> 0 б'-модель дает правильные значения коэффициента теплопроводности для одноатомного газа (число Прандтля Рг = 2/3). В этом состоит ее преимущество перед моделью Крука.

Сравнения с расчетами по точному уравнению Больцмана, методу Монте-Карло и с экспериментальными данными показали хорошую точность 5-модели и ее преимущества перед другими моделями [16, 41, 46, 60, 61, 62, 67]. Успешное применение 5-модели стимулировало ее обобщение на течения двухатомного газа с учетом вращательных степеней свободы, и развитие соответствующего численного метода [33, 34, 41, 42].

В целом, модельное кинетическое уравнение значительно проще для численного решения на вычислительных машинах, чем точное уравнение Больцмана. Во-первых, модельное уравнение не содержит пятимерный интеграл столкновений. Во-вторых, для пространственно-од-нои двухмерных задач специальная структура модельного интеграла столкновений позволяет понизить размерность задачи, проводя интегрирование по одной или нескольким компонентам вектора молекулярной скорости? и вводя взамен полной функции распределения / пару редуцированных функций распределения фиф [74].

Тем не менее, несмотря на относительную простоту,^{-модель} является сложным интегро-дифференциальным уравнением высокой размерности. Заметим, что модельный интеграл столкновений является функцией от макропараметров газа, которые определяются интегрированием функции распределения в пространстве скоростей.

Остановимся коротко на численных методах решения модельных уравнений.

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

Для решения стационарных задач часто применяется метод стационарных итераций, который основан на неявной схеме следующего вида [67]:

Йг = ^-D^t-I) — /") (R9) где к — номер итерации. При использовании данной схемы не требуется запоминание значений функции распределения, что ведет к существенной экономии машинной памяти. После задания граничных условий на границах области и поверхности тела для каждого значения молекулярной скорости из квадратурной сетки осуществляется маршевый расчет значений функции распределения в соответствие с направлени-—&diamsем вектора? и вычисляется вклад в интегральные суммы. Примеры решения плоских и осесимметричных задач по методу стационарных итераций при больших и умеренных числах Кнудсена могут быть найдены, например, в работах [20, 37, 41, 42, 63, 64, 65, 66, 43, 91, 92,93, 95]. Метод стационарных итераций также применим к уравнению Больцмана с точным оператором столкновений [47].

Заметим, что при расчетах течений с большими числами Кнудсена для достижения удовлетворительной точности часто необходимо учитывать разрывы функции распределения [67], что существенно усложняет численный алгоритм и его программную реализацию. Особенно чувствительны к учету разрывов области пониженной плотности газа. Подавляющее большинство расчетов проводилось по схемам первого порядка аппроксимации [55, 56, 67]. Учет разрывов в конечно-разностной схеме производился, например, в работах [37, 63, 69]- при этом вычисление интегралов по молекулярной скорости производилось без учета разрывов.

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

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

Третья трудность численного решения кинетического уравнения связана с необходимостью использования полностью консервативных методов решения. Если при больших числах Кнудсена для уменьшения ошибки в выполнении законов сохранения достаточно сгустить сетку по молекулярной скорости, то при малых числах Кнудсена такой подход требует неоправданно подробных сеток.

Один из методов расчета стационарных течений при малых числах Кнудсена был развит в работах [11, 12, 17, 21]. Основная идея предложенного метода состоит в том, чтобы при малых числах Кнудсена основные макропараметры находить из системы законов сохранения, с поправками на добавочные напряжения и потоки, определяемые из кинетического уравнения. Другой способ обойти проблему плохой сходимости метода итераций при малых чисел Кнудсена состоит в использовании нестационарного уравнениястационарные решения при этом получаются методом установления по времени.

Методы второго порядка аппроксимации предложены в работах [19, 28]. Особое внимание уделялось расчету граничных условий со вторым порядком точности. Заметим, что в цитированных работах численный алгоритм включает в себя линейную конечно-разностную схему. Конечно-разностная схема называется линейной, если она имеет постоянные коэффициенты для линейного уравнения с постоянными коэффициентами. Как следует из знаменитой теоремы Годунова [18], не существует монотонных линейных схем второго и выше порядка аппроксимации. Следовательно, расчет течений с большими градиентами по данным схемам может приводить к появлению нефизических осцилляций.

Для получения решений свободных от осцилляций необходимо использовать нелинейные методы. Такие методы для уравнений газодинамики построены во множестве работ, например [26, 72, 77, 81, 87]. В основном они представляют собой обобщение известной схемы Годунова [18]. Применение нелинейных методов в расчетах по кинетическому уравнению пока не получило широкого распространения. Заметим, что прямое применение данных методов к решению кинетического уравнения затруднено из-за наличия сложного нелинейного источника в правой части и другого способа постановки граничных условий. Первые примеры таких расчетов приведены в работах [32, 39, 44, 45, 51, 86]. Использовавшиеся схемы обеспечивали второй порядок аппроксимации внутри расчетной области и по-видимому только первый порядок на границе.

Построение консервативных численных алгоритмов представляет из себя сложную задачу вычислительной математики из-за высокой размерности кинетического уравнения и наличия сложного нелинейного источника в правой части. Для точного уравнения Больцмана задача построения консервативных алгоритмов была недавно решена в работе [13, 14, 15, 50]. При решении модельных уравнений в настоящее время используются различные процедуры коррекции граничных условий, функции распределения и модельного оператора столкновений [6, 19, 35, 36]. Также упомянем работу [31, 32], где процедура коррекции применялась при решении линеаризованного уравнения Больцмана. Несмотря на хорошие результаты, полученные при использовании данных коррекций, они все же представляют собой в значительной степени искусственный прием обеспечения консервативности метода.

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

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

Метод исследования. Применяются методы механики разреженных газов и вычислительной математики.

Научная новизна.

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

Предложена процедура учета разрывов функции распределения при вычислении интегралов по молекулярной скорости для двухмерного стационарного кинетического уравнения.

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

Рассчитано поперечное обтекание пластины с различными граничными условиями гиперзвуковым потоком разреженного газа для чисел Маха 2 < М < 30 и числа Кнудсена Kn > 1.

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

Построен полностью консервативный численный метод третьего и выше порядка аппроксимации для двухмерного нестационарного кинетического уравнения с модельным оператором столкновений в виде S-модели.

Рассчитано поперечное обтекание пластины сверхзвуковым потоком разреженного газа при числе Кнудсена Кп = 5 х Ю-3 и чисел Маха набегающего потока М = 5, 10 для теплоизолированной пластины и М = 5 для холодной пластины.

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

На защиту выносятся следующие положения.

1. Расчет нестационарного испарения с плоской поверхности конденсированной фазы для всех режимов испарения произвольной интенсивности.

2. Численный алгоритм решения двухмерного стационарного уравнения с учетом разрывов функции распределения не только в конечно-разностной схеме второго порядка аппроксимации, но и в процедуре вычисления интегралов от функции распределения по молекулярной скорости.

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

4. Расчет поперечного обтекания пластины гиперзвуковым потоком разреженного газа при граничных условиях теплоизолированной и холодных пластин.

5. Метод построения консервативных численных методов для трехмерного кинетического уравнения с оператором столкновений в виде S-модели без коррекции функции распределения.

6. Консервативная схема повышенного (от третьего и выше) порядка аппроксимации для двухмерного кинетического уравнения с оператором столкновений в виде S-модели.

7. Расчет поперечного обтекания пластины сверхзвуковым потоком разреженного газа для числа Кнудсена Кп = 5 х Ю-3 при граничных условиях теплоизолированной (числа Маха 5 и 10) и холодной (число Маха 5) пластины.

Аппробация работы.

Результаты диссертационной работы докладывались на:

• Российской национальной конференции по теплообмену, Москва, МЭИ, 1998 г.

• Пятом Международном совещании-семинаре 'Инженерные проблемы новой техники', Москва, МГТУ им. Н. Э. Баумана, 1998 г.

• Второй международной конференции по неравновесным процессам в соплах и струях, Санкт-Петербург, 1998 г.

• II Всероссийской научно-технической конференции молодых ученых 'Современные проблемы аэрокосмической техники', Жуковский, ЦАГИ им. Н. Е. Жуковского, 1999 г.

• Международной научно-технической конференции молодых ученых 'Современные проблемы аэрокосмической науки и техники', Жуковский, ЦАГИ им. Н. Е. Жуковского. 2000 г.

• Семинаре Сектора кинетической теории ВЦ им. А.А. Дородни-цина РАН, Москва, 2003 г.

• Семинаре по математической физике ИПМ им. М. В. Келдыша РАН, Москва, 2003 г.

• Семинаре ФАЛТ МФТИ, Жуковский, 2003 г.

Публикации. Численный метод решения во всех разделах разработан автором. Постановка задач во всех разделах принадлежит Е. М. Шахову. По теме диссертации опубликовано 10 работ [88−97].

Структура и объем работы. Диссертация состоит из введения, трех разделов, заключения и списка литературы. Диссертация изложена на 171 странице, содержит 47 рисунков. Библиография содержит 97 наименований.

Заключение

.

В диссертации для кинетического уравнения Больцмана с модельным оператором столкновений разработаны новые методы повышенной точности. С применением данных методов решен ряд одномерных и двухмерных задач.

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

2. Для двухмерного стационарного кинетического уравнения реа-^ лизован численный алгоритм решения с полным учетом разрывов функции распределения. Процедура учета разрывов функции распределения при интегрировании по молекулярной скорости является новой.

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

4. В задаче о поперечном обтекании пластины показано, что течение в зоне донного вакуума определяется граничным условием на пластине. В частности, плотность газа за холодной пластиной по крайней мере на порядок выше, чем за теплоизолированной пластиной. Для холодной пластины затекание газа в область за пластину имеет место даже для очень больших скоростей потока и чисел Кнудсена. Т.

5. В третьей главе диссертации для нестационарного трехмерного кинетического уравнения предложен новый метод построения полностью консервативных численных методов, который обеспечивает консервативность расчета для любых чисел Кнудсена и граничных условий диффузного отражения, а также правильную аппроксимацию вектора потока тепла.

6. Для двухмерного нестационарного модельного уравнения предложен численный алгоритм повышенного (от третьего и выше) порядка аппроксимации как про пространству, так и по времени.

7. Построенным консервативным методом рассчитано поперечное обтекание пластины потоком разреженного газа для граничного условия теплоизолированной (числа Маха 5 и 10) и холодной (число Маха 5) пластины для числа Кнудсена Кп = 5 х Ю-3. Исследовано влияние граничных условий на картину течения. Показано, что в целом кинетическое решение хорошо согласуется с решением уравнений Навье-Стокса везде кроме области донного вакуума за пластиной.