Разработка и исследование торцевых поверхностей магнитноразрядного измерителя плотности

Выпускная работа бакалавра

тема: Разработка и анализ торцевых поверхностей магнитноразрядного измерителя плотности

Перечень условных обозначений и сокращений

МИП - магнитноразрядный измеритель плотности;

СВА - собственная внешняя атмосфера;

МС - магнитной системой

КА - космический аппарат

Введение

Работы по созданию аппаратуры для измерения параметров разреженной атмосферы в ЦНИИ РТК ведутся более 30 лет. Первый прибор "Альфа-1М" был создан в середине 70-х годов. За ним последовала разработка комплекса ДВЛС для изделия "Буран" и комплекса "Индикатор" для ОК "Мир", созданные по ТЗ НПО "Энергия". Основное назначение такой аппаратуры - измерение плотности разреженной атмосферы, как СВА, так и набегающего потока. Актуальность исследований, направленных на создание методов контроля герметичности орбитальных космических объектов путем разработки базовых решений по принципам построения, составу и аппаратной реализации средств контроля параметров разреженной атмосферы, основывается на опыте обеспечения жизнедеятельности долговременных орбитальных комплексов и особенно актуальна при возникновении нештатных ситуаций, подобных возникавшим на борту ОК “Мир”, в части обеспечения контроля негерметичности объектов в полете.

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

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

1.      
Анализ и исследование существующих математических моделей МИП

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

Как видно из рисунка 1.1 МИП представляет собой двухэлектродную цилиндрическую систему, состоящую из анода и катода. Между ними создается электрическое поле с напряженностью Е »105 В/м. Вдоль оси прибора магнитной системой создается постоянное магнитное поле напряженностью Н»5*104А/м.

Механизм возникновения самостоятельного разряда следующий. Электрон, появившийся в анодно-катодном промежутке, под действием электрического поля начинает двигаться к аноду по удлиненной из-за воздействия магнитного поля траектории. Из-за удлинения траектории повышается вероятность соударения электронов с молекулами разреженного газа и, следовательно, их ионизации. Обладая значительной кинетической энергией, положительные ионы при столкновении с катодом выбивают из его поверхности вторичные электроны, которые, двигаясь к аноду, также участвуют в ионизации газа. Ток положительных ионов на катод и ток электронов на анод в сумме численно равны величине тока в измерительной цепи. В результате ионизации газа возникает самостоятельный электрический разряд, ток которого в достаточно широком диапазоне имеет линейную зависимость от плотности газовой среды.

Диапазон давлений, при которых используются МИП, от 10-3 до 10-12 мм рт.ст. в значительной степени перекрывает требования по диапазону измерения на КА. Высокая эффективность преобразования, определяемая по степени ионизации разреженной газовой среды, которая в 100 раз превышает аналогичный показатель измерителей, использующих термоэлектронный метод ионизации, отсутствие накальных элементов или других внешних источников ионизации, отсутствие прецизионных элементов и деталей определяют высокие эксплуатационные характеристики МИП. Они устойчивы к внешним воздействиям, надежны, просты в эксплуатации, допускают длительную непрерывную работу. Благодаря этим характеристикам, и обеспечивается возможность их применения на КА.

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

В общем случае результаты измерений в натурных условиях не могут быть однозначно интерпретированы как плотность разреженного газа в зоне измерения МИП. Необходимо введение поправок, учитывающих отличие реальных условий проведения измерения от условий при наземной градуировке. К этим отличиям могут быть отнесены:

-      наличие направленных потоков разреженного газа, обусловленных как движением КА по орбите, так и собственным газовыделением систем ОК, причем указанные потоки являются свободномолекулярными, то есть частота столкновений молекул в объеме, соизмеримом с размерами МИП, пренебрежимо мала (для азота при T = 300 K и давлении 10-4 мм рт.ст. длина свободного пробега составляет 0,5 м при максимальном размере анодно-катодного промежутка МИП не более 0,05 м);

-    различием температур набегающего потока и МИП.

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

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

Первым шагом по этому пути является моделирование аэродинамического взаимодействия свободномолекулярных потоков нейтрального газа с конструктивными элементами МИП. Целью моделирования является определение распределения концентрации разреженного газа внутри МИП при различных параметрах набегающего потока. В качестве нейтрального газа при моделировании используется молекулярный азот.

1.1 Анализ программы MODMD05

Определение концентрации молекул разряженного газа в произвольном объеме.

Пусть на некоторый произвольный объем W воздействует свободномолекулярный поток, причем количество молекул, влетающих в объем в единицу времени, равен Kv (рис. 1.2).

Рис. 1.2 Произвольный объем со свободномолекулярным потоком

Необходимо определить количество Kw и среднюю концентрацию молекул Nw в объеме W.

При взаимодействии свободномолекулярного потока с объемом имеют место два следующих процесса:

-      увеличение числа молекул в объеме со скоростью Kv;

-      уменьшение числа молекул в объеме со скоростью x·Kw, где: x - коэффициент, показывающий какая часть молекул вылетает из объема в единицу времени.

Тогда можно записать:

,                                                                        (1)

где    t - время.

Решая это дифференциальное уравнение при начальных условиях: t=0, Kw=0, получим :

                                                                   (2)

В установившемся режиме, при t ® ¥:

                                                                                        (3)

Выясним физический смысл коэффициента x. Для этого рассмотрим решение уравнения (1) при отсутствии внешнего потока (Kv=0)

и начальных условиях t=0, Kw = Kw0. Тогда:

= Kw0 × exp(-x× t)                                                                     (4)

При t ® ¥: Kw ® 0, то есть все молекулы вылетают из объема W. Среднее время нахождения молекулы в объеме W может быть определено следующим образом:

                                                                (5)

Тогда выражения для определения Kw и Nw могут быть записаны:

= Kv × Tw и

                                                                     (6)

Внутри объема W распределение концентрации молекул (Nj) может быть определено следующим образом:

,                                                                                     (7)

где Wj - часть объема W, задаваемая параметром j (либо совокупностью параметров), которые определяют размеры объема Wj и его положение в объеме W (W= åWj);- среднее время нахождения молекул в объеме Wj (Tw= åTj).

Таким образом, средняя концентрация молекул (Nj) в объеме, задаваемом параметром j либо совокупностью параметров и входящем в объем W, определяется величиной объема (Wj), средним временем нахождения молекулы в объеме (Tj) и количеством молекул, влетающих в объем (W) в единицу времени.

Задача определения распределения молекул разреженного газа в объеме W может быть решена путем математического моделирования с использованием вероятностного численного метода (метода Монте-Карло).

Для этого необходимо:

-      задать параметры, определяющие представление объема W как совокупности объемов Wj;

-      смоделировать входной поток молекул разреженного газа, который в частном случае может быть представлен в виде направленного потока молекул разреженного газа, имеющих кроме того тепловую составляющую скорости, и определить величину Kv для заданных параметров потока и объема W;

-      смоделировать движение молекул внутри объема W, определить время нахождения каждой молекулы в составных частях объема Wj и рассчитать Tj - среднее время нахождения молекул в объеме Wj.

-      в соответствии с выражением (7) рассчитать распределение концентрации молекул (Nj) внутри объема W.

Моделирование объема

При моделировании объема W ограничимся рассмотрением круглого прямоугольного цилиндра, определяемого радиусом R0 и высотой L, у которого боковые стенки для молекул газа непроницаемы, а проницаемость основания может быть задана тем или иным образом. Размеры цилиндра при этом таковы, что потоки разреженного газа, воздействующие на цилиндр являются свободномолекулярными. В общем случае внутри этого цилиндра расположен второй цилиндр с радиусом RA < R0, непроницаемый для молекул газа и ограничивающий размеры объема W, в котором моделируется движение молекул. Этот цилиндр является анодом датчика. Такая модель объема W близка к реальной геометрии МИП.

Прямоугольную систему координат будем размещать в центре основания цилиндра так, чтобы ось X была направлена вдоль оси симметрии (рис.1.3).

Рис. 1.3 Модель МИП относительно осей координат

Объем W будем представлять в виде совокупности объемов Wpk, для которых:

-      1£ p £ps, где ps - количество равных частей (долей), на которые делится объем W вдоль оси X; при этом индексу 1 соответствует часть объема, включающая начало координат, нарастание индекса в обозначении соответствует увеличению величины x;

-      1 £ k £2 ks, где 2·ks - количество равных цилиндрических секторов, на которые объем W делится радиусами, проведенными из начала координат, в плоскости YZ; при этом отсчет секторов ведется по часовой стрелке вокруг положительного направления оси X (нумерация секторов ведется от положительного направления оси Y);

Для нахождения распределения концентрации разреженного газа в объеме W в направлении радиуса оснований цилиндра введем R1 и R2 - заданные радиусы (R0 > R2 > R1 > RA), ограничивающие объем, в котором рассчитывается концентрация молекул, частью цилиндрического сектора частный объем с основанием в виде кольцевого сектора с радиусами R1 и R2. Они представлены на рисунке 1.4.

Рис. 1.4 Вид боковых поверхностей с заданными радиусами

Частный объем рассчитывается по следующей формуле:

  ,                            (8)

где r1 = R1/R0, r2 = R2/R0, l = L/R0 - относительные геометрические размеры, отнесенные к радиусу R0.

Влет молекул в объем W и вылет из него происходят через торцевые стенки, которые в общем случае могут иметь зоны прозрачности и непрозрачности. Зоны прозрачности помимо геометрических размеров могут характеризоваться коэффициентом прозрачности (fs), причем 0 £ fs (y,z) £ 1. Принимаем, что непрозрачной части торцевой стенки соответствует fs(y,z) = 0, а полностью прозрачной части соответствует fs(y,z) = 1.

Кроме указанных геометрических характеристик объем W характеризуется абсолютной температурой стенок - Т, от которых происходит отражение молекул разряженного газа.

Моделирование набегающего потока

При моделировании набегающего потока определяется вектор скорости молекул, влетающих в объем W. Рассмотрим случай, когда набегающий поток может быть представлен в виде направленного потока молекул разреженного газа с заданной величиной скорости (V0) и углами относительно системы координат XYZ: jv (угол между вектором скорости V0 и осью X) и Yv (угол между проекцией вектора скорости V0 на плоскость YZ и осью Y), которые кроме того имеют случайную составляющую скорости, определяемую абсолютной температурой разреженного газа Tm. Тогда вектор скорости молекулы набегающего потока может быть представлен:

Модуль наиболее вероятной скорости (VT0) может быть определен следующим образом:

  ,                                                                            (10)

где    RT »590 м2/(с2×град)-соответствует молекулярному азоту;

jT и YT - углы, определяющие направление VT и представляющие собой случайные величины, распределенные по равномерному закону (0 £jT £ p и 0 £ YT <2. p)

В общем случае Vx может быть как больше, так и меньше 0. Первому случаю соответствует влет молекул в объем W через переднюю торцевую стенку (x = 0), а второму - через заднюю (x = l). При этом набегающий поток может быть рассмотрен как сумма двух аддитивных потоков, направленных навстречу друг другу.

Формирование набегающего потока в соответствии с выражением (9) позволяет смоделировать практически все виды реально существующих потоков нейтральных молекул. Действительно:

-      при V0 >> VT0 и jv ®0 (Vxv >> VxT) моделируется воздействие разреженной земной атмосферы для высот полета КА на МИП, ось X которого совпадает с направлением движения КА;

-      при V0 >> VT0 и jv =p/2 (Vxv = VxT) моделируется воздействие разреженной земной атмосферы для высот полета КА на МИП, ось X которого перпендикулярна направлению движения КА;

-      при V0=0 моделируется воздействие на МИП стационарной составляющей собственной внешней атмосферы КА;

-      при Tm=0 и V0¹0 моделируется воздействие на МИП потоков собственных газовыделений систем КА (возможно дополнительное моделирование параметров указанных потоков).

Моделирование движения молекулы внутри объема

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

Таким образом моделирование движения молекулы внутри объема включает:

-    расчет траектории движения молекулы;

-    анализ условий вылета молекулы;

-    формирование вектора скорости при отражении.

Исходными данными для расчета траектории движения молекулы являются координаты начала движения и проекции вектора скорости молекулы на оси координат X,Y,Z. Целью расчета является определение координат точки пересечения молекулой внутренней поверхности объема.

Анализ условий вылета молекулы заключается в сравнении координат точки пересечения молекулой внутренней поверхности объема с координатами зон прозрачности и непрозрачности и принятии решения либо о вылете молекулы из объема, либо об ее отражении от внутренней поверхности объема.

Для определения вектора скорости молекулы, отраженной от внутренней поверхности объема, необходимо смоделировать механизм ее взаимодействия с внутренней поверхностью.

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

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

Величиной, характеризующей тип взаимодействия и свойства поверхности (механические и химические) твердого тела, является коэффициент аккомодации энергии молекулы, определяемый выражением [1,с.215]:

 ,                                                                                   (11)

где    Ei - энергия падающей молекулы;- энергия отраженной молекулы;d - энергия, соответствующая температуре поверхности T.

Так как

 ,

то для модулей скоростей можно записать:

,                                                         (12)

где     Vr - модуль скорости отраженных молекул;модуль скорости падающих молекул;d - модуль скорости молекул, имеющих температуру равную температуре поверхности.

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

Величина коэффициента аккомодации ac для условий, сходных с условиями, при которых проводится моделирование, составляет от 0.5 до 1 [2, с.95].

Значения проекций вектора скорости отраженной молекулы определяются принятым для моделирования законом рассеяния молекулы от внутренней поверхности объема. Крайними случаями рассеяния являются зеркальное отражение и диффузное рассеяние, подобное отражению света от шероховатой белой поверхности [2, с.85].

В общем случае, сумму этих видов отражений описывает зеркально-диффузная функция распределения [3, п.1.7], [1, п.5.4]. Однако, из результатов экспериментов, представленных в [2, с.89] следует, что для полированных металлических поверхностей при всех углах падения довольно точно соблюдается рассеяние по диффузному закону, не учитывающему зеркальную составляющую отражения. При этом отраженные молекулы рассеиваются в пределе полусферы таким образом, что интенсивность потока отраженных молекул в телесном угле dw пропорциональна косинусу b между нормалью к поверхности, от которой происходит отражение (ось X0) и направлением рассеяния и не зависит от угла падения (рис.1.5).

Рис.1.5 Интенсивность потока отраженных молекул в телесном угле dw

Закон диффузного рассеяния (или закон косинуса), как функция распределения интенсивности потока отраженных молекул, был принят при моделировании отражения молекул от внутренних поверхностей объема.

Распределение концентрации молекул внутри объема

Количество молекул, влетающих в объем в единицу времени, может быть рассчитано следующим образом:

= Kv1 + Kv2 = n × Nv × Vxср1 × Sэкв1 + (1 -n) × Nv×Vxср2 × Sэкв2,

где     Kv1 и Kv2 - количество молекул, влетающих в объем через переднюю и заднюю стенки, соответственно;- концентрация набегающего потока;

£ n = n1/n £ 1 - коэффициент, равный отношению n1 - количества молекул, влетевших в объем через переднюю стенку (Vx > 0), к n - общему количеству молекул (общему количеству испытаний);ср1 =(1/n1) ×åVxi1 - средняя скорость влета молекул в объем W через переднюю стенку за n1 испытаний;ср2 = (1/n2)×å Vxi2- средняя скорость влета молекул в бъем W через заднюю стенку за n2= n - n1 испытаний.

,

где Sэкв1 - эквивалентная площадь влета молекул через переднюю стенку,

 ,

где Sэкв2 - эквивалентная площадь влета молекул через заднюю стенку.

Если Sэкв1 = Sэкв2 = Sэкв, то выражение (7) может быть записано:

Тогда относительное распределение концентрации молекул:

,                                                            (14)

где   

 

средняя скорость влета молекул в объем.

Из выражения (14) следует, что распределение концентрации молекул в объеме W не зависит от концентрации набегающего потока, разумеется, если сохраняются условия, при которых поток является свободномолекулярным.

Пусть размеры зон прозрачности: S1 - для передней стенки, S2 - для задней стенки и коэффициенты прозрачности, соответственно, fs1 и fs2 равны между собой. Тогда эквивалентная площадь влета:

экв= S × fs = q ×p×R02×fs = p × R02 × as,                                              (15)

где    S = S1 = S2, fs = fs1 = fs2;

q = S/ p× R02- относительный размер зоны прозрачности передней и задней стенок;= q × fs - интегральный коэффициент прозрачности для передней и задней стенок (при fs = 1 as = q).

Среднее время нахождения молекул в объеме Wpk может быть рассчитано по следующей формуле:

                                                                            (16)

где     dpk - расстояние, которое в i-ом испытании молекула пролетает в объеме Wpk;- скорость пролета объема Wpk в i-ом испытании;- число испытаний.

Подставляя в выражение (13) значения Vxср, Sэкв, получаем:

Kv = Nv Vxср p×r02 as                                                                       (17)

Подставляя в выражение (14) значения Vxср, Sэкв, Tpk, Wpk, получаем:

Или:

                                            (18)

где     Vxср0 = Vxср/Vст0 - средняя скорость влета молекул в объем W, отнормированная по Vст0;

 

суммарное время нахождения молекул в объеме за n испытаний, пронормированное по R0/Vст0;

 

модуль наиболее вероятной скорости молекулы, имеющей температуру стенок объема W.

Выражение (18) описывает основную расчетную формулу моделирования аэродинамического взаимодействия набегающего потока с МИП и определяет аналитическую зависимость между относительным распределением концентрации молекул разреженного газа в объеме W с одной стороны и исходными данными и результатами моделирования с другой.

Алгоритм моделирования.

На рисунке 1.6 представлен алгоритм объединяющий локальные задачи в единую систему моделирования аэродинамического взаимодействия свободномолекулярного потока с объемом.

Рис. 1.6 Алгоритм моделирования

Таким образом моделирование включает:

а)   формирование исходных данных. Исходными данными для моделирования являются: R, T, Tm, l, r1, r2, rA, V0, jv, Yv, n, ps, ks, as1 ,as2 ,ac.

б)       проведение n независимых испытаний. В каждом испытании в соответствии с исходными данными генерируется вектор скорости молекулы и координаты точки влета молекулы в объем W, рассчитывается траектория движения молекулы от точки влета до границы объема и определяются времена нахождения молекулы в объемах Wpk с учетом r1, r2 и координаты новой точки пересечения траектории движения молекулы с границей объема, а так же анализируются условия вылета молекулы из объема. Если условия вылета не выполняются, то моделируется отражение молекулы с учетом коэффициента аккомодации ac и вновь рассчитывается траектория ее движения до новой точки на границе объема; так повторяется до тех пор, пока условия вылета молекулы из объема W не будут выполнены, причем для каждой траектории определяются времена нахождения молекулы в объемах Wpk с учетом r1, r2 и координаты новой точки пересечения траектории движения молекулы с границей объема.

в)       по результатам n испытаний подсчитывается tpk - суммарное время нахождения n молекул в каждом из объемов Wpk, при расчете которого используются расстояния, отнормированные по величине R0 и скорости, отнормированные по величине Vст0 и Vxср0 - средняя скорость влета n молекул в объем W, отнормированная по Vст0.

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

Описание алгоритма моделирования.

Алгоритм, представленный на рисунке 1.6, программно реализуется при помощи трех файлов: рабочего файла, содержащего написанную на языке Pasccal программу моделирования, файла исходных данных и файла результатов. Моделируется аэродинамическое взаимодействие со свободномолекулярным потоком разреженного газа МИП, конструкция которого, близка к описанной в пукте 1.1.2. Эскиз МИП (дифференциальный МИП) представлен на рисунке 1.4.

Дифференциальный МИП имеет следующие размеры (в скобках приведены соответствующие относительные величины, отнесенные к R0):

-    радиус цилиндра: R0=14.5 мм;

-         радиус анода: RA=4 мм (ra=0,276);

-         длина цилиндра: L=32.6 мм (l=2,25);

-         радиус внутренней кромки отверстия: R01=4.5 мм (r01=0,31);

-         радиус внешней кромки отверстия: R02=14 мм (r02=0,965).

Формирование исходных данных.

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

Величины R, T, Tm, l, rA, ps, ks, а также r01, r02, необходимые для расчета as1, as2, вводятся как константы рабочей программы, остальные - через файл исходных данных. Информация в файле исходных данных записываются в одну или несколько строк, а в каждой строке в следующем порядке:- количество испытаний;- внутренний радиус расчетного объема, отнесенный к радиусу

исследуемого объема;- внешний радиус расчетного объема, отнесенный к радиусу

исследуемого объема;- направленная скорость молекул, влетающих в исследуемый объем;

ас - коэффициент аккомодации энергии молекул;_x (jv) - угол между направлением вектора направленной скорости влетающих молекул и осью Х;_у (Yv) - угол между проекцией вектора направленной скорости влетающих молекул на плоскость YZ и осью Y;- метка, завершающая строку исходных данных.

Запись в файле исходных данных имеет следующий вид:

0.3 1.0 8000 0.95 0 90 0

0.3 1.0 8000 0.8 0 90 1

В случае m1=0 после цикла расчета с исходными данными, записанными в рассматриваемой строке, программа снова обращается к файлу исходных данных и считывает следующую строку исходных данных.

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

Из величин, входящих в выражение (18) для расчета распределения относительной концентрации молекул (Npk0) внутри исследуемого объема, значения n, r1, r2, l, ps, ks являются исходными данными для моделирования, as рассчитывается до начала цикла испытаний на основании исходных данных о размерах и коэффициентах прозрачности торцевых стенок исследуемого объема, а величины Vxср0, tpk определяются в процессе моделирования.

Так как входные отверстия в переднем и заднем торцах одинаковы, то их коэффициенты прозрачности определяются следующим выражением:

.

Таким образом, выражение (18) может быть представлено в виде:

 .                                                                        (19)

Кроме as ещё ряд величин рассчитывается до начала испытаний.

Направленная и тепловая скорости молекул нормируются на скорость, соответствующую температуре стенок исследуемого объема:

 

                                              (20)

Проекции вектора направленной нормированной скорости W0 на координатные оси X,Y,Z с учетом углов влета, в соответствии с (9) определяются следующим образом:

,

,                          (21)

.

Генерация вектора скорости молекулы и координат точки влета

Основные принципы моделирования вектора скорости молекул набегающего потока изложены в пункте 1.1.3. Проекции вектора направленной скорости постоянны и одинаковы для всех молекул набегающего потока и определяются по формуле (21).

Модуль температурной составляющей скорости VT представляет собой случайную величину, распределенную по закону Максвелла.

Максвелловское распределение задается следующим образом:

а)   При помощи генератора случайных чисел задается значение вектора тепловой скорости VT в диапазоне от 0 до 5 (скорость VT нормирована относительно скорости молекулы, имеющей температуру стенок исследуемого объема).

В соответствии с законом Максвелла рассчитывается вероятность y того, что тепловая скорость молекулы будет иметь значение VT: y = VT2×exp (1-VT2).

б)   При помощи генератора случайных чисел выбирается произвольное значение вероятности y0 в диапазоне (0,1). Если y £ 0, то заданное значение вектора скорости VT используется для дальнейших расчетов. Если y0>y, процедура выбора VT повторяется.

Таким образом, скорости VT, вероятность появления которых в соответствии с распределением Максвелла мала (y<<1) будут встречаться в последовательности реже, чем скорости, для которых(y £ 1).

Проекции вектора температурной скорости VT на оси X,Y,Z определяются следующим образом:

С помощью генератора случайных чисел определяется величина c в диапазоне (-1,1), которая является косинусом угла между направлением вектора те-пловой скорости и положительным направлением оси Х. Выражение:

= VT×c                                                                                               (22)

определяет значение проекции вектора скорости на ось Х.

а)   Проекция вектора температурной скорости на плоскость YZ определяется выражением:

= VT sin(arccos(c)).

Учитывая тригонометрическую формулу sin2c + cos2c= 1, получаем:

                                                                        (23)

В плоскости YZ направление вектора скорости равновероятно и определяется случайным углом e , который откладывается от положительного направ-ления оси Y. Угол e выбирается при помощи генератора случайных чисел в диапазоне (0,2p). Выражения для проекций вектора скорости на оси Y и Z имеют вид:

vyT = vyzT × cos(e),= vyzT × sin(e)                                                                                 (24)

Описанный метод задания проекций вектора тепловой скорости обеспечивает его равновероятное распределение в выбранной системе координат.

Итоговые проекции вектора скорости на координатные оси X, Y, Z

в соответствии с (9) представляют собой суммы соответствующих проекций векторов направленной и тепловой скоростей:

vx = vxT + vx00,= vyT + vy00,                                                                             (25)

vz = vzT + vz00.

При каждом i-ом испытании необходимо задать точку влета молекулы в исследуемый объем, то есть координаты начала движения молекулы.

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

Торец, через который влетает i-ая молекула, определяет проекция вектора скорости vx. В случае vx>0 молекула влетает через отверстие в переднем торце (х1 = 0), в случае vx<0 - через отверстие в заднем торце (х1 = l). Для расчета Vxср0 - проекции на ось X средней скорости влета молекул в объем W, отнормированной по Vст0, рассчитывается сумма проекций на ось X скоростей влета молекул в объем W: vn:= vn + abs ( vx ).

Координаты y1, z1 точек, в которых молекулы влетают в исследуемый объем, должны быть равномерно распределены по площади входных отверстий. Учитывая это, координаты точки влета i - ой молекулы в исследуемый объем выбираются следующим образом:

а)   С помощью генератора случайных чисел определяется пара координат:

- координата z1 в диапазоне от -1 до 1;

координата y1 в диапазоне от -1 до 0 для (х1 = 0) или в диапазоне от 0 до 1 для (х1 =l).

б)   Анализируется выполнение условия нахождения выбранной пары координат внутри полукольца отверстия. :

> y12 + z12 > r012.                                                                            (26)

Это говорит о том, что выбранная пара находится внутри полукольца отверстия. Если условие (26) не выполняются, то генерируется новая пара чисел.

Таким образом, при большом количестве испытаний (n®¥) координаты точек влета молекул равномерно "засеивают" площади входных отверстий.

Расчет траектории движения молекулы включает в себя расчет времен пересечения молекулой границ расчетных объемов Wpk и определение координаты конечной точки прямолинейного участка траектории молекулы в исследуемом объеме W.

Для расчета требуемых времен и координат используются следующие начальные условия движения:

-     начальная точка прямолинейной траектории движения молекулы с координатами x1, y1, z1, которая может быть либо точкой влета молекулы в иссле-дуемый объем, либо конечной точкой предыдущей прямолинейной траектории движения молекулы (то есть точкой отражения);

-    проекции вектора скорости молекулы, задаваемые либо при влете молекулы в исследуемый объем, либо при отражении молекулы от стенок объема.

Прямолинейная траектория движения молекулы описывается уравнением движения:

.                                                                  (27)

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

Теперь рассмотрим процедуру поиска решения для случаев пересечения траектории движения молекулы с боковыми поверхностями цилиндров с радиусами rA, r1, r2, r0. Для этого рассмотрим движение молекулы в плоскости YZ.

Проекцией боковой поверхности цилиндра на плоскость YZ является окружность, и следовательно, в общем случае система уравнений может быть записана следующим образом:

Решая систему (28), получаем:

(y1+vyT)2+(z1+vzT)2=r2                                                                   (29)

Решив это уравнение относительно Т, получаем:

               (30)

Если значения tmax и tmin - комплексные числа, то траектория движения молекулы не пересекает рассматриваемую поверхность; положительные значения tmax и tmin означают, что молекула движется извне по направлению к рассматриваемой поверхности; отрицательные значения tmax и tmin означают, что молекула находится вне рассматриваемого цилиндра и удаляется от него; если tmax >0, а

tmin <0, то начальная точка движения находится внутри рассматриваемого цилиндра.

Времена пересечения траектории движения молекулы с поверхностями торцов цилиндра определяются по проекции траектории движения молекулы на ось Х.

Координатами торцов по оси X являются: х = 0 (для переднего торца) и

х = l (для заднего торца), тогда времена пересечения траектории движения молекулы с торцами t0и t1 определяются выражениями:

 

 .                                                      (32)

Для нахождения конечной точки прямолинейной траектории движения молекулы необходимо из времен tmax(r0), tmin(rA), t0 и t1 выбрать наименьшее положительное время -tкон. И тогда:

= vx × tкон + x1= vy × tкон + y1                                                                           (33)= vz × tкон + z1.

По результатам выбора tкон определяется поверхность, на которую попадает молекула. При этом задаются значения переменных Mх и My, определяющих эту поверхность: если молекула попадает на передний торец (tкон = t0), то Mх = 1, если на задний (tкон = t1), то Mх = -1 ; для боковой поверхности и анода - Mх = 2, причем при попадании молекулы на боковую стенку исследуемого объема (tкон =tmax(r0)) - My = 1, а при попадании молекулы на анод (tкон = tmin(rA)) - My = -1.

Для определения начального и конечного времени нахождения молекулы в расчетном объеме, ограниченном радиусами r1 и r2 и длиной l, необходимо рассчитать вещественные времена tmax(r2), tmin(r2), tmax(r1) и tmin(r1) и время tкон.

Значение Tmin- начального времени нахождения молекулы в расчетном объеме определяется местоположением начальной точки, знаком величин tmax(r1,r2), tmin(r1,r2), а также условием не превышения значениями tmax(r1,r2) или tmin(r1,r2) времени tкон. Так:

-     если r12<y12+z12<r22 , то Tmin=0;

-         если y12+ z12>r22 и 0< tmin(r2)< tкон, то Tmin=tmin(r2);

-     если tmin(rА) <0 или мнимое, 0<tmax(r1) <tкон, то Tmin=tmax(r1).

Значение Tmax - конечного времени нахождения молекулы в расчетном объеме - определяется из значений tкон , tmin(r1)< tкон, tmax(r2)< tкон. Так:

-     если значение tmin(rA)- мнимое, а tmax(r1,r2) - реальное и 0<tmin(r1)< tmax(r1)< tкон, то траектория молекулы пересекает расчетный объем дважды; в этом случае выби-раются две пары значений Tmin, Tmax;

-         если tmin(r2) мнимое или >tкон,то это означает, что молекула не попадает в расчетным объем на данном прямолинейном участке траектории.

Результатом описанной процедуры является одна или две пары значений времен влета Tmin и вылета Tmax для расчетного объема - в случае, если молекула в него попадает на данном прямолинейном участке траектории, а также координаты конечной точки прямолинейного участка движения x2, y2, z2 и значения

параметров Mx и My, определяющие поверхность, на которую попадает молекула.

Расчет времен пребывания молекулы в частных объемах Wpk

Расчет времен пребывания молекулы в частном объеме необходим для расчета относительного распределения концентраций молекул внутри цилиндра. Как показано в разделе 2, расчетный объем поделен на 2 × ks × ps частных объемов Wpk, времена нахождения молекулы в которых, есть разность времен влета и вылета.

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

Рассмотрим процедуру нахождения координат и времен пересечения траектории с образующими частный объем Wpk плоскостями. В проекции на плоскость YZ частный объем представляется в виде кольцевого сектора, внешняя и внутренняя дуги которого - части окружностей, соответствующие границам расчетного объема. Координаты пересечения траектории движения молекулы с радиальными границами кольцевого сектора рассчитываются из системы уравнений движения и образующей сектора, описываемой выражением:

×sinl - z×cosl = 0,                                                                              (34)

где    l - угол между осью Y и радиальной границей сектора.

Решением системы уравнений являются координаты:

Время пересечения траекторией боковой поверхности частного объема рассчитывается из выражения:

.                                                 (37)

Так как расчетный объем поделен в плоскости YZ на четное число секторов, то будем для секторов kR Î [1,ks] (правая половина проекции) отсчитывать углы a от 0 до p, а для секторов kL Î [ks+1,2ks] (левая половина) - углы l Î [0,-p] с интервалом p/ks. Тогда из выражения (37) получаем:

 ,      (38)

,

где    k0 Î [0,ks].

Зная значения tR и tL можно рассчитать времена пересечения траектории движения молекулы со всеми радиальными границами секторов. Времена пересечения траектории движения с плоскостями, разделяющими расчетный объем на ps равных долей по оси Х, рассчитываются из выражений:

 ,                                                                                  (39)

где    х - координата плоскостей, являющихся границами долей, на которые разбит расчетный объем вдоль оси Х.

Координата х принимает значения: (p0 × l)/ps; p0 Î[0,ps].

Определение значений времен нахождения молекулы в частном объеме производится следующим образом:

а)   Задается пара значений p0 и p0+1. Для них рассчитываются времена

t0¸1 и t0¸1+1 . Далее:

-         если хотя бы одно из них больше Tmin и меньше Tmax, или одно из них больше Tmax, а другое меньше Tmшт, то из времен t0¸1, t0¸1+1 ,Tmin, Tmax выбираются времена влета t0 и вылета t00 из долевого объема - объема, ограниченного поверхностями цилиндров с r1, r2 и плоскостями х'=(p0×l)/ps и х"=((p0+1) × l)/ps;

а)   если оба значения t0¸1 и t0¸1+1 меньше Tmin или больше Tmax, то берется следующая пара значений p0 и p0+1, и расчет повторяется.

б)       Для пары значений k0 и k0+1 рассчитываются времена tR и tR+1 , tL и tL+1 - как времена пересечения радиальных границ секторов kR и kL=2 × ks+1-kR. Из времен t0 и t00 и времен tR, tR+1 , tL и tL+1 выбираются времена влета и вылета для частных объемов в случае, когда хотя бы одно из времен tR, tR+1, tL и tL+1 лежит в пределах от t0 до t00, или же времена t0 и t00 лежат в пределах от tL,R до tL+1,R+1. Время нахождения молекулы в объеме Wpk, как разница времен влета и вылета, прибавляется к содержимому ячейки nk[k,p] и записывается в ячейку nk[k,p] матрицы времен. Расчет происходит в цикле от k0=0 до k0=ks-1 - когда будут рассчитаны времена tL,R ¸ tL+1,R+1 для секторов kR=ks и kL=ks+1.

в)       Задается следующая пара значений p0и p0+1, вплоть до значений p0=ps-1 и p0=ps.

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

Анализ условий вылета молекул из исследуемого объема

i-тое испытание заканчивается, когда молекула вылетает из исследуемого объема, то есть когда координата конечной точки прямолинейного участка траектории оказывается в области отверстия.

После расчета каждой конечной точки траектории при движении молекулы от стенки к стенке осуществляется проверка координаты y2 моле-кулы. При этом также анализируется переменная Mx, которой присвоены различные значения в зависимости от вида поверхности, на которую попала молекула. При Mx=1, что соответствует переднему торцу, при y2<0 и выполнении условия (26) молекула считается вылетевшей в отверстие на переднем торце; при Mx=-1, что соответствует заднему торцу, y2>0 и выполнении условия (26) молекула считается вылетевшей в отверстие на заднем торце.

После того, как i-тая молекула покидает исследуемый объем, начинается i+1ое испытание.

Формирование вектора скорости при отражении.

При столкновении молекулы, движущейся в исследуемом объеме, со стенкой происходит ее отражение. При этом меняется как модуль вектора скорости VV, так и направление ее движения, а следовательно и проекции скорости движения молекулы на оси X, Y, Z. Основные теоретические положения процесса отражения изложены в р.4.

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

Модуль вектора отраженной скорости в соответствии с (12) определяется из выражения:

                                             (40)

где    VT - величина скорости, соответствующая температуре стенки, выбранная из распределения Максвелла.

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

Для удобства моделирования распределения сформированной величины вектора скорости по проекциям вводится условная система координат X0Y0Z0, ориентированная в пространстве таким образом, что плоскость Y0Z0 - касательная к поверхности в точке отражения, а ось X0 направлена по нормали к этой поверхности. Началом координат в новой системе является точка отражения.

Проекции вектора скорости VV на оси координат X0Y0Z0 определяются выражениями:

vx0=VV×cosb;=VV×sinb;=vyz0×cosa;                                                                              (41)=vyz0×sina,

где     b - угол, который моделируется в соответствии с распределением вероятности P(b) по закону косинуса и задает отклонение вектора скорости от нормали к поверхности отражения в диапазоне [0,p/2];

a - угол между осью Y0 и проекцией вектора отраженной скорости на плоскость Y0Z0, выбираемый при помощи генератора случайных чисел в диапазоне [0,2p].

Интенсивность рассеяния молекул в направлении b пропорциональна cosb, следовательно:

 ,                                                                              (42)

где    db - телесный угол рассеяния;- площадь рассеяния.

В свою очередь:

dS = 2p × r0 × sinb × r0 × db = 2p × r02 × sinb × db,

тогда:

Полученную функцию отнормируем:

(b) = m × P0(b),     где òP0(b) = 1, тогда:

p×r02

Таким образом, нормированная функция распределения P0(b):

                                                             (43)

В соответствии с выведенным соотношением угол b выбирается с использованием генератора случайных чисел по методу аналогичному методу, примененному при моделировании распределения Максвелла (см.п.7.2.1.).

После определения проекций вектора скорости в системе координат X0Y0Z0 необходимо вернуться в исходную систему координат XYZ.

В случае отражения молекулы от торцов исследуемого объема эта процедура не представляет сложности, поскольку при отражении от переднего торца системы координат XYZ и X0Y0Z0 совпадают, и, следовательно:

vx = vx0,     vy = vy0,     vz = vz0.                                                     (44)

При отражении молекулы от заднего торца координатные оси Y и Y0 совпадают, а Х и Х0, Z и Z0 взаимно противоположны, таким образом:

vx = - vx0,   vy = vy0,   vz = - vz0.                                                  (45)

При отражении молекулы от боковой поверхности исследуемого объема и от анода для определения значения проекции вектора скорости на оси Y и Z, необходимо определить угол g между осью Y и проекцией вектора скорости отраженной молекулы на плоскость YZ. Для этого необходимо совместить начала координат систем XYZ и X0Y0Z0, то есть перенести точку отражения в центр координат плоскости YZ (рис.1.7). При отражении молекулы от верхней части боковой стенки исследуемого объема (y2>0) или нижней части анода (y2<0) угол g составит:

g = W1 + p - W3.

При отражении молекулы от боковой стенки при y2<0 или от анода

при y2>0:

g = W1 - W3

Причем:

,

Поскольку в случае отражения от поверхности анода или боковой стенки система координат X0Y0Z0 ориентирована таким образом, что координатные оси Z0 и X взаимно противоположны, а плоскости X0Y0 и YZ параллельны, то можно записать:

= -vz0,                                                               (46)

Тогда проекции вектора скорости vy и vz при отражении от анода или боковой стенки определяются выражениями:

vy =vyz × cosg,                                                                                  (47)= vyz × sing                                                               

В случае, если отражение происходит в точке, одна из координат которой (y2, z2) равна нулю, то выражение для проекций вектора скорости на оси YZ принимают вид:

а)  При отражении от анода:

-    если y2=0, то:          

,

,                                  (48)

-    если z2=0, то:

,

,                               (49)

б)  При отражении от боковой стенки:

-     если y2=0, то:

,

,                                (50)

-    
если z2=0, то:      

,

,                                  (51)

Определение положения начала системы координат X0Y0Z0 и, следовательно, выбор выражений для переноса проекций вектора скорости vx0, vy0, vz0 на оси X, Y, Z происходит при анализе переменных Mx и My, значения которых определяются при расчете прямолинейного участка траектории движения молекулы и соответствуют той поверхности исследуемого объема, от которой отражается молекула.

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

Расчет относительного распределения концентрации молекул в исследуемом объеме. Формирование матрицы результатов. Печать в файл результатов.

В течение цикла, состоящего из n испытаний, программа формирует матрицу nk[2 × ks,ps], элементы которой представляют собой суммы времен пребывания молекул (tpk) в каждом из частных объемов, на которые разбит расчетный объем, и сумму модулей проекций на ось X скорости влета молекул в объем W. После окончания испытаний, при i=n необходимо рассчитать распределение относительных концентраций молекул в рассматриваемых объемах и вывести результаты расчетов в удобном для последующей обработки виде.

Относительная концентрация молекул (Npk0) в частном объеме (Wpk) в соответствии с выражением (19) определяется: tpk - средним временем нахождения молекулы в частном объеме, Vxср0 - средней скоростью влета молекул в исследуемый объем (W) и величинами g и n, определяемыми на основании задаваемых исходных данных. Значения относительной концентрации молекул рассчитываются для каждого из 2 × ks × ps частных объемов.

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

в следующем порядке:- количество испытаний;

|Vx| - проекция на ось X средней скорости влета молекул в исследуемый объем;- внутренний радиус расчетного объема;- внешний радиус расчетного объема;- направленная скорость влета молекул;

ас - коэффициент аккомодации;_x - угол между направлением вектора направленной скорости

влетающих молекул и осью Х;_y - угол между проекцией вектора направленной скорости на плоскость YZ и осью Y;

|Vвх| - средняя скорость влета молекул в исследуемый объем;

|Vвых| - средняя скорость вылета молекул из исследуемого объма;ср. - среднее время пребывания молекул в расчетном объеме;ср. - средняя концентрация молекул в расчетном объеме;

В следующих строках представлены значения относительных концентраций молекул в частных объемах:

-    в нечетных строках - для правой половины расчетного объема

для секторов с номерами k Î [1,ks];

-    в четных строках - для левой половины расчетного объема для

секторов с номерами k Î [ks+1, 2 × ks].

Каждая пара строк значений концентраций соответствует одной издолей, на которые разбит частный объем вдоль оси Х.

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

Описание программы моделирования MODMD05

Программа modmd05.pas загружает в среду Паскаль файл исходных данных "inp2.txt" и файл результатов "concentr", устанавливает тип связи программы и файлов и закрывает их после отработки программы.

Все неоднократно повторяющиеся процедуры выделены в программе

в отдельные блоки - подпрограммы - для обеспечения большей компактности программы и повышения ее наглядности.

Формирование вектора скорости при отражении молекулы от поверхности осуществляется в подпрограмме scattering. Расчет времен пребывания молекулы в частных объемах, на которые разбит расчетный объем, выполняет подпрограмма countmolec. Расчет траектории движения молекулы в исследуемом объеме выполняет подпрограмма counttrack. Генерация точек, распределенных в соответствии с законом Максвелла, осуществляется в подпрограмме maxwell.

1.2 Описание программы MODMD24

На базе существующей программы моделирования взаимодействия молекул набегающего потока и СВА с конструктивными элементами MODMD05, была создана программа моделирования MODMD24, в которой введен дополнительный коэффициент, характеризующий исследуемый поток молекул: Tmm - температура направленного потока молекул с Максвелловским законом распределения модуля скоростей. При этом появляется возможность дополнительно моделировать набегающий поток, создаваемый собственными газовыделениями КА.

Также введены были дополнительные параметры, характеризующие входные отверстия датчика - центральный угол ugol, задающий размер части кольца, через которую осуществляется влет молекул в исследуемый объем и вылет из него, и коэффициент прозрачности отверстий kpr.

Моделирование потока собственных газовыделений.КА

В программе MODMD24 выражение для расчета модуля направленной скорости молекул, которое в MODMD05 описывается выражением (20), реализуется в виде:

,                                                         (52)

где     W0 - направленная скорость молекул, нормированная на модуль наиболее вероятной скорости молекулы, имеющей температуру стенок исследуемого объема: Vст02 = RT × T;- модуль направленной скорости молекул, влетающих в исследуемый объем;- температура стенок исследуемого объема;- значение газовой постоянной для молекулярного азота;- величина, распределенная по закону Максвелла;- температура направленного потока молекул с Максвелловским законом распределения модуля скорости.

При значении 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 - радиусы, определяющие размеры отверстий в переднем и заднем торцах МИП;[град.] - угол, определяющий размеры отверстий МИП (см.рис.1);- коэффициент прозрачности  отверстий МИП.

Рис.1.8 Модель датчика с учетом параметра ugol

Выбор точек влета-вылета молекул в исследуемый объем определяется геометрическим положением открытой области. Розыгрыш точек влета происходит в несколько этапов:

а)  В зависимости от знака компоненты скорости vx, определяется торец, через который молекула влетает в исследуемый объем, а, следовательно, координата x2.

б)   При помощи генератора случайных чисел выбирается случайное

число xkpr в диапазоне (0,1), которое характеризует вероятность того, что молекула, оказавшись в области отверстия, проникнет сквозь него в исследуемый объем.

в)   При помощи генератора случайных чисел выбираются значения координат y2 и z2 точки влета в диапазоне (-1,1) до тех пор, пока выбранная точка не окажется в области одного из отверстий. Для этого необходимо, чтобы совместно выполнялись условия:

+ z22 < r022,                                                                                      (54)

y22 + z22 > r012,                                                                                (55)

или

(при sk = 1)                                       (57) < kpr

разряд моделирование аэродинамический поток

где    sk - параметр, принимающий значение "-1" при розыгрыше влета в передний торец, и значение "1" при розыгрыше влета в задний торец.

Первое слагаемое в выражениях (56) задает косинус половины угла, определяющего размер входного отверстия (в радианах), второе -косинус угла, под которым предполагаемая точка влета молекулы расположена по отношению к координатной оси Y.

При одновременном соблюдении условий (54),(55),(56),(57) молекула считается влетевшей в исследуемый объем.

Условия (54),(55),(56),(57) должны выполняться и при розыгрыше вылета молекулы из объема.

Изменение структуры файла исходных данных

В связи с описанными дополнениями изменяется структура соответствующего файла исходных данных. Исходные данные для программы MODMD24 содержатся в файле inp24.txt и включают три дополнительных параметра - Tm, Tmm и ugol.

Исходные данные записываются в файле inp24.txt в следующем порядке:- количество испытаний;- внутренний радиус расчетного объема, отнесенный к радиусу исследуемого объема;- внешний радиус расчетного объема, отнесенный к радиусу исследуемого объема;[град.] - угол, определяющий размер отверстия в переднем (заднем) торце;

ас - коэффициент аккомодации энергии молекул;- направленная скорость молекул, влетающих в исследуемый объем;абсолютная температура разреженного газа;- температура направленного потока молекул, скорости которых распределены по закону Максвелла;_x (jv) - угол между направлением вектора направленной скорости влетающих молекул и продольной осью МИП (Х);_у (Yv) - угол между проекцией вектора направленной скорости влетающих молекул на плоскость YZ, перпендикулярную продольной оси X, и осью Y [1];- метка, завершающая строку исходных данных.

Таким образом, строка исходных данных имеет следующий примерный вид:

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 Общие положения

При проведении испытаний принимаются постоянными:

= 2,25, rA = 0,276, r01 = 0,36, r02 = 0,9

(соответствуют геометрическим размерам МИП),

= 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 - удвоенное значение удельной газовой постоянной;- абсолютная температура газа при моделировании набегающего потока как равновесного газа, К;- абсолютная температура потока 2-го вида (максвелловское распределение модуля скоростей, все молекулы движутся сонаправленно), К;- абсолютная температура стенок объема, К;- коэффициент аккомодации;- метка, завершающая строку исходных данных.

Запись в файле исходных данных имеет следующий вид:

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

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.

.         "Моделирования аэродинамического взаимодействия свободномолекулярных потоков с конструктивными элементами" (отчет об исследованиях), 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;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;;scattering (Mx,My:real); var vx0,vyz0,vy0,vz0,gamma,beta,alfa:real;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;;counttrack;2,3; var x1,y1,z1,tmin1,tmax1,tmin2,tmax2,tmin,tmax,taumin,taumax,taukon,tAnod:real; f0,f1,f2,fr1,fr2,fa2:real;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;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;:end;;:=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;

: 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;

: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;;byben;y0,y,sum:real;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;;

{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);(a[2,2]:7:4);(a[3,2]:7:4); close(out4); end;}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;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;(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)(-cos(pi*ugol/360)+y2/sqrt(y2*y2+z2*z2)>0))((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;: 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*('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);.

Приложение 2

Текст программы modmd82krug

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;Maxwell;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;;scattering(Mx,My:real); var vx0,vyz0,vy0,vz0,gamma,beta,alfa:real;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;;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;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;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;:end;;:=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;

: 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;

: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;;byben;y0,y,sum:real;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;;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;(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;: 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);.