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

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

.

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

После того, как вычислена функция распределения f_x для найденного случайного значения x, вычисляется функция . Затем генерирует очередное случайное число . Выбранное значение x либо принимается, либо отвергается в зависимости от того, больше или меньше , чем . Если < , то значение x принимается, как удовлетворяющее принятому распределению. Если >, то значение x не принимается, и операция повторяется: генерируется новое значение (значения) x, вычисляется f_x, вычисляется , и т. д.

Для функции распределения Максвелла , где , T – температура молекулы, v₁ – одна из трех компонент скорости . Равномерное распределение значений v задается формулой . При этом для a и b выбираются некоторые конечные значения вместо реальных пределов от до + . Если положить a и b равными и соответственно, то часть значений, лежащих за этими пределами, равна 1-erf (10), т. е. очень мала. Таким образом, и . Генерируется очередное значение ξ, если , то v₁ принимается. Если , то значение v₁ отвергается, и весь процесс повторяется до тех пор, пока очередное значение не будет принято.

Затем аналогичные процедуры повторяются для второй и третьей компоненты скорости v₂ и v₃. Процедура генерации трех компонент скорости по распределению Максвелла показана ниже:

Procedure GetMaxwellDistr(Var V1, V2, V3, Betta: Extended);
Var
FsV: Extended;
Label
V1R;
Label
V2R;
Label
V3R;
begin
V1R:
V1:=(-10+20*Random)/Betta;
FsV:=Exp(-V1*V1*Betta*Betta);
if (FsV<Random) then goto V1R;
V2R:
V2:=(-10+20*Random)/Betta;
FsV:=Exp(-V2*V2*Betta*Betta);
if (FsV<Random) then goto V2R;
V3R:
V3:=(-10+20*Random)/Betta;
FsV:=Exp(-V3*V3*Betta*Betta);
if (FsV<Random) then goto V3R;
end;

II.3.8. Определение распределения концентрации и давления

Для определения распределения давления внутри анализируемой системы существует два основных подхода [2,3].

Первый состоит в нахождении числа соударений частиц с фрагментом стенки вакуумной системы. Затем, на основании этой информации находится распределение давления. Главным недостатком такого подхода является то, что полученный результат характеризует распределение давления только вблизи стенок вакуумной системы. Однако такой метод вполне применим, если необходимо предсказать значение давления, которое будет показано датчиком, прикрепленным к стенке, например датчиком, прикрепленным к стенке трубопровода. Достоинством данного подхода является крайняя простота его реализации. Анализируемая система делится на несколько зон. Например, трубопровод делится на 10 одинаковых частей. Затем в процессе расчета проводимости подсчитывается количество ударов (взаимодействий) частиц с каждой из 10 частей. Удобно сделать этот расчет в зависимости от продольной координаты z: если a₁ < z < b₁ – первая часть, a₂ < z < b₂ – вторая, и т. д. Зная общее количество пробных частиц (например, 10000), получаем следующий массив , где N_i – количество взаимодействий (соударений) в i-й зоне; N_общ – общее количество пробных частиц; F_i – площадь i-й части. Затем, умножив каждый элемент этого массива на значение общего потока молекул Q [молекул/с], получим значения числа соударений молекул в секунду о единицу площади для каждой части . Зная молярную массу молекул и температуру газа, можно получить значения давления вблизи стенок для каждой части , где k – постоянная Больцмана; μ – молярная масса; T – температура газа; R – универсальная газовая постоянная. Полученный массив значений давления p_i для нескольких зон характеризует распределение давления вблизи стенок в данных зонах. Подсчет числа соударений удобно вести в том месте алгоритма, где происходит анализ взаимодействия частицы с поверхностью.

Реализация второго подхода более сложна, однако он дает возможность определять трехмерное поле концентрации. Суть данного подхода состоит в следующем. Задается некоторое значение времени Δt, с условием, чтобы , где l – характерный линейный размер анализируемой вакуумной системы; v – скорость молекул. Причем лучше, если будет выполняться . После этого вся система так же, как и в предыдущем методе, разбивается на части, с той лишь разницей, что теперь эти части объемные, т. е. характеризуются тремя координатами. В начале расчета задается скорость частиц. После запуска частицы начинается подсчет пройденного ею расстояния. Для расчета пройденного расстояния удобно пользоваться параметром t, поскольку модуль его значения – это расстояние от точки старта частицы до точки финиша. Так, если частица вылетела с поверхности из точки M₀(x₀, y₀, z₀) и перелетела на другую поверхность в точку M(x, y, z), то пройденное ею расстояние от M₀ до M численно равно . Затем в реализацию алгоритма расчета необходимо включить часть, отвечающую за определение координат частицы, которые она будет иметь через время Δt после начала отсчета времени. Для этого необходимо считать время, в течение которого частица находится в полете. Это удобно сделать, используя значение модуля параметра t – . Так, время, которое частица находится в полете, рассчитывается по формуле: , где начальное значение . Эта процедура выполняется до тех пор, пока . Перед каждым суммированием выполняется проверка. Если , то выполняется суммирование и значение увеличивается, а если , то значение параметра подбирается таким образом, чтобы , и для этого значения параметра t определяются координаты частицы. При подборе необходимого значения параметра t нужно учитывать то, что знак нового параметра должен быть таким же, как и у старого. Подбор нового значения параметра t_n выполняется по соотношению: . Здесь Δt – время, через которое фиксируются координаты частицы; – общее время полета частицы; v – скорость частицы;
– соотношение, характеризующее знак старого значения параметра t (данное соотношение может принимать значение +1 или –1). После определения значения параметра t_n, рассчитываются координаты частицы, которые она будет иметь через время Δt после начала полета:

Используя эти координаты, можно определить в какой из частей системы находится частица и, аналогично предыдущей методике, увеличить счетчик частиц в этой части (зоне): N_i=N_i+1. Далее с использованием сохраненного старого значения параметра t находятся истинные координаты, соответствующие поверхности, на которую попала частица:

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

По окончании процесса расчета полученные значения N_i для каждой части используются для определения концентрации, [молекул/м³]: , где Q – значение потока [молекул/с]; Δt – значение отрезка времени, через которое фиксируются координаты частиц; N_общ – общее число пробных частиц (обычно 10000 или 100000); V_i – объем i-й части, [м³]; N_i – число частиц в i-й части. Поскольку каждая частица может несколько раз зафиксироваться через промежутки времени Δt, то фактически получается, что одна частица моделирует себя и несколько других частиц, вылетевших по данной траектории через время Δt. Данный подход позволяет моделировать квазипостоянный напуск частиц в систему, уменьшая значение Δt, т. е. если в традиционном методе пробной частицы происходит единственный мгновенный напуск, например, 10000 частиц, и затем производится анализ их полета, то при использовании данного подхода происходит установка 10000 траекторий, по которым через равные промежутки времени Δt запускаются 10000 частиц в систему. Поэтому, с уменьшением значения Δt, уменьшается промежуток между мгновенными запусками и процесс напуска становится квазипостоянным.

II.3.9. Пример расчета параметров коаксиального трубопровода

В качестве иллюстрации применения метода Монте-Карло для анализа сложных вакуумных систем рассмотрим следующую задачу. Имеется вакуумная система, состоящая из двух коаксиальных цилиндров – один из которых полый – внешний, имеет температуру, соответствующую тепловой скорости молекул v₂, а второй – внутренний, литой (рис. II.3.5), имеет температуру, соответствующую скорости молекул v₁. Стрелками показаны направления влета и вылета частиц из данной системы.

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