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

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

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

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

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

МС — магнитной системой КА — космический аппарат

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

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

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

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

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

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

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

Диапазон давлений, при которых используются МИП, от 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;

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

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

(1)

где t — время.

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

(2)

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

(3)

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

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

Kw = Kw0 exp (- t) (4)

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

(5)

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

Kw = Kv Tw и

(6)

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

(7)

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

Tj — среднее время нахождения молекул в объеме 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, где 2ks — количество равных цилиндрических секторов, на которые объем 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: v (угол между вектором скорости V0 и осью X) и v (угол между проекцией вектора скорости V0 на плоскость YZ и осью Y), которые кроме того имеют случайную составляющую скорости, определяемую абсолютной температурой разреженного газа Tm. Тогда вектор скорости молекулы набегающего потока может быть представлен:

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

(10)

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

(11)

где Ei — энергия падающей молекулы;

Er — энергия отраженной молекулы;

E — энергия, соответствующая температуре поверхности T.

Так как

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

(12)

где Vr — модуль скорости отраженных молекул;

Viмодуль скорости падающих молекул;

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

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

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

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

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

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

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

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

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

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

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

Nv — концентрация набегающего потока;

0 = n1/n 1 — коэффициент, равный отношению n1 — количества молекул, влетевших в объем через переднюю стенку (Vx > 0), к n — общему количеству молекул (общему количеству испытаний);

Vxср1 =(1/n1) Vxi1 — средняя скорость влета молекул в объем W через переднюю стенку за n1 испытаний;

Vxср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экв= S fs = q R02fs = R02 as, (15)

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

q = S/ R02- относительный размер зоны прозрачности передней и задней стенок;

as = q fs — интегральный коэффициент прозрачности для передней и задней стенок (при fs = 1 as = q).

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

(16)

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

Vpk — скорость пролета объема Wpk в i-ом испытании;

n — число испытаний.

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

Kv = Nv Vxср 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, v, v, 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, вводятся как константы рабочей программы, остальные — через файл исходных данных. Информация в файле исходных данных записываются в одну или несколько строк, а в каждой строке в следующем порядке:

n — количество испытаний;

r1 — внутренний радиус расчетного объема, отнесенный к радиусу

исследуемого объема;

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

исследуемого объема;

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

ас — коэффициент аккомодации энергии молекул;

ugol_x (v) — угол между направлением вектора направленной скорости влетающих молекул и осью Х;

ugol_у (v) — угол между проекцией вектора направленной скорости влетающих молекул на плоскость YZ и осью Y;

m1 — метка, завершающая строку исходных данных.

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

10 000 0.3 1.0 8000 0.95 0 90 0

10 000 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 = VT2exp (1-VT2).

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

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

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

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

vxT = VT (22)

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

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

vyzT = VT sin (arccos ()).

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

(23)

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

vyT = vyzT cos (),

vzT = vyzT sin () (24)

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

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

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

vx = vxT + vx00,

vy = 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).

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

r022 > 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)

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

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

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

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

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

. (32)

Для нахождения конечной точки прямолинейной траектории движения молекулы необходимо из времен max (r0), min (rA), 0 и 1 выбрать наименьшее положительное времякон. И тогда:

x2 = vx кон + x1

y2 = vy кон + y1 (33)

z2 = vz кон + z1.

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

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

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

— если r12

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

— если min (rА) <0 или мнимое, 0<�кон, то Tmin=max (r1).

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

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

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

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

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

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

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

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

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

ysin — zcos = 0, (34)

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

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

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

. (37)

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

(38)

где k0 [0,ks].

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

(39)

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

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

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

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

и +1. Далее:

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

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

б) Для пары значений k0 и k0+1 рассчитываются времена R и R+1, L и L+1 — как времена пересечения радиальных границ секторов kR и kL=2 ks+1-kR. Из времен t0 и t00 и времен R, R+1, L и L+1 выбираются времена влета и вылета для частных объемов в случае, когда хотя бы одно из времен R, R+1, L и L+1 лежит в пределах от t0 до t00, или же времена t0 и t00 лежат в пределах от L, R до L+1,R+1. Время нахождения молекулы в объеме Wpk, как разница времен влета и вылета, прибавляется к содержимому ячейки nk[k, p] и записывается в ячейку nk[k, p] матрицы времен. Расчет происходит в цикле от k0=0 до k0=ks-1 — когда будут рассчитаны времена L, R L+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=VVcos;

vyz0=VVsin;

vy0=vyz0cos; (41)

vz0=vyz0sin,

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

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

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

(42)

где d — телесный угол рассеяния;

dS — площадь рассеяния.

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

dS = 2 r0 sin r0 d = 2 r02 sin d,

тогда:

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

P () = m P0(), где P0() = 1, тогда:

r02

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

(43)

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

= 1 + - 3.

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

при y2>0:

= 1 — 3

Причем:

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

vx = -vz0, (46)

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

vy =vyz cos, (47)

vz = vyz sin

В случае, если отражение происходит в точке, одна из координат которой (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 частных объемов.

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

в следующем порядке:

n — количество испытаний;

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

r1 — внутренний радиус расчетного объема;

r2 — внешний радиус расчетного объема;

V0 — направленная скорость влета молекул;

ас — коэффициент аккомодации;

ugol_x — угол между направлением вектора направленной скорости

влетающих молекул и осью Х;

ugol_y — угол между проекцией вектора направленной скорости на плоскость YZ и осью Y;

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

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

Tср. — среднее время пребывания молекул в расчетном объеме;

Nср. — средняя концентрация молекул в расчетном объеме;

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

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

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

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

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

Каждая пара строк значений концентраций соответствует одной из

ps долей, на которые разбит частный объем вдоль оси Х.

В заключительной строке выводятся значения фактического количества испытаний 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;