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

Рис. II.3.5. Коаксиальный трубопровод

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

Прежде всего, пронумеруем поверхности. Входное сечение представляет собой плоскость, границами которой являются две окружности с радиусами r₁ и r₂, соответствующими радиусам внутреннего и внешнего цилиндров. Входное сечение – это поверхность № 1 системы. Внутренний цилиндр – литой, и, поэтому через него частицы не летят, он имеет радиус r₁. Внутренний цилиндр – поверхность № 2. Внешний цилиндр имеет радиус r_2, и он будет поверхность № 3. Последняя в данной системе поверхность № 4 – выходное сечение, которое представляет собой плоскость, аналогичную по форме поверхности № 1. Длина обоих цилиндров одинаковая и равна a. Для расчета распределения давления необходимо знать скорости частиц, значение входного потока Q [молекул/с], и промежуток времени Δt, через который будет происходить фиксация координат частиц. Примем значение коэффициента аккомодации равным 1, т. е. при ударе о поверхность частица сразу приобретает скорость, соответствующую температуре этой поверхности. Тогда, в нашем случае при ударе о поверхность № 2 частица приобретает скорость v₁, а при ударе о поверхность № 3 – v₂. При влете в структуру, частицы имеют скорость v₁. Значение Δt примем равным . В программе для краткости примем, что текущая скорость частицы выражается переменной RSpeed, а указанное значение
Δt – переменной dt.

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

Первым этапом алгоритма расчетного модуля является этап розыгрыша координаты частицы на входе в систему (см. параграф II.3.2). В нашем случае входное сечение представляет собой не круг, а кольцо, поэтому часть координат, попавших во внутреннюю окружность (Sqrt(x*x+y*y)<r₁) будет отброшена. Реализация этого этапа показана далее:

XSc:=0;
YSc:=0;
ZSc:=0;
RS:=r2;
MainLabel8:
GetStartL(CurrentX,CurrentY,CurrentZ,RS,XSc,YSc,ZSc);
if (Sqrt(CurrentX*CurrentX+CurrentY*CurrentY)<r1) then goto MainLabel8;

Здесь CurrentX, CurrentY, CurrentZ – текущие координаты частицы (на данном этапе – координаты частицы во входном сечении системы); XSc, Ysc, ZSc – координаты центра входного сечения; RS – радиус входного сечения (в нашем случае RS:=r₂).

Задаемся текущими значениями служебных переменных:

PossibleTValue:=1e-8;
PossibleTSign:=0;
CurrentSurface:=1;
RSpeed:=v1;
CurrentDistanceS:=0;

Здесь PossibleTValue – возможное значение параметра t; PossibleTSign – знак возможного значения параметра t; CurrentSurface – номер текущей поверхности; Rspeed – текущее значение скорости частицы.

Далее создается массив значений параметров t для случаев пересечения с разными поверхностями:

MainLa:
CreateTParameter;

Процедура CreateTParameter оперирует с глобальными переменными, для нашей задачи она показана ниже:

procedure CreateTParameter;
var
SurfaceType: Integer;
i,j: Integer;
Alfa,Gamma,D: Extended;
begin
if ((CurrentSurface=1) or (CurrentSurface=4) then
SurfaceType:=1;
if ((CurrentSurface=2) or (CurrentSurface=3) then
SurfaceType:=2;
GetAlfaGamma(Alfa,Gamma,CurrentSurface);
UpdateDirect(SurfaceType,Alfa,Gamma);
D:=0;
TParameter[1]:=GetLinearT(D,CurrentZ,NCos);
TParameter[2]:=0;
TParameter[3]:=GetCylinderTM(CurrentX,CurrentY,LCos,
MCos,r1);
TParameter[4]:=GetCylinderTP(CurrentX,CurrentY,LCos,
MCos,r1);
TParameter[5]:=GetCylinderTM(CurrentX,CurrentY,LCos,
MCos,r2);
TParameter[6]:=GetCylinderTP(CurrentX,CurrentY,LCos,
MCos,r3);
D:=a; \\ a – длина цилиндров
TParameter[7]:=GetLinearT(D,CurrentZ,NCos);
TParameter[8]:=0;
end;

Суть данной процедуры в том, что в зависимости от типа текущей поверхности (SurfaceType, 1 – плоскость, 2 – цилиндр) определяются направляющие косинусы (глобальные переменные LCos, MCos, Ncos): вначале определяются углы Alfa и Gamma, характеризующие положение текущей координаты (функция GetAlfaGamma), а потом определяются направляющие косинуса в зависимости от типа текущей поверхности. Так как у поверхностей второго порядка (конус, цилиндр, сфера, и т. д.) есть две точки пересечения с прямой (два значения параметра t), вторые значения параметра t для плоскостей приравниваются к нулю.

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

MainLabel2:
PossibleTIndex:=FindMinimumT(PossibleTValue,PossibleTSign);
\\ поиск минимального по модулю значения t заданного
\\знака
if (PossibleTIndex=0) then
begin
ErrorsNumber:=ErrorsNumber+1;
goto MainLabel1;
end;
PossibleSurface:=GetSurfaceNumber(PossibleTIndex); \\ определения номера возможной поверхности перелета по \\номеру
PossibleTValue:=TParameter[PossibleTIndex]; \\ элемента в массиве значений t
PossibleX:=CurrentX+LCos*PossibleTValue;
PossibleY:=CurrentY+MCos*PossibleTValue;
\\ определение возможных координат следующей точки \\соударения
PossibleZ:=CurrentZ+NCos*PossibleTValue;
if (TestCoordinate(PossibleX,PossibleY,PossibleZ,PossibleSurface)=2)
\\ проверка принадлежности данных возможных координат
then \\ соответствующей поверхности
begin
if (GetTClone(CurrentSurface)=PossibleSurface) then \\ проверка: не является ли данная поверхность «парной»
begin
PossibleTSign:=-GetSign(PossibleTValue);
\\ и если является, то задается знак поиска в массиве t
goto MainLabel2;
end;

Необходимо отметить один важный момент, показанный в данном фрагменте. После проверки на принадлежность возможной точки взаимодействия, соответствующей выбранному значению параметра t поверхности, производится проверка на «парность» поверхности, с которой частица вылетает, по отношении к поверхности, на которую частица может попасть. Это необходимо делать в случаях, когда необходимо, чтобы частица не могла пролететь сквозь какую-либо поверхность. Например, поверхность № 2 – литой, непроницаемый цилиндр, и поэтому через него частицы пролетать не могут. Соответственно, он является «парной» поверхностью для самого себя. Парные поверхности также могут быть противоположными стенками непроницаемой (например, металлической) детали. Ниже показаны функции для определения «парной» поверхности, проверки принадлежности возможной координаты соответствующей поверхности и определения номера поверхности по номеру значения параметра t в массиве для данной задачи:

function GetTClone(var N:LongInt): LongInt;
begin
if (N=2) then GetTClone:=2 else GetTClone:=0;
end;
function TestCoordinate(var X,Y,Z:Extended; SN:LongInt): Integer;
begin
TestCoordinate:=1;
if ((SN=1) or (SN=4) then
begin
if ((Sqrt(X*X+Y*Y)<=r2) and (Sqrt(X*X+Y*Y)>=r1)) then TestCoordinate:=2;
end;
if ((SN=2) or (SN=3) then
begin
if ((Z<=a) and (Z>=0)) then TestCoordinate:=2;
end;
end;
function GetSurfaceNumber(var P:LongInt): LongInt;
begin
if ((P mod 2)=1) then GetSurfaceNumber:=(P+1) div 2 else GetSurfaceNumber:=P div 2;
end;