Главная
База знаний "Allbest"
Физика и энергетика
Разработка и анализ торцевых поверхностей магнитноразрядного измерителя плотности
Разработка и анализ торцевых поверхностей магнитноразрядного измерителя плотности
Принцип действия магнитноразрядного измерителя плотности, механизм возникновения самостоятельного разряда. Разработка модернизированной математической модели моделирования аэродинамического взаимодействия набегающего потока с заданными параметрами.
| Рубрика | Физика и энергетика |
| Вид | дипломная работа |
| Язык | русский |
| Дата добавления | 03.02.2012 |
| Размер файла | 798,2 K |
Отправить свою хорошую работу в базу знаний просто. Используйте форму, расположенную ниже
Студенты, аспиранты, молодые ученые, использующие базу знаний в своей учебе и работе, будут вам очень благодарны.
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 моделируется воздействие разреженной земной атмосферы для высот полета КА на МИП.
Моделирование влета (вылета) молекул в датчик
Программа для моделирования распределения концентрации аэродинамического потока молекул в датчике 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) (57)
xkpr < kpr
разряд моделирование аэродинамический поток
где sk - параметр, принимающий значение "-1" при розыгрыше влета в передний торец, и значение "1" при розыгрыше влета в задний торец.
Первое слагаемое в выражениях (56) задает косинус половины угла, определяющего размер входного отверстия (в радианах), второе -косинус угла, под которым предполагаемая точка влета молекулы расположена по отношению к координатной оси Y.
При одновременном соблюдении условий (54),(55),(56),(57) молекула считается влетевшей в исследуемый объем.
Условия (54),(55),(56),(57) должны выполняться и при розыгрыше вылета молекулы из объема.
Изменение структуры файла исходных данных
В связи с описанными дополнениями изменяется структура соответствующего файла исходных данных. Исходные данные для программы 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.3 Описание программы MODMD79
Разработка модернизированной математической модели
Основной идеей, положенной в основу модернизации программы 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 лишь небольшие изменения.
Разработка программы моделирования
В соответствии с изложенными теоретическими идеями была создана программа 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.90
В случае 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го рода при условии отклонения оси датчика (оси Х) от оси направления потока на 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 г.
Приложение 1
Текст программы 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); writelnwriteln('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*
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); writelnwriteln('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); 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)); writelnwriteln('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); writelnwriteln(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.
Размещено на Allbest.ru
Подобные документы
Определение концентрации молекул разряженного газа в произвольном объеме. Моделирование набегающего потока, движения молекулы внутри объема. Генерация вектора скорости молекулы и координат точки влета. Моделирование потока собственных газовыделений.
дипломная работа [1,8 M], добавлен 06.07.2011Создание аппаратуры для измерения параметров разреженной атмосферы. Механизм возникновения самостоятельного газового разряда в скрещенных электрическом и магнитном полях. Алгоритм моделирования, разработка и описание программы. Испытания и анализ данных.
дипломная работа [1,4 M], добавлен 10.11.2011История и основное энергетическое понятие фотометрии; визуальные и физические методы. Разработка оптико-механической схемы лазерного измерителя скорости на основе спекл-полей; расчет оптических параметров, чувствительности; описание установки в динамике.
курсовая работа [123,9 K], добавлен 19.05.2013Определение линейных, фазных токов, размеров и витков обмоток. Среднее значение плотности тока в обмотках. Расчет обмотки и площади поверхностей охлаждения обмоток. Определение плотности теплового потока. Расчет стоимости трансформатора и электрозатрат.
курсовая работа [1,8 M], добавлен 23.01.2011Условия возникновения электрического разряда в газах. Принцип ионизации газов. Механизм электропроводности газов. Несамостоятельный газовый разряд. Самостоятельный газовый разряд. Различные типы самостоятельного разряда и их техническое применние.
реферат [32,3 K], добавлен 21.05.2008Изучение методики обработки результатов измерений. Определение плотности металлической пластинки с заданной массой вещества. Расчет относительной и абсолютной погрешности определения плотности материала. Методика расчета погрешности вычислений плотности.
лабораторная работа [102,4 K], добавлен 24.10.2022Решение экспериментальных задач по определению плотности твердых веществ и растворов, с различной массовой долей растворенного вещества. Измерение плотности веществ, оценка границ погрешностей. Установление зависимости плотности растворов от концентрации.
курсовая работа [922,0 K], добавлен 17.01.2014Анализ основных форм самостоятельного разряда в газе. Исследование влияния относительной плотности воздуха на электрическую прочность разрядного промежутка. Определение значения расстояния между электродами, радиуса их кривизны для электрического поля.
лабораторная работа [164,5 K], добавлен 07.02.2015Исходные данные и расчетные формулы для определения плотности твердых тел правильной формы. Средства измерений, их характеристики. Оценка границы относительной, абсолютной погрешностей результата измерения плотности по причине неровности поверхности тела.
лабораторная работа [26,9 K], добавлен 30.12.2010Причины возникновения подъемной силы летательного аппарата. Заслуги Жуковского в развитии аэродинамики. Понятие турбулентности и процесс возникновения зоны повышенной плотности на передней части снаряда. Принципы всасывания потока воздуха в двигатель.
реферат [2,2 M], добавлен 01.06.2013
