Содержание к диссертации
Введение
1. Уравнение Больцмана. Модели взаимодействия 9
1.1 Уравнение Больцмана и его модели 9
1.2 Упругое взаимодействие 11
1.3 Неупругое взаимодействие 18
2. Численные методы решения уравнения Больцмана 41
2.1 Регулярные методы 43
2.2 Методы Монте-Карло 45
2.3 Метод прямого статистического моделирования 53
3. Приближенные методы решения задач динамики разреженного газа 65
3.1 Локальные методы 67
3.2 Мостовая схема 70
3.3 Самоподобная интерполяция 72
4. Аналитические решения задач динамики разреженного газа 76
4.1 Свободномолекулярное течение по коротким трубам 77
4.2 Течение разреженного газа между параллельными пластинами 81
4.3 Теплопередача в плоском и цилиндрическом течении Куэтта 86
5. Гиперзвуковая аэротермодинамика
5.1 Плоские задачи 102
5.2 Пространственные течения 111
5.3 Теплопередача в критической точке 125
6. Некоторые приложения
6.1 Падение тел на Землю из дальнего космоса 131
6.2 Падение сферического тела на Марс 146
6.3 Проект “Фобос-грунт” 154
6.4 Проект “Экзо Марс” 167
Заключение 184
Литература 1
- Упругое взаимодействие
- Методы Монте-Карло
- Мостовая схема
- Течение разреженного газа между параллельными пластинами
Упругое взаимодействие
Интеграл столкновений j содержит: во-первых, интегрирование по всем возможным скоростям той частицы, с которой сталкивается данная частица, имеющая скорость \; во-вторых, интегрирование проводится по всем значениям телесного угла рассеяния (%- угол рассеяния или меридиональный угол, є - азимутальный или широтный угол). Интеграл столкновений имеет вид: Кроме вычислительных трудностей, связанных с многократным интегрированием в интеграле столкновений всегда существуют трудности моделирования взаимодействия молекул. При этом требуется выполнить несколько, в общем случае противоречивых условий: 1) взаимодействие молекул должно быть “похоже на настоящее”; 2) скорости молекул после столкновения должны рассчитываться достаточно просто; 3) параметры модели взаимодействия должны быть таковы, чтобы они соответствовали измеряемым в эксперименте параметрам газа. Рассмотрим эти пункты более подробно.
Если решается задача обтекания, и требуется найти силовые характеристики обтекания тела, то взаимодействие молекул газа можно моделировать как взаимодействие упругих частиц. Если же требуется вычислить тепловые характеристики, особенно при высокой температуре, как это бывает в гиперзвуковых течениях газа, и нужно учитывать возбуждение внутренних степеней свободы молекул, то модельное взаимодействие должно быть неупругим. Обычно в реальных расчетах столкновений молекул может быть порядка нескольких миллионов. Если модель взаимодействия такова, что скорости молекул после столкновения вычисляются по сложным формулам (или находятся из решения дифференциальных уравнений), то на вычисление каждого столкновения требуется много машинного времени, и затраты времени на решение задачи становятся неприемлемы (хотя с ростом мощностей ЭВМ условие простоты модели взаимодействия становится все менее критичным).
Модели взаимодействия содержат, как правило, некоторое количество параметров, которые можно определить, косвенно, из экспериментальных данных (например, зависимость модельного коэффициента вязкости от температуры должна соответствовать экспериментальной зависимости для реального газа).
Существует лишь несколько моделей взаимодействия, которые используются в расчетах по методу прямого статистического моделирования, так как они в той или иной степени удовлетворяют требованиям перечисленным выше.
Рассмотрим упругое взаимодействие молекул. Обычно считается, что одна молекула неподвижна (она находится в точке 0), а другая движется со скоростью g. После столкновения она приобретает скорость g . Из закона сохранения энергии следует, что g = g и можно записать , = ge, где е -единичный вектор вдоль g . Поскольку упругое столкновение происходит в одной плоскости, то существенная его характеристика - это угол рассеяния % , который является углом между векторами g и g . Угол є - это угол между некоторой плоскостью и плоскостью столкновения. В общем случае сечение столкновения а зависит от g и %, то есть r = r{g,x)-Рассмотрим несколько частных случаев Твердые сферы
Пусть диаметр сферы d, % = ті - 2, Ъ = d sin у/, db = d cos ц/dy/, o = nd2 Интеграл столкновений J(fJ) = —](ffi-ff1)g sin xdxd d2 Скорости после столкновения 2V i; 2 l 2V i; 2 v ; (1.2) В методе Монте-Карло координаты единичного вектора, равномерно распределенного по сфере вычисляется следующим образом. Если Д случайное число равномерно распределенное в интервале (0,1), то компоненты вектора Є выражаются по следующим формулам cos z = 1-21 ; sin % = 2 (1-%); 8 = 2 (1.3)
Как уже отмечалось, для получения параметров модели, необходимо вычислить, например, коэффициенты переноса и сравнить модельные результаты с результатами экспериментов. Чаще всего для этих целей в случае простого газа используется коэффициент вязкости. Методы вычисления коэффициентов переноса хорошо известны из кинетической теории газов.
Методы Монте-Карло
Параметр а2 может быть выбран произвольно. Как показали расчеты, его величина не оказывает влияния на колебательную релаксацию. Целесообразно определить а2 таким образом, чтобы увеличить температурный диапазон применимости модели, ограничением на применение которой служит условие существования решения уравнения (1.20) А212 _ итп = Ф(Т, g, V, а9, с(Т)) (\ + А2 /2) 2 которое разрешимо при Ф 0.25. Предлагается определять а2 из уравнения, являющегося условием минимума функции Ф по а2 для максимального значения этой функции по g где g является решением уравнения: дФ/dg =0. Квантовомеханический гармонический осциллятор Колебательный импульс J2 вычисляется из уравнения (1.16) J 2 = -m(kg) + [m 2 (kgf - 2т(Щ + AE 2 )f2 (1.28) Для нахождения изменения колебательной энергии молекул АЕа используется формула Ландау-Теллера [25-27] для вероятностей перехода гармонического осциллятора /L+, = (« + !«,„ =«PU, = exp(-f) (1.29)
Эти вероятности перехода получены в предположении адиабатичности столкновений и являются справедливыми при условии Ро,1 1 ( Po,i -вероятность перехода с нулевого на первый уровень ). В этих же предположениях для средней величины энергии, которой обмениваются частицы при условии, что в начальный момент времени все молекулы находятся на нулевом колебательном уровне, можно записать AE =hwP0l (1.30) С другой стороны, для этой же величины из уравнения релаксации получаем АЕ = (1.31) где определена в (1.25). Приравнивая (1.30) и (1.31), находим P0,i(T) Рол(Т) = Z(T)exp(e/T)-\ Осциллятор, находящийся на уровне п, с вероятностью Рп,п+1 переходит на уровень п+1 и его энергия увеличивается на величину rlCQ , с вероятностью Ри,и-7 переходит на уровень и-1 и его энергия уменьшается на величину МО , а с вероятностью Pn,n=l-(Pn,n+i+Pn,n-i) остается на том же колебательном уровне и его энергия не меняется.
Однако, в данном варианте, модель имеет существенный недостаток - не выполняется условие детального баланса, так как подкоренное выражение в (1.28) может быть отрицательным. Если в этом случае предположить, что внутренняя энергия не меняется, то есть увеличивается вероятность Рп п за счет вероятности Рп,п+1, так как подкоренное выражение в (1.28) может стать отрицательным лишь при переходе колебательной энергии на более высокий уровень, то нарушается баланс между прямыми и обратными переходами. Чтобы этого не произошло, необходимо увеличить вероятность Рп,п+1 на некоторый коэффициент q, компенсирующий изменение баланса. Пусть Р\л -новое значение вероятности перехода с 0-го на 1-й колебательный уровень с учетом выполнения условия детального баланса. Вероятность P o,i cвязана со старым значением вероятности Po,i параметром q
Пространственно-однородная релаксация азота. Символы – модель ангармонического осциллятора, линии – классический гармонический осциллятор. 1,4 – поступательная; 2,5 – вращательная; 3,6 колебательная температуры. Доля частиц, or для которых подкоренное выражение в (1.28) отрицательно находится по формуле (1.32). ж g a = 2\a((p)sin((p)cos((p)d(p, ос{ф)= \ fgvdg, g = J/2 кв 1 7W COS здесь - угол между векторами к и g, /=1, если одна частица переходит на более высокий колебательный уровень, /=2, если обе частицы переходят на более высокий колебательный уровень при столкновении, fgv - распределение относительной скорости молекул после столкновения. Вероятность перехода Р оЛ должна удовлетворять условию: P0i=P oi-aP oi, которое с учетом (1.29) дает: q=l/(l-a). Ангармонический осциллятор. Модель диссоциации
Обе выше предложенные модели колебательной релаксации имеют ограничения по максимальной температуре, что делает невозможным их применение при достаточно высоких температурах, когда начинают происходить физико-химические процессы в ударных волнах. Колебательная энергия гармонического осциллятора может принимать сколь угодно большие значения, что не соответствует процессу диссоциации, когда спектр колебательной энергии ограничен сверху. Поэтому возникает необходимость рассмотрения колебательной релаксации на основе более реального межатомного потенциала, например потенциала Морзе: U(r) = D(1-e-/3r)2 (1.33) здесь D - энергия диссоциации, г - отклонение межмолекулярной координаты от равновесного положения. Энергия Еп колебательного уровня п определяется выражением Еп =ha n(1-xe(n + 1)) (1.34) где хе=6кв/40 - параметр ангармоничности. Процесс столкновения молекул будем описывать, так же как и для гармонического осциллятора, но с учетом того, что колебательная энергия молекулы принимает значения (1.34). Это позволяет проделать процедуру коррекции вероятностей перехода для ангармонического осциллятора и использовать для вычисления импульса после столкновения формулу (1.28). Вероятности переходов для ангармонического осциллятора определяются на основе формулы Ландау
Мостовая схема
Таким образом, для расчета методом мажорантной частоты не требуется задание сетки по координатам и шага по времени At, как в схемах 1 и 2. Отметим, что величина объема W0 задается во всех методах своим линейным размером h (в схемах 1 и 2 это размер шага сетки), который должен быть меньше минимальной длины свободного пробега в расчетной области. Число частиц, необходимое для расчета, заранее не известно, поэтому требуется повторные расчеты для разных N. Отметим, также, что время расчета растет линейно по числу частиц, что усложняет повторные расчеты. Тем не менее, метод может найти применение в тех случаях, когда задание сетки затруднено или отсутствуют критерии выбора At.
Метод - счетчик столкновений В [100] на основе рассмотрения стохастических дифференциальных уравнений, которые при определенных условиях переходят в уравнение Больцмана, предложен метод аналогичный методу мажорантной частоты. Вводится шаг по времени At. Передвижение частиц и их столкновения разделены. Алгоритм выглядит следующим образом: 1) Вычисляется количество столкновений за время At Nc=-N(N-\)AmAt; 2) Выбирается случайно пара частиц с номерами і и j; 3) Вычисляются р.. и Лу/Лт . Если р.. h и R Л Лт, то столкновение истинно (то есть частицы меняют свои скорости на скорости после столкновения), в противном случае столкновение фиктивно (то есть частицы своих скоростей не меняют). Полное число столкновений равно N (то есть сумма чисел истинных и фиктивных столкновений ); 4) Все частицы передвигаются в течение времени At и переход на 1). Достоинствами этого метода является существенно более простой алгоритм, чем в методе мажорантной частоты и независимость времени расчета от числа частиц N. Недостатком этого алгоритма является необходимость ввода, кроме параметров схемы h и N - параметра At. Модифицированный метод Берда
В методе Берда (схема 1) расчетная область разбивается на ячейки, линейный размер которых h меньше местной длины свободного пробега. В области задается N частиц в соответствии с начальной функцией распределения. Задается шаг по времени At. Столкновения частиц и передвижения за время At проводятся последовательно. На этапе столкновений выбирается случайно пара частиц / и j. Столкновение проводится в том случае, если R AV/Am. В счетчик времени прибавляется величина (для модели молекул в виде твердых сфер) tc = 2W(aglJNm(Nm-1)) Здесь N - число частиц в ячейке на данном шаге, Я - средняя величина А, т т ч в данной ячейке. Столкновения продолжаются до тех пор, пока сумматор t не станет равным At. В [100] было предложено заменить в формуле для tc gtj на gm, а в остальном алгоритм оставить прежним. В этом случае получается алгоритм в точности соответствующий алгоритму метода - счетчика столкновений. Действительно, если N - число частиц в ячейке на данном шаге, то число столкновений в данной ячейке будет N =-N (N -1)ЯАі = Аф с 2 т \ т / т 1с Алгоритм имеет вид: 1) Задается в расчетной области N частиц в соответствии с начальной функцией распределения; 2) В каждой ячейке (например, с номером к) определяется число столкновений по формуле Nkc=1 Nk(Nk-1%At; 3) Выбирается из Nk в ячейке случайно пара частиц с номерами і и j и с вероятностью Яу/Я происходит истинное столкновение, а с вероятностью 1 - Я Якт фиктивное. Суммарное число столкновений равно Nkc; 4) Все частицы передвигаются со своими скоростями в течение времени At и далее переход на 2). В этом методе к параметрам схемы At и N прибавляется еще размер ячейки h . Этот метод более эффективен, чем первые два, так как более точное задание Лкт уменьшает количество фиктивных столкновений.
Для оптимального нахождения параметров схемы предлагается следующий алгоритм [103]. Для пробного расчета задаются At, N и h произвольно. В процессе расчета вычисляется распределение длин пробега и количества столкновений по расчетной области. Далее строится такая координатная сетка, чтобы в ячейках было одинаковое число столкновений, а размер ячеек, должен быть меньше местной длины свободного пробега. Что касается количества частиц в ячейках, то этот вопрос является достаточно тонким для всех алгоритмов метода ПСМ. С одной стороны, в ячейках должно находиться достаточное число частиц, чтобы произошло столкновение (но не слишком много, так как с ростом числа частиц увеличивается время расчета), то есть не менее двух. С другой стороны, этих частиц должно быть много, чтобы при случайном выборе пары в число партнеров не попали частицы, которые уже участвовали в столкновениях.
Течение разреженного газа между параллельными пластинами
В данной главе решается задача Куэтта с различными температурами стенок. Подробные постановки этой задачи приводятся в [128, 129]. Интерес к данной задаче возрос после установления её принципиальных особенностей [130]: немонотонность потока тепла; изменение направления потока тепла при изменении числа Кнудсена (Кп) в некотором диапазоне отношения температур и скорости движения пластин. Указанные особенности вынуждают вносить некоторые изменения в стандартный метод самоподобной интерполяции.
Одной из простейших задач, для которой, однако, до сих пор не получено точное решение уравнения Больцмана, является задача Куэтта о течении и теплопередаче между параллельными бесконечными пластинами, движущимися друг относительно друга. На этом сравнительно простом течении опробованы почти все известные методы решения уравнения Больцмана. С другой стороны, задача имеет самостоятельный интерес, так как позволяет прояснить характер течения вблизи поверхности тел, обтекаемых разреженным газом.
Сильно разреженный газ
Рассмотрим стационарное плоское течение Куэтта газа в горизонтальной плоскости, которую совместим с плоскостью ху. Ось X направлена вдоль пластин. Газ движется между двумя параллельными плоскими стенками на расстоянии h друг от друга, одна из которых покоится, а другая движется в собственной плоскости с постоянной скоростью W. Таким образом, в задаче задается средняя числовая плотность газа П, температуры стенок Т0 при у = 0, Тх при у = h и скорость стенки W при у = И. Требуется найти распределение плотности, скорости, температуры газа между пластинами, а также напряжение трения и тепловой поток к пластинам. При Кп « 1 частота столкновений между частицами газа существенно меньше, чем частота их столкновений со стенками. В этом случае столкновениями между частицами пренебрегается и для описания движения газа применяется свободномолекулярная модель. Тогда функции распределения отраженных молекул по скоростям будут максвелловскими с температурами соответствующих стенок.
Здесь ,,, , z—декартовы компоненты скорости молекул, причем направлена перпендикулярно граничным поверхностям, = :++ ,z — полная скорость молекул. 7"0,7j—температуры нижней и верхней пластин, W — скорость движения верхней пластины. Параметры п0 и щ определяются из условия непротекания газа на стенках (равенства нулю нормальной составляющей скорости газа на стенках). Вычисляя плотность газа и приравнивая нормальную составляющую скорости газа нулю получаем уравнения для определения п0 ипх.
Тепловой поток направлен от нижней пластины к верхней при t 1 + S2/2,(q+cm q-vm), а от верхней пластины к нижней при ґ 1 + S2/2, (q+cm q cm). Величина теплового потока - немонотонна в зависимости от отношения температур. Вязкий газ При Кп«1 в качестве модели сплошной среды примем модель вязкой несжимаемой жидкости. Для течения Куэтта в той же постановке как и в случае свободномолекулярного течения система уравнений вязкой несжимаемой жидкости имеет вид [131]
Таким образом, тепловой поток складывается из диффузионной части (газ покоится, = 0) и конвективной части (температуры стенок одинаковы, Т0 = Т1). Если W = 0 и Т1 Т0, газ покоится, тепловой поток направлен к стенке у = 0, то есть от горячей стенки к холодной. Если Т0 = Т1, то тепловой поток всегда направлен к соответствующей стенке. При у = 0 тепловой поток отрицателен, а при y = h тепловой поток положителен, то есть, стенки нагреваются за счет трения. Так как на стенке y = h диффузионный и конвективный тепловые потоки направлены в разные стороны, можно подобрать такие отношения температур и скорости, что тепловой поток станет равным нулю или поменяет свой знак. Для стенки у = 0 этого сделать нельзя, так как диффузионный и конвективный потоки направлены в одну сторону(к стенке). Запишем уравнение (4.27) в безразмерном виде q = L = Kn —(t-1)--S 2 =Kn(q+ -q ) (4.28) J-HC 8\/2W \ 1 не 1нс J где 7+ =—(7-1), я =— S 2 и z/ 2 0 TTr m „ 15 к T0 hnkT1\ m w \2kT1 4 m T1 m—масса молекулы, к— постоянная Больцмана, п— числовая плотность, Кп— число Кнудсена Величина теплового потока может иметь разный знак в зависимости от отношения температуры пластин. При ґ 1 + 452/15 тепловой поток положителен, а при t 1 + 4S2 /15 тепловой поток отрицателен.
Для напряжения трения можно записать ди__ W /7x =-JU— = -JU — ду И й (4.29) p =p /nkT=-KnS IT НС XV 1 "VI ху Самоподобная интерполяция Поскольку в сильно разреженном газе и в сплошной среде тепловой поток при определенных соотношениях скорости пластины и температур стенок может быть как положительным, так и отрицательным, можно предположить, что тепловой поток может менять знак при изменении числа Кнудсена. Для проверки этого предположения построим приближенное аналитическое решение задачи используя метод самоподобной интерполяции.