II.3 Метод Монте-Карло пробной частицы для свободномолекулярного режима (Нестреров С.Б., Васильев Ю.К., Андросов А.В. Методы расчета вакуумных систем), страница 3

//Нахождение значения второго параметра t при пересечении
//конуса, заданного уравнением
//(z-c)²=ctg²(x²+y²)
//Аргументы: X, Y, Z– координаты точки вылета
//(x₀, y₀, z₀);
// L, M, N – направляющие косинусы l, //m, n;
//C – соответствующий член уравнения конуса c; G – угол 
//в уравнении конуса.
function GetConusTP(var X,Y,Z,L,M,N,C,G:Extended): Extended;
begin
if (((X*L+Y*M-Tan(G)*Tan(G)*N*(Z-C))*(X*L+Y*M-Tan(G)*Tan(G)*N*(Z-C))-(L*L+M*M-Tan(G)*Tan(G)*N*N)*(X*X+Y*Y-Tan(G)*Tan(G)*(Z-C)*(Z-C)))<0)
then GetConusTP:=0
else GetConusTP:=(-(X*L+Y*M-Tan(G)*Tan(G)*N*(Z-C))+sqrt(((X*L+Y*M-Tan(G)*Tan(G)*N*(Z-C))*(X*L+Y*M-Tan(G)*Tan(G)*N*(Z-C))-(L*L+M*M-Tan(G)*Tan(G)*N*N)*(X*X+Y*Y-Tan(G)*Tan(G)*(Z-C)*(Z-C)))))/(L*L+M*M-Tan(G)*Tan(G)*N*N);
end;

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

1. Поиск параметра, соответствующего точке столкновения

Далее, после формирования массива параметров t для всех поверхностей системы, необходимо найти параметр, соответствующий точке столкновения. Для этого нужно найти минимальный по модулю параметр t, причем значения x, y, z, найденные по формулам:

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

Поиск минимального, удовлетворяющего всем вышеописанным условиям, значения ведется следующим образом:

1. Задается некоторое начальное значение t

Например t =1^-10 м, т. е. радиус молекулы в метрах. (cм. также п. I «инициализация служебных параметров алгоритма»).

1. Ищется параметр t, значение, которого по модулю минимально, но больше чем значение (по модулю) текущего t (задается аргументом P). Код этой функции представлен ниже:

function FindMinimumT(var P:Extended; S:Integer): LongInt;
var
M:Extended;
K,i:LongInt;
begin
K:=0;
M:=1e24;
for i:=1 to SizeTParameter do
begin
if ((Abs(TParameter[i])<M) and (Abs(TParameter[i])>Abs(P)) and ((TParameter[i]*S)>=0)) then
begin
K:=i;
M:=Abs(TParameter[i]);
end;
end;
FindMinimumT:=K;
end;

Здесь нужно обратить внимание на аргумент S (может принимать следующие значения: –1, 0, 1). Он задает необходимый для поиска знак параметра, т. е. можно искать параметр t любого знака (S = 0), меньше нуля (S = -1) или больше нуля (S = 1). Обычно используется значение S = 0, но бывают случаи, когда необходимо отсекать значения параметра t определенного знака. Эти задачи будут рассмотрены позже.

1. Анализируются координаты возможной точки столкновения.

Для этого по формулам

находятся значения возможных координат x, y, z и проверяются на принадлежность поверхности, которой соответствует найденный выше параметр t. Например, значение координаты z для цилиндра должно быть больше 0 и меньше длины цилиндра. Если найденная точка не попадает в пределы соответствующей поверхности, то процесс повторяется с п. IV.2.2, причем в качестве текущего значения t выбирается только что найденный. Если же найденная точка попадает в пределы, то переходим к следующему пункту (IV.2.4).

1. Координаты точки столкновения частицы со стенкой системы x, y, z найдены.

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

1. Анализ взаимодействия частицы с поверхностью

Взаимодействие молекулы с поверхностью моделируется следующим образом. С помощью датчика случайных чисел генерируется случайное число μ. И если μ оказывается меньше или равно принятому коэффициенту прилипания ( ) на i-й поверхности, то молекула считается отраженной. Далее в зависимости от результатов анализа возможен переход на любую из следующих ветвей:

• Частица не осталась на поверхности.

Тогда (x₀, y₀, z₀) = (x, y, z). Переход в п. III (с. 11).

• Частица прилипла к поверхности.

На данном этапе происходит анализ различных параметров данной частицы, таких как координаты, номер поверхности на которой осталась частица, значения направляющих косинусов, время полета и т. д. Это необходимо для пополнения некоего «банка данных» анализируемой системы, для дальнейшей статистической обработки. Такими «данными» могут служить координаты точки, значения направляющих косинусов, с которыми частица пришла в данную точку, пройденное частицей расстояние и т. д. Далее проверяется условие окончания расчета – последняя ли это испытанная частица, если нет, то берется следующая частица (счетчик частиц увеличивается на 1) и производится переход в п. I (см. с. 7).

1. Статистический анализ и вывод окончательных результатов

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

,

где – количество поглощенных системой частиц; – общее количество испытанных частиц; – количество вылетевших из системы частиц. Коэффициент Клаузинга .

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

II.3.3. Нахождение полярных диаграмм скоростей частиц

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

Индикатрисы частиц строятся следующим образом. Выбирается поверхность, для которой необходимо построить полярные диаграммы скоростей частиц. Причем частицы должны пересечь эту поверхность или «перелететь» на нее. Обычно такой поверхностью служит входное или выходное сечение вакуумной системы – вход в крионасос, выход из трубопровода. При перелете на поверхность в процессе работы алгоритма уже найдены направляющие косинусы, которые определяют направление полета частицы к этой поверхности. Направляющие косинусы для перелета находились
в п. III (с. 11), а то, что частица перелетела на поверхность, для которой будут находиться индикатрисы определяется в п. IV.2.4. И именно данные из п. III (т. е. уже известные) используются для построения индикатрисы.

После выбора поверхности, для которой будут строится полярные диаграммы, определяется ось, параллельная нормали к этой поверхности. Наиболее часто индикатрисы определяются для поверхностей, перпендикулярных продольной оси OZ, поэтому именно эта ось в основном используется в качестве оси, вдоль которой анализируется структура потока. Для построения индикатрис, в процессе работы алгоритма необходимо накапливать значения направляющих косинусов частиц, которые пересекают заданную поверхность. Если анализируется поток вдоль оси OZ, то надо накапливать данные о направляющем косинусе n, который определяет угол между направлением полета и осью OZ. Для оси OX это будет направляющий косинус l, а для оси OY соответственно направляющий косинус m.

Далее весь угол от 0^о до 90^о делится на равные промежутки, например, на 90 равных частей по 1^о. Потом, в процессе работы алгоритма определяются направляющие косинусы и соответствующие им углы, и подсчитывается количество частиц, которые при пересечении заданной поверхности имели угол вылета от 0^о до 1^о, от 1^о до 2^о, от 2^о до 3^о и т. д. Таким образом, в результате получается некоторый массив, который содержит количество частиц, пересекших заданную поверхность под определенным углом. Фрагмент кода на Паскале, реализующий процесс такого накопления показан ниже:

CosineIndex:=Trunc(Abs(ArcCos(Ncos))*CosinesSize)+1;
Cosines[1,CosineIndex]:=Cosines[1,CosineIndex]+1;

Здесь Ncos – значение направляющего косинуса n; CosinesSize – количество частей, на которые дробится угол вылета от 0° до 90° (в данном случае CosinesSize:=90); Cosines – массив, который в зависимости от номера элемента будет содержать количество частиц, вылетающих под соответствующим углом: первый элемент содержит количество частиц, вылетевших под углом от 0° до 1°, второй – количество частиц, вылетевших под углом от 1° до 2°, третий – от 2° до 3°, и т. д.

Имея подобный массив, построение сводится к следующему. Нормаль к поверхности определяет угол в 0°. Касательная к поверхности определяет угол в 90°. В соответствии с найденным массивом от нормали к поверхности откладываются углы и на сторонах этих углов откладываются соответствующие значения количества частиц, вылетевших под этим углом (рис. II.3.3). На рисунке жирными линиями показано количество частиц для соответствующего угла. Далее конечные точки жирных линий соединяются и эта соединительная линия и образует индикатрису частиц, пролетающих через заданную поверхность ( см. рис. II.3.3 – штриховая линия).

Для построения индикатрисы частиц, пролетающих через заданную поверхность в единицах потока, не «число частиц под углом», а «число частиц под углом на единицу площади», необходимо элементы полученного массива умножить на площадь соответствующего телесного угла. Площадь телесного угла пропорциональна его синусу, поэтому можно умножить элементы массива на синус соответствующего угла:

For CosineIndex:=1 to CosinesSize do
Cosines[1,CosineIndex]:=Cosines[1,CosineIndex]*
Sin(CosineIndex*(Pi/2)/CosinesSize).

Здесь CosineIndex – номер элемента массива от 1 до CosinesSize.