Разработка и анализ торцевых поверхностей магнитноразрядного измерителя плотности
Определение концентрации молекул разряженного газа в произвольном объеме. Моделирование набегающего потока, движения молекулы внутри объема. Генерация вектора скорости молекулы и координат точки влета. Моделирование потока собственных газовыделений.
| Рубрика | Физика и энергетика |
| Вид | дипломная работа |
| Язык | русский |
| Дата добавления | 06.07.2011 |
| Размер файла | 1,8 M |
Отправить свою хорошую работу в базу знаний просто. Используйте форму, расположенную ниже
Студенты, аспиранты, молодые ученые, использующие базу знаний в своей учебе и работе, будут вам очень благодарны.
1.2.1 Описание программы моделирования MODMD05
Программа modmd05.pas загружает в среду Паскаль файл исходных данных "inp2.txt" и файл результатов "concentr", устанавливает тип связи программы и файлов и закрывает их после отработки программы.
Все неоднократно повторяющиеся процедуры выделены в программе
в отдельные блоки - подпрограммы - для обеспечения большей компактности программы и повышения ее наглядности.
В основном теле программы находится процедура считывания данных из исходного файла, реализован механизм розыгрыша точек влета в исследуемый объем, расчет коэффициентов, необходимых для формирования матрицы результатов, содержатся условия выполнения некоторых подпрограмм и процедура формирования матрицы результатов с записью ее в файл результатов.
Формирование вектора скорости при отражении молекулы от поверхности осуществляется в подпрограмме scattering. Расчет времен пребывания молекулы в частных объемах, на которые разбит расчетный объем, выполняет подпрограмма countmolec. Расчет траектории движения молекулы в исследуемом объеме выполняет подпрограмма counttrack. Генерация точек, распределенных в соответствии с законом Максвелла, осуществляется в подпрограмме maxwell.
1.3 Описание программы MODMD24
На базе существующей программы моделирования взаимодействия молекул набегающего потока и СВА с конструктивными элементами MODMD05, была создана программа моделирования MODMD24, в которой введен дополнительный коэффициент, характеризующий исследуемый поток молекул: Tmm - температура направленного потока молекул с Максвелловским законом распределения модуля скоростей. При этом появляется возможность дополнительно моделировать набегающий поток, создаваемый собственными газовыделениями КА.
Также введены были дополнительные параметры, характеризующие входные отверстия датчика - центральный угол ugol, задающий размер части кольца, через которую осуществляется влет молекул в исследуемый объем и вылет из него, и коэффициент прозрачности отверстий kpr.
1.3.1 Моделирование потока собственных газовыделений.КА
В программе MODMD24 выражение для расчета модуля направленной скорости молекул, которое в MODMD05 описывается выражением (20), реализуется в виде:
, (52)
где W0 - направленная скорость молекул, нормированная на модуль наиболее вероятной скорости молекулы, имеющей температуру стенок исследуемого объема: Vст02 = RT T;
V0 - модуль направленной скорости молекул, влетающих в исследуемый объем;
T - температура стенок исследуемого объема;
RT - значение газовой постоянной для молекулярного азота;
V - величина, распределенная по закону Максвелла;
Tmm - температура направленного потока молекул с Максвелловским законом распределения модуля скорости.
При значении Tmm=0 выражение (52) сводится к исходному выражению (20).
В результате, новая модификация программы MODMD24 дает возможность моделировать все виды потоков нейтральных молекул, воздействующих на исследуемый МИП:
- при Tmm 0, Tm=0 и V0=0 моделируется воздействие на МИП потоков собственных газовыделений систем КА;
- при Tmm=0, Tm 0 и V0=0 моделируется воздействие на МИП стационарной составляющей собственной внешней атмосферы КА;
- при Tmm=0, Tm 0, V0 > 0 моделируется воздействие разреженной земной атмосферы для высот полета КА на МИП.
1.3.2 Моделирование влета (вылета) молекул в датчик
Программа для моделирования распределения концентрации аэродинамического потока молекул в датчике MODMD05 была модифицирована таким образом, чтобы можно было произвольно изменять размеры отверстий в переднем и заднем торцах. Изменение размеров отверстий регулируется заданием значения центрального угла ugol (см. рис.1.8) в соответствующем файле исходных данных inp24.txt. В связи с изменением размеров входных отверстий в программе MODMD24 изменены формула для расчета площадей отверстий в торцах и механизм задания координат точек влета молекул в исследуемый объем.
Так как отверстия в переднем и заднем торцах одинаковы по форме и площади, интегральный коэффициент прозрачности торцевых стенок определяется выражением:
(53)
где r01, r02 - радиусы, определяющие размеры отверстий в переднем и заднем торцах МИП;
ugol [град.] - угол, определяющий размеры отверстий МИП (см.рис.1);
kpr - коэффициент прозрачности отверстий МИП.
Рис.1.8 Модель датчика с учетом параметра ugol
Выбор точек влета-вылета молекул в исследуемый объем определяется геометрическим положением открытой области. Розыгрыш точек влета происходит в несколько этапов:
а) В зависимости от знака компоненты скорости vx, определяется торец, через который молекула влетает в исследуемый объем, а, следовательно, координата x2.
б) При помощи генератора случайных чисел выбирается случайное число xkpr в диапазоне (0,1), которое характеризует вероятность того, что молекула, оказавшись в области отверстия, проникнет сквозь него в исследуемый объем.
в) При помощи генератора случайных чисел выбираются значения координат y2 и z2 точки влета в диапазоне (-1,1) до тех пор, пока выбранная точка не окажется в области одного из отверстий. Для этого необходимо, чтобы совместно выполнялись условия:
y22 + z22 < r022, (54)
y22 + z22 > r012, (55)
(при sk = -1) (56)
или
(при sk = 1) (57)
xkpr < kpr
где sk - параметр, принимающий значение "-1" при розыгрыше влета в передний торец, и значение "1" при розыгрыше влета в задний торец.
Первое слагаемое в выражениях (56) задает косинус половины угла, определяющего размер входного отверстия (в радианах), второе -косинус угла, под которым предполагаемая точка влета молекулы расположена по отношению к координатной оси Y.
При одновременном соблюдении условий (54),(55),(56),(57) молекула считается влетевшей в исследуемый объем.
Условия (54),(55),(56),(57) должны выполняться и при розыгрыше вылета молекулы из объема.
1.3.3 Изменение структуры файла исходных данных
В связи с описанными дополнениями изменяется структура соответствующего файла исходных данных. Исходные данные для программы MODMD24 содержатся в файле inp24.txt и включают три дополнительных параметра - Tm, Tmm и ugol.
Исходные данные записываются в файле inp24.txt в следующем порядке:
n - количество испытаний;
r1 - внутренний радиус расчетного объема, отнесенный к радиусу исследуемого объема;
r2 - внешний радиус расчетного объема, отнесенный к радиусу исследуемого объема;
ugol [град.] - угол, определяющий размер отверстия в переднем (заднем) торце;
ас - коэффициент аккомодации энергии молекул;
V0 - направленная скорость молекул, влетающих в исследуемый объем;
Tm- абсолютная температура разреженного газа;
Tmm - температура направленного потока молекул, скорости которых распределены по закону Максвелла;
ugol_x (v) - угол между направлением вектора направленной скорости влетающих молекул и продольной осью МИП (Х);
ugol_у (v) - угол между проекцией вектора направленной скорости влетающих молекул на плоскость YZ, перпендикулярную продольной оси X, и осью Y [1];
m1 - метка, завершающая строку исходных данных.
Таким образом, строка исходных данных имеет следующий примерный вид:
10000 0.3 1.0 180 0.95 8000 300 0 0 90 1
1.4 Описание программы MODMD79
1.4.1 Разработка модернизированной математической модели
Основной идеей, положенной в основу модернизации программы MODMD24, является гипотеза о пропорциональности относительной концентрации в объеме сумме отношений времен пролета молекул в исследуемом объеме и свободном пространстве.
Рассмотрим два одинаковых объема (1 и 2), в которые попадают одинаковые потоки частиц. Разница в том, что в первом объеме рассеяние на стенках диффузное, а во втором - зеркальное. Известно, что при решении задач свободномолекулярной аэродинамики фиктивные стенки можно заменять на зеркально отражающие, поэтому концентрация в первом объеме будет равной концентрации окружающей среды. Таким образом, сравнивая между собой времена пролета молекул в первом и во втором объеме, можно сделать вывод об относительной концентрации в объеме с диффузноотражающими стенками.
Разобьем все влетающие частицы на классы: в класс j входят все частицы с определенной точкой влета, определенной скоростью влета, одинаково рассеиваемые в объеме 1. Пусть при зеркальном отражении в рассматриваемом объеме таких частиц будет Nj, тогда при диффузном отражении число частиц класса j в объеме 1 будет равно:
,
где tреалj- среднее время нахождения частиц класса j в объеме 1;
tидеалj - время нахождения частицы класса j в объеме 2.
Очевидно, что
Тогда, сложив частицы всех классов, получим, что в объеме 1 количество частиц равно:
,
где N - количество рассматриваемых частиц(число частиц в объеме 2).
При зеркальном отражении плотность внутри объема равна внешней плотности n = N/V, где V - величина рассматриваемого объема.
Рассмотрим относительную плотность в объеме 1.
(58)
Таким образом (58) - расчетная формула для всего объема.
Если объем 1 разбит на несколько частей, то отношение количества частиц класса j в частном объеме к количеству частиц класса j в объеме 1 равно tчастj/tреалj, где tчастj - время нахождения в частном объеме частиц класса j. Если Nj - количество частиц класса j в объеме 1, то количество частиц класса j в частном объеме
Очевидно, что
Просуммируем по всем классам j, тогда
Плотность в частном объеме естественно задавать как
,
тогда относительная плотность в частном объеме
.(59)
Таким образом, (59) - расчетная формула для частного объема.
Покажем, что (59) соответствует формуле (58). Пусть все частные объемы одинаковы по величине, тогда
По определению:
При равных частных объемах формула (60) очевидна.
Достоинством данного подхода является то, что для его реализации достаточно внести в текст программы MODMD24 лишь небольшие изменения.
1.4.2 Разработка программы моделирования
В соответствии с изложенными теоретическими идеями была создана программа MODMD79. Модификации в текст MODMD24 вносились таким образом, чтобы основное расчетное ядро оставалось без ощутимых изменений:
а) добавлены новые переменные:
- tid - время пролета молекулы в зеркально отражающей трубе;
- q1 - новое значение вспомогательного коэффициента q;
- Tsss - старое значение переменной Tss, теперь она вычисляется по другой формуле;
- NNN - значение относительной концентрации во всем объеме, рассчитанное по новому методу;
- nk1[] - массив, хранящий значения относительной концентрации в частных объемах, рассчитанные по новому методу, nk[] теперь вспомогательный массив, содержит периоды пролета текущей молекулы в частных объемах;
б) изменения в расчетных процедурах: в конце процедуры counttrack в Tss теперь накапливается величина taukon/tid, в Tsss сохранено старое значение Tss;
в) изменения в теле основной программы:
- значение q1 вычисляется по формуле: (r12 - r22)(1 - rA2)/(2 ks ps a);
- для каждой конкретной молекулы после моделирования влета определяется величина tid = |L/vx|, обнуляется массив nk[];
- после вылета текущей молекулы значения nk[] делятся на tid, полученный nk[] суммируется с nk1[];
- после отработки всех n траекторий вычисляется NNN по сумме всех элементов nk1[] и вспомогательным коэффициентам (a, q1);
- выходные данные содержат интегральные значения относительной концентрации, рассчитанные как по новому, так и по старому методу (распределение по частным объемам - только по новому методу).
2. ДОРАБОТКА МАТЕМАТИЧЕСКОЙ МОДЕЛИ ТОРЦЕВЫХ СТЕНОК
МИП имеет входные отверстия, изображенные на рис.1.8а и рис.1.8б. Отверстия в переднем и заднем торцах одинаковы по форме, размерам и, следовательно, площади. Программа расчета распределения концентрации в объеме МИП MODMD82 отличается от описанной в программе MODMD79 формулой для расчета площадей отверстий и механизмом задания координат точек влета молекул в исследуемый объем.
Рис. 1.8а Торцевые стенки с овальными отверстиями
Рис. 1.8б Торцевые стенки с круглыми отверстиями
Коэффициент прозрачности торцевых сеток, определяемый суммой площадей простых фигур: одной трети кольца с радиусами (r11+r12) и (r11 - r12) и 4 кругов с радиусом r12 для рис.1.8а и 4 кругов с радиусом r22 для рис.1.8б, вычисляется по следующей формуле:
(60)
Механизм задания координат точек влета в исследуемый объем связан с вероятностью влета. Для оценки вероятности попадания молекулы с координатами y2, z2 в объем, будем использовать параметр - радиус влета молекулы, равный расстоянию от оси X до точки влета на плоскости YZ:
(61)
Вероятность попадания молекулы с в объем зависит от значения . Радиус влета при заданной геометрии отверстий определяет отношение длины дуги с радиусом , проходящей в области отверстий, к длине всей окружности. Это отношение и соответствует вероятности попадания молекулы с радиусом влета в объем.
Из рис.4а видно, что не равная нулю вероятность попадания в объем существует только для случаев, когда:
r11 - r12 < < r11 + r12, r21 - r22 < < r21 + r22. (61)
Выражение, определяющее вероятность влета молекулы в область
(r11 - r12, r11 + r12) имеет вид:
(62)
где sum = 0,5 ( r11 + r12 + ).
Первое слагаемое в (62) соответствует отношению совокупной длины дуг в области отверстий к длине окружности с , без учета полукругов с радиусами r12. Второе слагаемое учитывает изменение длины дуг на полукругах в зависимости от .
Выражение, определяющее вероятность влета молекулы в область
(r21 - r22, r21 + r22) имеет вид:
(63)
где sum = 0.5 (r21 + r22 + ).
Розыгрыш точек влета происходит в несколько этапов:
а) В зависимости от значения компоненты скорости vx, определяется торец, через который молекула влетает в исследуемый объем, а, следовательно, координата x2.
б) При помощи генератора случайных чисел выбираются значения y2 и z2, а затем рассчитывается r0=rвл в диапазоне (0,1) до тех пор, пока полученное значение не удовлетворит условию (61).
в) При помощи генератора случайных чисел определяется y0 в диапазоне (0,1), причем условием попадания молекулы в объем является выполнение неравенства: y0 y, где y - рассчитывается либо по формуле (62), либо - (63), в зависимости какой торец используется.
При розыгрыше вылета молекулы из объема решается аналогичная задача, в которой выполняется сравнение y0 и y при известном (если выполняется условие (61)) и, в зависимости от результата, молекула либо вылетает из объема, либо отражается от торцевой стенки.
3. ИСПЫТАНИЯ И АНАЛИЗ ДАННЫХ
3.1 Цель испытаний
Целью испытаний является подтверждение правильности функционирования программ MODMD82 и MODMD82krug и последующий анализ газового потока при разгерметизации КА.
3.2 Краткие сведения о рабочей программе MODMD82 и MODMD82krug
молекула газ моделирование поток
Программы MODMD82 и MODMD82krug предназначены для моделирования аэродинамического взаимодействия набегающего потока с заданными параметрами и МИП, торцевые стенки которого могут изменять в зависимости от угла раскрытия, основным отличием программ MODMD82 и MODMD82krug от MODMD79 является наличие торцевых заслонок.
Исходными данными для моделирования являются:
а) Величины, описывающие геометрические размеры и основные свойства МИП: L, rA, r01, r02, ugol, kpr, T, ac.
б) Величины, описывающие набегающий поток: R, Tm,Tmm, V0,ugol_x ugol_y.
в) Величины, определяющие параметры моделирования: n,ps, ks, r1, r2.
Выходными величинами, полученными при моделировании, являются: распределение относительной концентрации разреженного газа по частным объемам (определяются значениями ps, ks, r1, r2) и ряд расчетных величин: Время, Конц., n1, n2, |Vx|, |Vвх|, |Vвых|, Tss.
Для нас большей интерес представляют данные с концентрацией, полученные в центральной части МИП, т.к. в ходе практических исследований было доказано, что рабочей поверхностью МИП является центральная часть внутренней поверхности датчика.
3.3 Общие положения
При проведении испытаний принимаются постоянными:
L = 2,25, rA = 0,276, r01 = 0,36, r02 = 0,9 (соответствуют геометрическим размерам МИП), T = 300 К, R = 590 м2/(с2 град), ps = 3, ks = 12, r1 = 0,276, r2 = 1,000, ugol_x = (0..180) град, ugol_y = (0..180) .
Остальные исходные данные имеют следующие начальные значения:
kpr = 1, ac = 0,9, Tm = 0 К, Tmm = 300, V0 = 0, n = 40000.
3.4 Формирование исходных данных
Исходные данные для моделирования формируются из данных, содержащихся в рабочем файле, и данных, которые считываются рабочей программой из файла исходных данных.
Информация в файле исходных данных записываются в одну или несколько строк, а в каждой строке в следующем порядке:
n - количество испытаний;
L - относительная длина МИП;
rA - относительный радиус анода;
r01 - максимальный относительный радиус входного отверстия;
r02 - минимальный относительный радиус входного отверстия;
ugol - угловая ширина входного отверстия, ;
kpr - коэффициент прозрачности входных отверстий;
r1 - верхний относительный радиус расчетного объема;
r2 - нижний относительный радиус расчетного объема;
V0 - средняя скорость набегающего потока, м/с;
ugol_x - угол между вектором средней скорости потока и осью X, ;
ugol_y - угол между проекцией вектора средней скорости потока на плоскость YOZ и осью Y, ;
R - удвоенное значение удельной газовой постоянной;
Tm - абсолютная температура газа при моделировании набегающего потока как равновесного газа, К;
Tmm - абсолютная температура потока 2-го вида (максвелловское распределение модуля скоростей, все молекулы движутся сонаправленно), К;
T - абсолютная температура стенок объема, К;
ac - коэффициент аккомодации;
m1 - метка, завершающая строку исходных данных.
Запись в файле исходных данных имеет следующий вид:
40000 2.25 0.276 0.31 0.965 90 1.0 0.276 1.0 8000 0 0 590 1800 0 300 0.9 1
......
......
40000 2.25 0.276 0.31 0.965 90 1.0 0.276 1.0 8000 90 0 590 1800 0 300 0.9 0
В случае m1=0 после цикла расчета с исходными данными, записанными в рассматриваемой строке, программа снова обращается к файлу исходных данных и считывает следующую строку исходных данных.
Расчет производится снова, но уже с другими исходными данными. Обновление исходных данных для моделирования проводится до ввода строки, в которой m1=1. Подобная запись исходных данных позволяет проводить целую серию расчетных экспериментов с различными исходными данными, не запуская каждый раз программу заново.
3.5 Методика испытаний
Исследование зависимости выходных величин, полученных при моделировании, от размера зон прозрачности включает проведение испытаний при следующих исходных данных:
Газовый поток 3го рода, который соответствует разгерметизации КА Датчик находится непосредственно перед течью в корпусе.
V0 = 0 м/с, Тmm = 300 K, kpr = 1, ugol_x = 0 град., ugol_y = 0 град
ugol = 1, 30, 60, 90, 120, 150, 180, 210, 240, 270, 300, 330, 360;
Количество испытаний n:
· 40000 для ugol = 30 ...360:
· 10000 для ugol = 1.
Результаты представлены в виде Графика 1 (MODMD82) и Графика 5 (MODMD82krug) (см.Приложение 3).
Газовый поток 3го рода при условии отклонения оси датчика (оси Х) от оси направления потока на 5 градусов
V0 = 0 м/с, Тmm = 300 K, kpr = 1, ugol_x = 5 град., ugol_y = 0 град
ugol = 1, 30, 60, 90, 120, 150, 180, 210, 240, 270, 300, 330, 360;
Количество испытаний n:
· 40000 для ugol = 30 ...360:
· 10000 для ugol = 1.
Результаты представлены в виде Графика 2 (MODMD82) и Графика 6 (MODMD82krug) (см.Приложение 3).
Газовый поток 3го рода при условии отклонения оси датчика (оси Х) от оси направления потока на 10 градусов
V0 = 0 м/с, Тmm = 300 K, kpr = 1, ugol_x = 10 град., ugol_y = 0 град
ugol = 1, 30, 60, 90, 120, 150, 180, 210, 240, 270, 300, 330, 360;
Количество испытаний n:
· 40000 для ugol = 30 ...360:
· 10000 для ugol = 1.
Результаты представлены в виде Графика 3 (MODMD82) и Графика 7 (MODMD82krug) (см.Приложение 3).
Газовый поток 3го рода при условии отклонения оси датчика (оси Х) от оси направления потока на 15 градусов
V0 = 0 м/с, Тmm = 300 K, kpr = 1, ugol_x = 10 град., ugol_y = 0 град
ugol = 1, 30, 60, 90, 120, 150, 180, 210, 240, 270, 300, 330, 360;
Количество испытаний n:
· 40000 для ugol = 30 ...360:
· 10000 для ugol = 1.
Результат представлен в виде Графика 4 (см.Приложение 3).
3.6 Интерпретация выходных данных
Концентрация разряженной газовой среды в каждой части является суммой концентраций, создаваемых в них каждым потоком. Благодаря неоднородной проницаемости сеток каждый поток будет создавать в каждой части камеры свою концентрацию разряженного газа. Концентрация будет больше в той части ионизационной камеры, которая для данного потока имеет большую проницаемость торцевой сетки. По результатам измерения токов частей катода можно судить о концентрации каждого потока.
Действительно, напряжение на выходе измерителей тока:
(64)
где и - коэффициенты преобразования верхней и нижний частей катода и подключенных к ним измерителей тока;
N1 и N2 - концентрация газа в нижней (Ndw) и верхней (Nup) частях ионизационных камеры, образованной частями катода.
При этом:
(65)
где Nc - концентрация равновесного газового потока;
Nv - концентрация направленного по оси динамического газового потока;
k1,l2,k2,l2 - коэффициенты пропорциональности концентрации в первой и второй частей ионизационный камеры для каждого вида газового потока соответственно.
Подставляем значения из (65) в (64), получаем:
(66)
Решение системы уравнений (3) имеет вид:
(67)
При симметричной конструкции ионизационной камеры и одинаковых характеристках измерителей тока можно принять ==, k1=k2=k, тогда:
(68)
где A,B,C - коэффициенты, определяемые конструтивными характеристиками устройства;
Uraz = U2 - U1 - разностный сигнал двуз измерителей тока;
Usum = U1 + U2 - суммарный сигнал двуз измерителей тока;
Коэффициенты A,B,C могут быть получены при градуировке устройства и расчетным путем. Для этого при отсутствии динамического газового потока Nv = 0, Uraz = 0 при известной концентрации равновесного газового потока Nc определяется коэффициент A. Затем при известной концентрации равновесного газового потока Nc и динамического газового потока Nv с заданными динамическими характеристиками определяются коэффициенты B и C.
Таким образом мы показали зависимость напряжения от концентрации внутри датчика.
Как можно предположить, то наиболее оптимальным вариантом является получение достаточного значения концентрации при обнаружении газового потока. Если датчик точно наведен на поток, то для получения точных значений
Как видно из построенных графиков, то чем меньше угол между осью Х и вектором потока, тем лучшие значения концентрации мы получаем.
Самым оптимальным случаем является, когда углы отклонения стремятся к нулю, тогда угол раскрытия торцов (при симметричном раскрытии) можно изменять в диапазоне от 90 до 190 градусов. В этом случаи N1 изменяется в пределах 10 процентов от максимального значения, а N2 остается практически неизменной, достигая своего максимум.
Отношение разностоноо сигнала к суммарному не дает точной картины, поэтому этот коэффициент нежелательно использовать в роли основного.
Моделирование при небольших отклонениях имтируют колебания КА, погрешности при фокусировке на поток, неточности в обратной связи и т.п.. Отклонение на угол в 5 градусов ведет к уменьшению максимума в пределах до 5 процентов и смещению оптимального диапозона раскрытия торцов. Для этого случая диапазон
ЗАКЛЮЧЕНИЕ
На текущей стадии были проанализированы несколько подходов к расчетам параметров среды и выявлен наиболее оптимальный.
Так же был разработан алгоритм моделирования, объединяющий локальные задачи в единую систему моделирования аэродинамического взаимодействия свободномолекулярного потока с объемом.
На следующем этапе планируется создать рабочую версию программы, основанной на языке программирования C.
СПИСОК ИСПОЛЬЗОВАННОЙ ЛИТЕРАТУРЫ
1. Г.Н. Паттерсон. Молекулярное течение газов. М. Государственное издательство физико-математической литературы. 1960г.
2. Газодинамика разреженных газов. Под редакцией М.Девиена. М. Издательство иностранной литературы. 1963г.
3. Ю.А. Кошмаров, Ю.А. Рыжов. Прикладная динамика разреженного газа. М. "Машиностроение". 1977г.
4. Берд Г. Молекулярная газовая динамика. - М.: Мир, 1981.
5. "Моделирования аэродинамического взаимодействия свободномолекулярных потоков с конструктивными элементами" (отчет об исследованиях), 1997 г.
ПРИЛОЖЕНИЕ
Текст программы modmd82
program dat(input,output);
label ll; const ks=12; ps=3; ProgID=82; r11=0.600;r12=0.240;s=3;
var a,q,vx,vy,vz,vyz,v,vT,vTm,vTmm,v0,w0,vv,y,y0,teta,xkpr,ksi,vn,ww0,q1,su1,sum1,su2,sum2:real;
x2,y2,z2,vx00,vyz00,vz00,vy00,Tm,Tmm,r0,r1,r2,ugol,ac,ugol_x,ugol_y,Tsr,Tss,Vsr,VVsr:real;
vxTmm,vxTmmv,vyzTmm,vyTmm,vyTmmv,vzTmm,vzTmmv,tid,Tsss,NNN, R,T,L,rA,r01,r02,kpr:real;
j,p,m1,k,kR,kL,f3,f4,sk:integer; n,i,n1,n2:longint; isp_No,m3,tekp:integer; defl:boolean;
inp1,out2,out3,out4:text; Fdir,conc,Nup,Ndn:real;
nk: array[1..2*ks,1..ps] of real;
nk1: array[1..2*ks,1..ps] of real; nc: array[1..ps,1..2] of real; d: array [1..ps] of real;
procedure Maxwell;
begin
repeat
v:=(random(5)+random);
y:=v*v*exp(1-v*v);
y0:=random;
until y0<=y;
if (v<=0.05) then v:=0.05;
end;
procedure scattering(Mx,My:real);
var vx0,vyz0,vy0,vz0,gamma,beta,alfa:real;
begin
repeat
beta:=pi/2*random;
y:=2*sin(beta)*cos(beta);
y0:=random;
until y0<=y;
Maxwell;
vv:=(1-ac)*vv*vv+ac*v*v+ac*(1-ac)*vv*vv*((v*v)/1.5-1);
if (vv<=0.0025) then vv:=0.05
else vv:=sqrt(vv);
alfa:=2*pi*random;
vx0:=vv*cos(beta); vyz0:=vv*sin(beta);
vy0:=vyz0*cos(alfa); vz0:=vyz0*sin(alfa);
if Mx<2 then begin vx:=Mx*vx0; vy:=vy0; vz:=Mx*vz0;end
else if not(y2*z2=0) then begin
if (y2*My)>0 then gamma:=arctan(z2/y2)+Pi-arctan(vy0/vx0)
else gamma:=arctan(z2/y2)-arctan(vy0/vx0);
vx:=-vz0; vy:=sqrt(vx0*vx0+vy0*vy0)*cos(gamma);
vz:=sqrt(vx0*vx0+vy0*vy0)*sin(gamma); end
else if y2=0 then begin
vx:=-vz0;vy:=-vy0*My*z2/abs(z2);vz:=-vx0*My*z2/abs(z2);end
else begin vx:=-vz0;vy:=-vx0*My*y2/abs(y2);vz:=vy0*My*y2/abs(y2);end;
end;
procedure counttrack;
label 2,3;
var x1,y1,z1,tmin1,tmax1,tmin2,tmax2,tmin,tmax,taumin,taumax,taukon,tAnod:real;
f0,f1,f2,fr1,fr2,fa2:real;
procedure countmolec;
var w: 0..ps; b1,b2,b3,b4,a1,a2:boolean;
d1,d2,d3,d4,tau1,tau2,tau3,tau4,tkon1,tkon2,t0,t00:real; label ii;
begin
Tsr:=Tsr+abs(tmax-tmin);
for w:=0 to ps-1 do
begin
tkon1:=((L*w/ps)-x1)/vx; tkon2:=((L*(w+1)/ps)-x1)/vx;
p:=w+1;
a1:=((tkon1>=tmin)and(tkon1<=tmax));
a2:=((tkon2>=tmin)and(tkon2<=tmax));
if a1 then if a2 then if vx>0 then begin t0:=tkon1;t00:=tkon2;end
else begin t0:=tkon2; t00:=tkon1;end
else if vx>0 then begin t0:=tkon1; t00:=tmax;end
else begin t0:=tmin; t00:=tkon1;end
else if a2 then if vx>0 then begin t0:=tmin; t00:=tkon2;end
else begin t0:=tkon2; t00:=tmax;end
else if (tkon1<tmin)and(tkon2>tmax) or (tkon1>tmax)and(tkon2<tmin)
then begin t0:=tmin; t00:=tmax;end
else goto ii;
d1:=y1-z1*vy/vz; d2:=d1;
tau1:=(d1-y1)/vy; tau2:=(d2-y1)/vy;
for kR:=1 to ks do
begin
d3:=vz*y1-vy*z1;
d4:=d3/(vz*cos(kR*Pi/ks)+vy*sin(kR*Pi/ks));
d3:=d3/(vz*cos(kR*Pi/ks)-vy*sin(kR*Pi/ks));
tau3:=(d3*cos(kR*Pi/ks)-y1)/vy; tau4:=(d4*cos(kR*Pi/ks)-y1)/vy;
b1:=((d1>0)and(tau1<=t00)and(tau1>=t0));
b2:=((d2>0)and(tau2<=t00)and(tau2>=t0));
b3:=((d3>0)and(tau3<=t00)and(tau3>=t0));
b4:=((d4>0)and(tau4<=t00)and(tau4>=t0));
if b1 then if b3 then nk[kR,p]:=nk[kR,p]+abs(tau1-tau3)
else if ((tau3-tau1)*d3)>0 then nk[kR,p]:=nk[kR,p]+abs(tau1-t00)
else nk[kR,p]:=nk[kR,p]+abs(tau1-t0)
else if b3 then if ((tau1-tau3)*d1)>0 then nk[kR,p]:=nk[kR,p]+abs(tau3-t00)
else nk[kR,p]:=nk[kR,p]+abs(tau3-t0)
else if (((d1>0)and(d3>0)and(((tau1>t00)and(tau3<t0))or((tau1<t0)and(tau3>t00))))
or ((d1*d3<0)and(((d3*(tau1-tau3)>0)and(tau1>t00)and(tau3>t00))
or((d1*(tau3-tau1)>0)and(tau1<t0)and(tau3<t0)))))
then nk[kR,p]:=nk[kR,p]+abs(t00-t0);
kL:=(2*ks+1)-kR;
if b2 then if b4 then nk[kL,p]:=nk[kL,p]+abs(tau2-tau4)
else if ((tau4-tau2)*d4)>0 then nk[kL,p]:=nk[kL,p]+abs(tau2-t00)
else nk[kL,p]:=nk[kL,p]+abs(tau2-t0)
else if b4 then if ((tau2-tau4)*d2)>0 then nk[kL,p]:=nk[kL,p]+abs(tau4-t00)
else nk[kL,p]:=nk[kL,p]+abs(tau4-t0)
else if (((d2>0)and(d4>0)and(((tau2>t00)and(tau4<t0))or((tau2<t0)and(tau4>t00))))
or ((d2*d4<0)and(((d4*(tau2-tau4)>0)and(tau2>t00)and(tau4>t00))
or((d2*(tau4-tau2)>0)and(tau2<t0)and(tau4<t0)))))
then nk[kL,p]:=nk[kL,p]+abs(t00-t0);
tau1:=tau3;tau2:=tau4;d1:=d3;d2:=d4;
end;
ii:end;
end;
begin
tmax:=0; tmin:=0;
x1:=x2;y1:=y2;z1:=z2;
f0:=vy*vy+vz*vz; f1:=(vy*y1+vz*z1)/f0;
f2:=f0-sqr(vz*y1-vy*z1);
fa2:=rA*rA*f0-sqr(vz*y1-vy*z1);
fr1:=r1*r1*f0-sqr(vz*y1-vy*z1);
fr2:=r2*r2*f0-sqr(vz*y1-vy*z1);
if f2>0 then
begin
f2:=sqrt(f2)/f0; taumin:=-f1-f2; taumax:=-f1+f2;
if not(vx=0) then if vx>0 then begin taukon:=(L-x1)/vx; f3:=-1;end
else begin taukon:=(0-x1)/vx; f3:=1;end
else taukon:=(2/vyz)+0.1;
tAnod:=-1;
if fr2>0 then
begin
fr2:=sqrt(fr2)/f0; tmin2:=-f1-fr2; tmax2:=-f1+fr2;
if fr1>0 then
begin
fr1:=sqrt(fr1)/f0; tmin1:=-f1-fr1; tmax1:=-f1+fr1;
if fa2>0 then begin fa2:=sqrt(fa2)/f0; tAnod:=-f1-fa2; end;
if (tmin2>0)and(tmin2<taukon) then Tmin:=tmin2
else if (tmin2<=0)and(tmin1>0) then Tmin:=0
else goto 2;
if (tmin1<taukon) then Tmax:=tmin1
else Tmax:=taukon;
countmolec;
2: if (tmax1>0)and(tmax1<taukon) then Tmin:=tmax1
else if (tmax1<=0)and(tmax2>0) then Tmin:=0
else goto 3;
if (tmax2<taukon) then Tmax:=tmax2
else Tmax:=taukon;
if (tAnod<0) then countmolec;
goto 3;
end;
if (tmin2>0)and(tmin2<taukon) then Tmin:=tmin2
else if (tmin2<=0)and(tmax2>0) then Tmin:=0
else goto 3;
if (tmax2>taukon) then Tmax:=taukon
else Tmax:=tmax2;
countmolec;
end;
3:end;
if taumax<taukon then begin taukon:=taumax; f3:=2;f4:=1;end;
if (tAnod>0)and(tAnod<taukon) then begin taukon:=tAnod; f3:=2;f4:=-1;end;
x2:=x1+taukon*vx; y2:=y1+taukon*vy; z2:=z1+taukon*vz;
Tss:=Tss+taukon/tid;Tsss:=Tsss+taukon; {ў Tss Є Ї«Ёў Ґ¬ taukon/tid; Tsss - бв ஥ § 票Ґ Tss}
if x2>L then x2:=L; if x2<0 then x2:=0;
end;
procedure byben;
var y0,y,sum:real;
begin
m3:=0;
if(r0>r11-r12)and(r0<r11+r12) then
begin
sum:=(r11+r12+r0)/2;
y0:=1/3+4/pi*2*arctan(sqrt((sum-r11)*(sum-r12)*(sum-r0)/sum)/(sum-r12));
y:=random;
if(y<=y0) then m3:=1;
end;
end;
{procedure vvod (var Nsu:mas);
begin
assign(out4,'c:\80Nsum.txt');
reset(out4);
for i1:=1 to s do
for j1:=1 to m do
read(out4,Nsu[i1,j1]);
{writeln(a[1,2]:7:4);
writeln(a[2,2]:7:4);
writeln(a[3,2]:7:4);
close(out4);
end;}
function summa:real;
var su1,sum1,su2,sum2,Koef:real;
begin
su1:=0;
for p:=1 to s do
su1:=su1+nc[p,1];
sum1:=su1/3;
su2:=0;
for p:=1 to s do
su2:=su2+nc[p,2];
sum2:=su2/3;
writeln(out2,'Nup=':8,'Ndn=':8);
writeln(out2,sum1:8:4,sum2:8:4);
Koef:=(sum2-sum1)/(sum1+sum2);
writeln(out2,'Koef=':8);
writeln(out2,Koef:8:4);
end;
begin
Randomize;
Assign(inp1,'c:\80inp.txt'); Reset(inp1);
writeln; writeln('Start, program ID=',ProgID:2); isp_No:=0;
repeat
Assign(out2,'c:\80all.txt'); Append(out2);
isp_No:=isp_No+1;
Read(inp1,n,L,rA,r01,r02,ugol,kpr,r1,r2,v0,ugol_x,ugol_y,R,Tm,Tmm,T,ac,m1);
writeln('---------------------------€б室лҐ-¤ лҐ------------------ь-ў-бҐаЁЁ:',isp_No:2);
writeln('n=':8,'L=':8,'rA=':8,'r01=':8,'r02=':8,'ugol=':8,'kpr=':8,'r1=':8,'r2=':8);
writeln(n:8,L:8:2,rA:8:3,r01:8:3,r02:8:3,ugol:8:1,kpr:8:3,r1:8:3,r2:8:3);
writeln('Vo=':8,'ug_x=':8,'ug_y=':8,'Rg=':8,'Tm=':8,'Tmm=':8,'T=':8,'ac=':8);
writeln(v0:8:0,ugol_x:8:1,ugol_y:8:1,R:8:1,Tm:8:1,Tmm:8:1,T:8:1,ac:8:3);
w0:=v0/sqrt(R*T); vTmm:=sqrt(Tmm/T); Fdir:=n;
vTm:=sqrt(Tm/T);vn:=0;Tsr:=0;Tss:=0;Vsr:=0;VVsr:=0;n1:=0;n2:=0;NNN:=0;
vx00:=w0*cos(ugol_x/180*pi);vyz00:=w0*sin(ugol_x/180*pi);
vy00:=vyz00*cos(ugol_y/180*pi);vz00:=vyz00*sin(ugol_y/180*pi);
vxTmm:=vTmm*cos(ugol_x/180*pi);vyzTmm:=vTmm*sin(ugol_x/180*pi)
vyTmm:=vyzTmm*cos(ugol_y/180*pi);vzTmm:=vyzTmm*sin(ugol_y/180*pi);
a:=(4*r12*r12 +(4/3)*r11*r12)*ugol/360;
q:=(r2*r2-r1*r1)*L/ps/(2*ks)/a;
q1:=(r2*r2-r1*r1)/ps/(2*ks)/a*(1-rA*rA); tekp:=0; {q1=(q/L)*(1-rA*rA), (1-rA^2) - Ї«®й ¤м ¬Ґ¦н«ҐЄвத®© з бвЁ в®а楢}
for k:=1 to 2*ks do for p:=1 to ps do nk1[k,p]:=0;
for i:=1 to n do
begin
if ((i/n*100) > (tekp+1)) then begin write(chr(8),chr(8),chr(8),chr(8),'%',tekp:3); tekp:=tekp+1 end;
Maxwell; defl:=false;
vxTmmv:=v*vxTmm;vyTmmv:=v*vyTmm;vzTmmv:=v*vzTmm;
Maxwell; v:=v*vTm+0.0000000001;
{teta:=pi*random;} teta:=2*random-1;
{vx:=v*cos(teta)+vx00+vxTmmv;} vx:=v*teta+vx00+vxTmmv;
if not(vx=0) then begin
if(vx>0) then begin sk:=-1;x2:=0; n1:=n1+1 end
else begin sk:=1; x2:=L; n2:=n2+1 end;
repeat y2:=2*Random-1; z2:=2*Random-1;
r0:=sqrt(y2*y2+z2*z2);
byben;
{xkpr:=Random;}
until (y2*y2+z2*z2<r02*r02)and(y2*y2+z2*z2>r01*r01)and(m3=1)and(((sk=1)and(-cos(pi*ugol/360)+y2/sqrt(y2*y2+z2*z2)>0))
or((sk=-1)and(cos(pi-pi*ugol/360)-y2/sqrt(y2*y2+z2*z2)>0)));
{vyz:=v*sin(teta);} vyz:=v*sqrt(1-teta*teta);
ksi:=2*pi*random;
vy:=vyz*cos(ksi)+vy00+vyTmmv; vz:=vyz*sin(ksi)+vz00+vzTmmv;
vv:=sqrt(vx*vx+vy*vy+vz*vz); VVsr:=VVsr+vv;
vn:=vn+abs(vx); j:=0; tid:=abs(L/vx); {<- (!) ¤«п Ї®«. ®вЄа. вагЎл, ®ЇҐа в®а ЇаЁбў. ¤«п tid в®«мЄ® §¤Ґбм!}
for k:=1 to 2*ks do for p:=1 to ps do nk[k,p]:=0;Tsss:=0;
repeat counttrack;
r0:=sqrt(y2*y2+z2*z2);
xkpr:=Random;
if (xkpr<kpr)and(abs(f3)=1)and(y2*y2+z2*z2>r01*r01)and(y2*y2+z2*z2<r02*r02)and(m3=1)and(((f3=-1)
and(-cos(pi*ugol/360)+y2/sqrt(y2*y2+z2*z2)>0))
or((f3=1)and(cos(pi-pi*ugol/360)-y2/sqrt(y2*y2+z2*z2)>0)))
then goto ll {Ґб«Ё Ї®Ї « ў §®г Їа®§а з®бвЁ, в® ўл«Ґв}
else begin scattering(f3,f4); defl:=true end {Ё зҐ, ®ва ¦ҐЁҐ}
until j=1;
ll: Vsr:=Vsr+vv; if (defl) then Fdir:=Fdir-1;
for k:=1 to 2*ks do for p:=1 to ps do nk[k,p]:=nk[k,p]/tid; { ®в®иҐЁҐ ўаҐ¬Ґ ॠ«.(¤«п Њ€Џ) Ё Ё¤Ґ «. (¤«п ®вЄалв®© вагЎл)}
for k:=1 to 2*ks do for p:=1 to ps do nk1[k,p]:=nk1[k,p]+nk[k,p];
end;
end;
for k:=1 to 2*ks do for p:=1 to ps do NNN:=NNN+nk1[k,p]; { c㬬Ёа㥬 ўбҐ ®в®иҐЁп ўаҐ¬Ґ ॠ«. Ё Ё¤Ґ «. Ї® з.®ЎкҐ¬ ¬}
ww0:=vn/n; Tss:=Tss*a/n*(1-rA*rA);NNN:=NNN/n/q1/ps/(2*ks); conc:=Tsr/n*ww0/q/ps/(2*ks); Fdir:=Fdir/n;
{conc - бв ஥ § з. ®в.Є®ж.ў Њ€Џ; ў Tss ¤® бЁе Ї®а Є Ї«Ёў «®бм ®в®иҐЁҐ taukon/tid, ⥯Ґам Tss ®в®а¬Ёа.; ®а¬Ёа㥬 NNN}
write(chr(8),chr(8),chr(8),chr(8));
writeln('------------------------------ђҐ§г«мв вл--------------------------------');
writeln('N_old=':8,'N_new=':8,'Tsr=':8,'n1=':8,'n2=':8,'|Vx|=':8,'|Vўе|=':8,'|Vўле|=':8,'Tss=':8,'Fdir=':7);
writeln(conc:8:4,NNN:8:4,Tsr/n:8:4,n1:8,n2:8,ww0:8:3,vvsr/n:8:3,vsr/n:8:3,Tss:8:4,Fdir:7:3);
writeln('------------------------------------------------------------------------');
writeln(out2,'ProgID=':8,'n=':8,'L=':8,'rA=':8,'r01=':8,'r02=':8,'ugol=':8,'kpr=':8,'r1=':8,'r2=':8);
writeln(out2,ProgID:8,n:8,L:8:2,rA:8:3,r01:8:3,r02:8:3,ugol:8:1,kpr:8:3,r1:8:3,r2:8:3);
writeln(out2,'V0=':8,'ug_x=':8,'ug_y=':8,'Rg=':8,'Tm=':8,'Tmm=':8,'T=':8,'ac=':8);
writeln(out2,v0:8:0,ugol_x:8:1,ugol_y:8:1,R:8:1,Tm:8:1,Tmm:8:1,T:8:1,ac:8:3);
writeln(out2,'N_old=':8,'N_new=':8,'Tsr=':8,'n1=':8,'n2=':8,'|Vx|=':8,'|Vin|=':8,'|Vout|=':8,'Tss=':8,'Fdir=':8);
writeln(out2,conc:8:4,NNN:8:4,Tsr/n:8:4,n1:8,n2:8,ww0:8:3,vvsr/n:8:3,vsr/n:8:3,Tss:8:4,Fdir:8:3);
writeln(out2,'Nup=':8,'Ndn=':8,'D=':8);
Nup:=0; Ndn:=0;
for p:=1 to ps do
begin
nc[p,1]:=0; nc[p,2]:=0; d[p]:=0;
for kR:=1 to 2*ks do
begin
if kR<=ks/2 then nc[p,1]:=nc[p,1]+nk1[kR,p];
if (kR>ks/2) and (kR<=1.5*ks) then nc[p,2]:=nc[p,2]+nk1[kR,p];
if kR>1.5*ks then nc[p,1]:=nc[p,1]+nk1[kR,p];
if kR<=ks then d[p]:=d[p]+nk1[kR,p];
if kR>ks then d[p]:=d[p]-nk1[kR,p];
end;
d[p]:=d[p]/n/q1/ks;
nc[p,1]:=nc[p,1]/n/q1/ks; nc[p,2]:=nc[p,2]/n/q1/ks;
write(out2,nc[p,1]:8:4,nc[p,2]:8:4,d[p]:8:4); writeln(out2);
Nup:=Nup+nc[p,1]; Ndn:=Ndn+nc[p,2];
end;
{ vvod(Nsu);}
Nup:=Nup/ps; Ndn:=Ndn/ps;
writeln(out2,'Conc':8,'priv':8,'volms':8);
for p:=1 to ps do
begin
write(out2,(nk1[1,p]+nk1[2*ks,p])/2/n/q1:8:3);
for kR:=2 to 2*ks do write(out2,(nk1[kR,p]+nk1[kR-1,p])/2/n/q1:8:3);
writeln(out2);
end;
summa;
writeln(out2);
close(out2);
Assign(out3,'c:\80int.txt'); Append(out3);
write(out3,ww0:9:4,vvsr/n:9:4,vsr/n:9:4,Tsr/n:9:4,Tss:9:4,Fdir:9:4,
conc:9:4,NNN:9:4,Ndn:9:4,Nup:9:4,Ndn-Nup:9:4);
write(out3,nc[2,2]-nc[2,1]:9:4); {changeble addon}
for p:=1 to 2 do
for kL:=1 to ps do write(out3,nc[kL,p]:9:4);
writeln(out3);
close(out3);
until m1=0;
writeln('Ready'); close(inp1);
end.
ПРИЛОЖЕНИЕ 2
Текст программы modmd82krug
program dat(input,output);
label ll; const ks=12; ps=3; ProgID=82; r21=0.600;r22=0.240;
var a,q,vx,vy,vz,vyz,v,vT,vTm,vTmm,v0,w0,vv,y,y0,teta,xkpr,ksi,vn,ww0,q1:real;
x2,y2,z2,vx00,vyz00,vz00,vy00,Tm,Tmm,r0,r1,r2,ugol,ac,ugol_x,ugol_y,Tsr,Tss,Vsr,VVsr:real;
vxTmm,vxTmmv,vyzTmm,vyTmm,vyTmmv,vzTmm,vzTmmv,tid,Tsss,NNN, R,T,L,rA,r01,r02,kpr:real;
j,p,m1,k,kR,kL,f3,f4,sk:integer; n,i,n1,n2:longint; isp_No,m3,tekp:integer; defl:boolean;
inp1,out2,out3:text; Fdir,conc,Nup,Ndn:real;
nk: array[1..2*ks,1..ps] of real;
nk1: array[1..2*ks,1..ps] of real; nc: array[1..ps,1..2] of real; d: array [1..ps] of real;
procedure Maxwell;
begin
repeat
v:=(random(5)+random);
y:=v*v*exp(1-v*v);
y0:=random;
until y0<=y;
if (v<=0.05) then v:=0.05;
end;
procedure scattering(Mx,My:real);
var vx0,vyz0,vy0,vz0,gamma,beta,alfa:real;
begin
repeat
beta:=pi/2*random;
y:=2*sin(beta)*cos(beta);
y0:=random;
until y0<=y;
Maxwell;
vv:=(1-ac)*vv*vv+ac*v*v+ac*(1-ac)*vv*vv*((v*v)/1.5-1);
if (vv<=0.0025) then vv:=0.05
else vv:=sqrt(vv);
alfa:=2*pi*random;
vx0:=vv*cos(beta); vyz0:=vv*sin(beta);
vy0:=vyz0*cos(alfa); vz0:=vyz0*sin(alfa);
if Mx<2 then begin vx:=Mx*vx0; vy:=vy0; vz:=Mx*vz0;end
else if not(y2*z2=0) then begin
if (y2*My)>0 then gamma:=arctan(z2/y2)+Pi-arctan(vy0/vx0)
else gamma:=arctan(z2/y2)-arctan(vy0/vx0);
vx:=-vz0; vy:=sqrt(vx0*vx0+vy0*vy0)*cos(gamma);
vz:=sqrt(vx0*vx0+vy0*vy0)*sin(gamma); end
else if y2=0 then begin
vx:=-vz0;vy:=-vy0*My*z2/abs(z2);vz:=-vx0*My*z2/abs(z2);end
else begin vx:=-vz0;vy:=-vx0*My*y2/abs(y2);vz:=vy0*My*y2/abs(y2);end;
end;
procedure counttrack;
label 2,3;
var x1,y1,z1,tmin1,tmax1,tmin2,tmax2,tmin,tmax,taumin,taumax,taukon,tAnod:real;
f0,f1,f2,fr1,fr2,fa2:real;
procedure countmolec;
var w: 0..ps; b1,b2,b3,b4,a1,a2:boolean;
d1,d2,d3,d4,tau1,tau2,tau3,tau4,tkon1,tkon2,t0,t00:real; label ii;
begin
Tsr:=Tsr+abs(tmax-tmin);
for w:=0 to ps-1 do
begin
tkon1:=((L*w/ps)-x1)/vx; tkon2:=((L*(w+1)/ps)-x1)/vx;
p:=w+1;
a1:=((tkon1>=tmin)and(tkon1<=tmax));
a2:=((tkon2>=tmin)and(tkon2<=tmax));
if a1 then if a2 then if vx>0 then begin t0:=tkon1;t00:=tkon2;end
else begin t0:=tkon2; t00:=tkon1;end
else if vx>0 then begin t0:=tkon1; t00:=tmax;end
else begin t0:=tmin; t00:=tkon1;end
else if a2 then if vx>0 then begin t0:=tmin; t00:=tkon2;end
else begin t0:=tkon2; t00:=tmax;end
else if (tkon1<tmin)and(tkon2>tmax) or (tkon1>tmax)and(tkon2<tmin)
then begin t0:=tmin; t00:=tmax;end
else goto ii;
d1:=y1-z1*vy/vz; d2:=d1;
tau1:=(d1-y1)/vy; tau2:=(d2-y1)/vy;
for kR:=1 to ks do
begin
d3:=vz*y1-vy*z1;
d4:=d3/(vz*cos(kR*Pi/ks)+vy*sin(kR*Pi/ks));
d3:=d3/(vz*cos(kR*Pi/ks)-vy*sin(kR*Pi/ks));
tau3:=(d3*cos(kR*Pi/ks)-y1)/vy; tau4:=(d4*cos(kR*Pi/ks)-y1)/vy;
b1:=((d1>0)and(tau1<=t00)and(tau1>=t0));
b2:=((d2>0)and(tau2<=t00)and(tau2>=t0));
b3:=((d3>0)and(tau3<=t00)and(tau3>=t0));
b4:=((d4>0)and(tau4<=t00)and(tau4>=t0));
if b1 then if b3 then nk[kR,p]:=nk[kR,p]+abs(tau1-tau3)
else if ((tau3-tau1)*d3)>0 then nk[kR,p]:=nk[kR,p]+abs(tau1-t00)
else nk[kR,p]:=nk[kR,p]+abs(tau1-t0)
else if b3 then if ((tau1-tau3)*d1)>0 then nk[kR,p]:=nk[kR,p]+abs(tau3-t00)
else nk[kR,p]:=nk[kR,p]+abs(tau3-t0)
else if (((d1>0)and(d3>0)and(((tau1>t00)and(tau3<t0))or((tau1<t0)and(tau3>t00))))
or ((d1*d3<0)and(((d3*(tau1-tau3)>0)and(tau1>t00)and(tau3>t00))
or((d1*(tau3-tau1)>0)and(tau1<t0)and(tau3<t0)))))
then nk[kR,p]:=nk[kR,p]+abs(t00-t0);
kL:=(2*ks+1)-kR;
if b2 then if b4 then nk[kL,p]:=nk[kL,p]+abs(tau2-tau4)
else if ((tau4-tau2)*d4)>0 then nk[kL,p]:=nk[kL,p]+abs(tau2-t00)
else nk[kL,p]:=nk[kL,p]+abs(tau2-t0)
else if b4 then if ((tau2-tau4)*d2)>0 then nk[kL,p]:=nk[kL,p]+abs(tau4-t00)
else nk[kL,p]:=nk[kL,p]+abs(tau4-t0)
else if (((d2>0)and(d4>0)and(((tau2>t00)and(tau4<t0))or((tau2<t0)and(tau4>t00))))
or ((d2*d4<0)and(((d4*(tau2-tau4)>0)and(tau2>t00)and(tau4>t00))
or((d2*(tau4-tau2)>0)and(tau2<t0)and(tau4<t0)))))
then nk[kL,p]:=nk[kL,p]+abs(t00-t0);
tau1:=tau3;tau2:=tau4;d1:=d3;d2:=d4;
end;
ii:end;
end;
begin
tmax:=0; tmin:=0;
x1:=x2;y1:=y2;z1:=z2;
f0:=vy*vy+vz*vz; f1:=(vy*y1+vz*z1)/f0;
f2:=f0-sqr(vz*y1-vy*z1);
fa2:=rA*rA*f0-sqr(vz*y1-vy*z1);
fr1:=r1*r1*f0-sqr(vz*y1-vy*z1);
fr2:=r2*r2*f0-sqr(vz*y1-vy*z1);
if f2>0 then
begin
f2:=sqrt(f2)/f0; taumin:=-f1-f2; taumax:=-f1+f2;
if not(vx=0) then if vx>0 then begin taukon:=(L-x1)/vx; f3:=-1;end
else begin taukon:=(0-x1)/vx; f3:=1;end
else taukon:=(2/vyz)+0.1;
tAnod:=-1;
if fr2>0 then
begin
fr2:=sqrt(fr2)/f0; tmin2:=-f1-fr2; tmax2:=-f1+fr2;
if fr1>0 then
begin
fr1:=sqrt(fr1)/f0; tmin1:=-f1-fr1; tmax1:=-f1+fr1;
if fa2>0 then begin fa2:=sqrt(fa2)/f0; tAnod:=-f1-fa2; end;
if (tmin2>0)and(tmin2<taukon) then Tmin:=tmin2
else if (tmin2<=0)and(tmin1>0) then Tmin:=0
else goto 2;
if (tmin1<taukon) then Tmax:=tmin1
else Tmax:=taukon;
countmolec;
2: if (tmax1>0)and(tmax1<taukon) then Tmin:=tmax1
else if (tmax1<=0)and(tmax2>0) then Tmin:=0
else goto 3;
if (tmax2<taukon) then Tmax:=tmax2
else Tmax:=taukon;
if (tAnod<0) then countmolec;
goto 3;
end;
if (tmin2>0)and(tmin2<taukon) then Tmin:=tmin2
else if (tmin2<=0)and(tmax2>0) then Tmin:=0
else goto 3;
if (tmax2>taukon) then Tmax:=taukon
else Tmax:=tmax2;
countmolec;
end;
3:end;
if taumax<taukon then begin taukon:=taumax; f3:=2;f4:=1;end;
if (tAnod>0)and(tAnod<taukon) then begin taukon:=tAnod; f3:=2;f4:=-1;end;
x2:=x1+taukon*vx; y2:=y1+taukon*vy; z2:=z1+taukon*vz;
Tss:=Tss+taukon/tid;Tsss:=Tsss+taukon; {ў Tss Є Ї«Ёў Ґ¬ taukon/tid; Tsss - бв ஥ § 票Ґ Tss}
if x2>L then x2:=L; if x2<0 then x2:=0;
end;
procedure byben;
var y0,y,sum:real;
begin
m3:=0;
if(r0>r21-r22)and(r0<r21+r22) then
begin
sum:=(r21+r22+r0)/2;
y0:=4/pi*2*arctan(sqrt((sum-r22)*(sum-r21)*(sum-r0)/sum)/(sum-r22));
y:=random;
if(y<=y0) then m3:=1;
end;
end;
begin
Randomize;
Assign(inp1,'c:\80inp.txt'); Reset(inp1);
writeln; writeln('Start, program ID=',ProgID:2); isp_No:=0;
repeat
Assign(out2,'c:\80all.txt'); Append(out2);
isp_No:=isp_No+1;
Read(inp1,n,L,rA,r01,r02,ugol,kpr,r1,r2,v0,ugol_x,ugol_y,R,Tm,Tmm,T,ac,m1);
writeln('---------------------------€б室лҐ-¤ лҐ------------------ь-ў-бҐаЁЁ:',isp_No:2);
writeln('n=':8,'L=':8,'rA=':8,'r01=':8,'r02=':8,'ugol=':8,'kpr=':8,'r1=':8,'r2=':8);
writeln(n:8,L:8:2,rA:8:3,r01:8:3,r02:8:3,ugol:8:1,kpr:8:3,r1:8:3,r2:8:3);
writeln('Vo=':8,'ug_x=':8,'ug_y=':8,'Rg=':8,'Tm=':8,'Tmm=':8,'T=':8,'ac=':8);
writeln(v0:8:0,ugol_x:8:1,ugol_y:8:1,R:8:1,Tm:8:1,Tmm:8:1,T:8:1,ac:8:3);
w0:=v0/sqrt(R*T); vTmm:=sqrt(Tmm/T); Fdir:=n;
vTm:=sqrt(Tm/T);vn:=0;Tsr:=0;Tss:=0;Vsr:=0;VVsr:=0;n1:=0;n2:=0;NNN:=0;
vx00:=w0*cos(ugol_x/180*pi);vyz00:=w0*sin(ugol_x/180*pi);
vy00:=vyz00*cos(ugol_y/180*pi);vz00:=vyz00*sin(ugol_y/180*pi);
vxTmm:=vTmm*cos(ugol_x/180*pi);vyzTmm:=vTmm*sin(ugol_x/180*pi);
vyTmm:=vyzTmm*cos(ugol_y/180*pi);vzTmm:=vyzTmm*sin(ugol_y/180*pi);
a:=(4*r22*r22)*ugol/360;
q:=(r2*r2-r1*r1)*L/ps/(2*ks)/a;
q1:=(r2*r2-r1*r1)/ps/(2*ks)/a*(1-rA*rA); tekp:=0; {q1=(q/L)*(1-rA*rA), (1-rA^2) - Ї«®й ¤м ¬Ґ¦н«ҐЄвத®© з бвЁ в®а楢}
for k:=1 to 2*ks do for p:=1 to ps do nk1[k,p]:=0;
for i:=1 to n do
begin
if ((i/n*100) > (tekp+1)) then begin write(chr(8),chr(8),chr(8),chr(8),'%',tekp:3); tekp:=tekp+1 end;
Maxwell; defl:=false;
vxTmmv:=v*vxTmm;vyTmmv:=v*vyTmm;vzTmmv:=v*vzTmm;
Maxwell; v:=v*vTm+0.0000000001;
{teta:=pi*random;} teta:=2*random-1;
{vx:=v*cos(teta)+vx00+vxTmmv;} vx:=v*teta+vx00+vxTmmv;
if not(vx=0) then begin
if(vx>0) then begin sk:=-1;x2:=0; n1:=n1+1 end
else begin sk:=1; x2:=L; n2:=n2+1 end;
repeat y2:=2*Random-1; z2:=2*Random-1;
r0:=sqrt(y2*y2+z2*z2);
byben;
{xkpr:=Random;}
until (y2*y2+z2*z2<r02*r02)and(y2*y2+z2*z2>r01*r01)and(m3=1);
{vyz:=v*sin(teta);} vyz:=v*sqrt(1-teta*teta);
ksi:=2*pi*random;
vy:=vyz*cos(ksi)+vy00+vyTmmv; vz:=vyz*sin(ksi)+vz00+vzTmmv;
vv:=sqrt(vx*vx+vy*vy+vz*vz); VVsr:=VVsr+vv;
vn:=vn+abs(vx); j:=0; tid:=abs(L/vx); {<- (!) ¤«п Ї®«. ®вЄа. вагЎл, ®ЇҐа в®а ЇаЁбў. ¤«п tid в®«мЄ® §¤Ґбм!}
for k:=1 to 2*ks do for p:=1 to ps do nk[k,p]:=0;Tsss:=0;
repeat counttrack;
r0:=sqrt(y2*y2+z2*z2);
xkpr:=Random;
if (xkpr<kpr)and(abs(f3)=1)and(y2*y2+z2*z2>r01*r01)and(y2*y2+z2*z2<r02*r02)and(m3=1)
then goto ll {Ґб«Ё Ї®Ї « ў §®г Їа®§а з®бвЁ, в® ўл«Ґв}
else begin scattering(f3,f4); defl:=true end {Ё зҐ, ®ва ¦ҐЁҐ}
until j=1;
ll: Vsr:=Vsr+vv; if (defl) then Fdir:=Fdir-1;
for k:=1 to 2*ks do for p:=1 to ps do nk[k,p]:=nk[k,p]/tid; { ®в®иҐЁҐ ўаҐ¬Ґ ॠ«.(¤«п Њ€Џ) Ё Ё¤Ґ «. (¤«п ®вЄалв®© вагЎл)}
for k:=1 to 2*ks do for p:=1 to ps do nk1[k,p]:=nk1[k,p]+nk[k,p];
end;
end;
for k:=1 to 2*ks do for p:=1 to ps do NNN:=NNN+nk1[k,p]; { c㬬Ёа㥬 ўбҐ ®в®иҐЁп ўаҐ¬Ґ ॠ«. Ё Ё¤Ґ «. Ї® з.®ЎкҐ¬ ¬}
ww0:=vn/n; Tss:=Tss*a/n*(1-rA*rA);NNN:=NNN/n/q1/ps/(2*ks); conc:=Tsr/n*ww0/q/ps/(2*ks); Fdir:=Fdir/n;
{conc - бв ஥ § з. ®в.Є®ж.ў Њ€Џ; ў Tss ¤® бЁе Ї®а Є Ї«Ёў «®бм ®в®иҐЁҐ taukon/tid, ⥯Ґам Tss ®в®а¬Ёа.; ®а¬Ёа㥬 NNN}
write(chr(8),chr(8),chr(8),chr(8));
writeln('------------------------------ђҐ§г«мв вл--------------------------------');
writeln('N_old=':8,'N_new=':8,'Tsr=':8,'n1=':8,'n2=':8,'|Vx|=':8,'|Vўе|=':8,'|Vўле|=':8,'Tss=':8,'Fdir=':7);
writeln(conc:8:4,NNN:8:4,Tsr/n:8:4,n1:8,n2:8,ww0:8:3,vvsr/n:8:3,vsr/n:8:3,Tss:8:4,Fdir:7:3);
writeln('------------------------------------------------------------------------');
writeln(out2,'ProgID=':8,'n=':8,'L=':8,'rA=':8,'r01=':8,'r02=':8,'ugol=':8,'kpr=':8,'r1=':8,'r2=':8);
writeln(out2,ProgID:8,n:8,L:8:2,rA:8:3,r01:8:3,r02:8:3,ugol:8:1,kpr:8:3,r1:8:3,r2:8:3);
writeln(out2,'V0=':8,'ug_x=':8,'ug_y=':8,'Rg=':8,'Tm=':8,'Tmm=':8,'T=':8,'ac=':8);
writeln(out2,v0:8:0,ugol_x:8:1,ugol_y:8:1,R:8:1,Tm:8:1,Tmm:8:1,T:8:1,ac:8:3);
writeln(out2,'N_old=':8,'N_new=':8,'Tsr=':8,'n1=':8,'n2=':8,'|Vx|=':8,'|Vin|=':8,'|Vout|=':8,'Tss=':8,'Fdir=':8);
writeln(out2,conc:8:4,NNN:8:4,Tsr/n:8:4,n1:8,n2:8,ww0:8:3,vvsr/n:8:3,vsr/n:8:3,Tss:8:4,Fdir:8:3);
writeln(out2,'Nup=':8,'Ndn=':8,'D=':8); Nup:=0; Ndn:=0;
for p:=1 to ps do
begin
nc[p,1]:=0; nc[p,2]:=0; d[p]:=0;
for kR:=1 to 2*ks do
begin
if kR<=ks/2 then nc[p,1]:=nc[p,1]+nk1[kR,p];
if (kR>ks/2) and (kR<=1.5*ks) then nc[p,2]:=nc[p,2]+nk1[kR,p];
if kR>1.5*ks then nc[p,1]:=nc[p,1]+nk1[kR,p];
if kR<=ks then d[p]:=d[p]+nk1[kR,p];
if kR>ks then d[p]:=d[p]-nk1[kR,p];
end;
d[p]:=d[p]/n/q1/ks;
nc[p,1]:=nc[p,1]/n/q1/ks; nc[p,2]:=nc[p,2]/n/q1/ks;
write(out2,nc[p,1]:8:4,nc[p,2]:8:4,d[p]:8:4); writeln(out2);
Nup:=Nup+nc[p,1]; Ndn:=Ndn+nc[p,2];
end;
Nup:=Nup/ps; Ndn:=Ndn/ps;
writeln(out2,'Conc':8,'priv':8,'volms':8);
for p:=1 to ps do
begin
write(out2,(nk1[1,p]+nk1[2*ks,p])/2/n/q1:8:3);
for kR:=2 to 2*ks do write(out2,(nk1[kR,p]+nk1[kR-1,p])/2/n/q1:8:3);
writeln(out2);
end;
writeln(out2);
close(out2);
Assign(out3,'c:\80int.txt'); Append(out3);
write(out3,ww0:9:4,vvsr/n:9:4,vsr/n:9:4,Tsr/n:9:4,Tss:9:4,Fdir:9:4,
conc:9:4,NNN:9:4,Ndn:9:4,Nup:9:4,Ndn-Nup:9:4);
write(out3,nc[2,2]-nc[2,1]:9:4); {changeble addon}
for p:=1 to 2 do
for kL:=1 to ps do write(out3,nc[kL,p]:9:4);
writeln(out3); close(out3);
until m1=0;
writeln('Ready'); close(inp1);
end.
ПРИЛОЖЕНИЕ 3
Графики полученные по результатам программы MODMD82
График 1. Распределение плотности при отклонении от оси Х на 5 градусов
График 3. Распределение плотности при отклонении от оси Х на 10 градусов
График 4. Распределение плотности при отклонении от оси Х на 15 градусов
ПРИЛОЖЕНИЕ 4
Графики полученные по результатам программы MODMD82krug
График 5. Распределение плотности при нахождении потока перед датчиком
График 6. Распределение плотности при нахождении потока перед датчиком
График 7. Распределение плотности при нахождении потока перед датчиком
Размещено на Allbest.ru
Подобные документы
Принцип действия магнитноразрядного измерителя плотности, механизм возникновения самостоятельного разряда. Разработка модернизированной математической модели моделирования аэродинамического взаимодействия набегающего потока с заданными параметрами.
дипломная работа [798,2 K], добавлен 03.02.2012Вычисление скорости молекул. Различия в скоростях молекул газа и жидкости. Экспериментальное определение скоростей молекул. Практические доказательства состоятельности молекулярно-кинетической теории строения вещества. Модуль скорости вращения.
презентация [336,7 K], добавлен 18.05.2011Определение структуры спектра атома, молекулы или образованной ими макросистемы их энергетическими уровнями. Спектры и структура атома водорода. Электронные состояния двухатомных молекул, электрические и оптические свойства. Молекулы с одинаковыми ядрами.
курсовая работа [52,0 K], добавлен 06.10.2009Определение центра тяжести молекулы и описание уравнения Шредингера для полной волновой функции молекулы. Расчет энергии молекулы и составление уравнения колебательной части молекулярной волновой функции. Движение электронов и молекулярная спектроскопия.
презентация [44,7 K], добавлен 19.02.2014Скорости газовых молекул. Понятие о распределении молекул газа по скоростям. Функция распределения Максвелла. Расчет среднеквадратичной скорости. Математическое определение вероятности. Распределение молекул идеального газа. Абсолютное значение скорости.
презентация [1,1 M], добавлен 13.02.2016Создание аппаратуры для измерения параметров разреженной атмосферы. Механизм возникновения самостоятельного газового разряда в скрещенных электрическом и магнитном полях. Алгоритм моделирования, разработка и описание программы. Испытания и анализ данных.
дипломная работа [1,4 M], добавлен 10.11.2011Молекулы идеального газа и скорости их движения. Упрyгoe стoлкнoвeниe мoлeкyлы сo стeнкoй. Опрeдeлeниe числа стoлкнoвeний мoлeкyл с плoщадкoй. Распрeдeлeниe мoлeкyл пo скoрoстям. Вывод формул для давления и энергии. Формула энергии идеального газа.
курсовая работа [48,6 K], добавлен 15.06.2009Гипотезы монополя Дирака. Магнитный заряд электрона, который тождественен кванту магнитного потока, наблюдаемого в условиях сверхпроводимости. Анализ эффекта квантования магнитного потока. Закон Кулона: взаимодействие электрического и магнитного заряда.
статья [205,4 K], добавлен 09.12.2010Кинематика точки. Способы задания движения. Определение понятия скорости точки и методы ее нахождения. Выявление ее значения при естественном способе задания равномерного движения. Способ графического представления скорости в декартовой системе координат.
презентация [2,3 M], добавлен 24.10.2013Сущность молекулы как наименьшей частицы вещества, обладающей всеми его химическими свойствами, экспериментальное доказательство их существования. Строение молекул, взаимосвязь атомов и их прочность. Методы измерения размеров молекул, их диаметра.
лабораторная работа [45,2 K], добавлен 11.02.2011
