Генерация вектора скорости молекулы и координат точки влета
Основные принципы моделирования вектора скорости молекул набегающего потока изложены в пункте 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), которая является косинусом угла между направлением вектора те-пловой скорости и положительным направлением оси Х. Выражение:
vxT = VT×c (22)
определяет значение проекции вектора скорости на ось Х.
б) Проекция вектора температурной скорости на плоскость YZ определяется выражением:
vyzT = VT sin(arccos(c)).
Учитывая тригонометрическую формулу sin2c + cos2c= 1, получаем:
(23)
В плоскости YZ направление вектора скорости равновероятно и определяется случайным углом e , который откладывается от положительного направ-ления оси Y. Угол e выбирается при помощи генератора случайных чисел в диапазоне (0,2p). Выражения для проекций вектора скорости на оси Y и Z имеют вид:
vyT = vyzT × cos(e),
vzT = vyzT × sin(e) (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)
,
где r - радиус соответствующей окружности.
Решая систему (28), получаем:
(y1+vyT)2+(z1+vzT)2=r2 (29)
Решив это уравнение относительно Т, получаем:
(30)
или :
|
,
где tmax и tmin - времена пересечения молекулой боковой поверхности цилиндра радиуса r.
Если значения tmax и tmin - комплексные числа, то траектория движения молекулы не пересекает рассматриваемую поверхность; положительные значения tmax и tmin означают, что молекула движется извне по направлению к рассматриваемой поверхности; отрицательные значения tmax и tmin означают, что молекула находится вне рассматриваемого цилиндра и удаляется от него; если tmax >0, а
tmin <0, то начальная точка движения находится внутри рассматриваемого цилиндра.
Времена пересечения траектории движения молекулы с поверхностями торцов цилиндра определяются по проекции траектории движения молекулы на ось Х.
Координатами торцов по оси X являются: х = 0 (для переднего торца) и
х = l (для заднего торца), тогда времена пересечения траектории движения молекулы с торцами t0и t1 определяются выражениями:
. (32)
Для нахождения конечной точки прямолинейной траектории движения молекулы необходимо из времен tmax(r0), tmin(rA), t0 и t1 выбрать наименьшее положительное время -tкон. И тогда:
x2 = vx × tкон + x1
y2 = vy × tкон + y1 (33)
z2 = 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, определяющие поверхность, на которую попадает молекула.
1.1.10. Расчет времен пребывания молекулы в частных объемах Wpk
Расчет времен пребывания молекулы в частном объеме необходим для расчета относительного распределения концентраций молекул внутри цилиндра. Как показано в разделе 2, расчетный объем поделен на 2 × ks × ps частных объемов Wpk, времена нахождения молекулы в которых, есть разность времен влета и вылета.
Координаты влета и вылета молекул из частного объема и соответствующие им времена рассчитываются путем решения системы из проекций на плоскость YZ и ось Х уравнения движения молекулы (27) и уравнений поверхностей, образующих частный объем. Времена влета и вылета из частного объема должны быть в пределах начального и конечного времени нахождения молекулы в расчетном объеме.
Рассмотрим процедуру нахождения координат и времен пересечения траектории с образующими частный объем Wpk плоскостями. В проекции на плоскость YZ частный объем представляется в виде кольцевого сектора, внешняя и внутренняя дуги которого - части окружностей, соответствующие границам расчетного объема. Координаты пересечения траектории движения молекулы с радиальными границами кольцевого сектора рассчитываются из системы уравнений движения и образующей сектора, описываемой выражением:
y×sinl - z×cosl = 0, (34)
где l - угол между осью Y и радиальной границей сектора.
В общем виде система уравнений может быть записана следующим обра-зом:
y × vz – z × vy = vz × y1 – vy × z1
(35)
y × sin l - z × cos l = 0
Решением системы уравнений являются координаты:
(36)
|
,
|
,
т.к. рассматриваем только положительную часть радиальной границы сектора.
Время пересечения траекторией боковой поверхности частного объема рассчитывается из выражения:
. (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.
Описанный механизм расчета и выбора повторяется для каждого прямолинейного участка траектории молекулы, и, суммируя времена для большого количества испытаний, получаем необходимые статистические данные для расчета распределения концентраций молекул в расчетном объеме.