1.1.14 Описание программы моделирования MODMD05

На базе существующей программы моделирования взаимодействия молекул набегающего потока и СВА с конструктивными элементами MODMD05, была создана программа моделирования MODMD24, в которой введен дополнительный коэффициент, характеризующий исследуемый поток молекул: T_mm - температура направленного потока молекул с Максвелловским законом распределения модуля скоростей. При этом появляется возможность дополнительно моделировать набегающий поток, создаваемый собственными газовыделениями КА.

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

1.2.1 Моделирование потока собственных газовыделений. КА

В программе MODMD24 выражение для расчета модуля направленной скорости молекул, которое в MODMD05 описывается выражением (20), реализуется в виде:

, (52)

где W₀ - направленная скорость молекул, нормированная на модуль наиболее вероятной скорости молекулы, имеющей температуру стенок исследуемого объема: V_ст0² = R_T T; V₀ - модуль направленной скорости молекул, влетающих в исследуемый объем; T - температура стенок исследуемого объема; R_T - значение газовой постоянной для молекулярного азота; V - величина, распределенная по закону Максвелла; T_mm - температура направленного потока молекул с Максвелловским законом распределения модуля скорости.

При значении T_mm=0 выражение (52) сводится к исходному выражению (20).

В результате, новая модификация программы MODMD24 дает возможность моделировать все виды потоков нейтральных молекул, воздействующих на исследуемый МИП:

- при T_mm 0, T_m=0 и V₀=0 моделируется воздействие на МИП потоков собственных газовыделений систем КА;

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

- при Tmm=0, Tm 0, V₀ > 0 моделируется воздействие разреженной земной атмосферы для высот полета КА на МИП.

1.2.2 Моделирование влета (вылета) молекул в датчик

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

Так как отверстия в переднем и заднем торцах одинаковы по форме и площади, интегральный коэффициент прозрачности торцевых стенок определяется выражением:

(53)

где r01, r02 - радиусы, определяющие размеры отверстий в переднем и заднем торцах МИП; ugol - угол, определяющий размеры отверстий МИП; kpr - коэффициент прозрачности отверстий МИП.

Рис.8 Модель датчика с учетом параметра ugol

Выбор точек влета-вылета молекул в исследуемый объем определяется геометрическим положением открытой области. Розыгрыш точек влета происходит в несколько этапов:

а) В зависимости от знака компоненты скорости vx, определяется торец, через который молекула влетает в исследуемый объем, а, следовательно, координата x₂.

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

в) При помощи генератора случайных чисел выбираются значения координат y₂ и z₂ точки влета в диапазоне (-1,1) до тех пор, пока выбранная точка не окажется в области одного из отверстий.

Для этого необходимо, чтобы совместно выполнялись условия:

y₂² + z₂² < r₀₂², (54)

y₂² + z₂² > r₀₁², (55)

(при sk = - 1) (56)

или

(при sk = 1) (57)

xkpr < kpr

где sk - параметр, принимающий значение "-1" при розыгрыше влета в передний торец, и значение "1" при розыгрыше влета в задний торец.

Первое слагаемое в выражениях (56) задает косинус половины угла, определяющего размер входного отверстия (в радианах), второе - косинус угла, под которым предполагаемая точка влета молекулы расположена по отношению к координатной оси Y.

При одновременном соблюдении условий (54), (55), (56), (57) молекула считается влетевшей в исследуемый объем.

Условия (54), (55), (56), (57) должны выполняться и при розыгрыше вылета молекулы из объема.

В связи с описанными дополнениями изменяется структура соответствующего файла исходных данных. Исходные данные для программы MODMD24 содержатся в файле inp24. txt и включают три дополнительных параметра - Tm, Tmm и ugol.

Исходные данные записываются в файле inp24. txt в следующем порядке:

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

r₁ - внутренний радиус расчетного объема, отнесенный к радиусу исследуемого объема;

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

ugol [град.] - угол, определяющий размер отверстия в переднем (заднем) торце;

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

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

Tm - абсолютная температура разреженного газа;

Tmm - температура направленного потока молекул, скорости которых распределены по закону Максвелла;

ugol_x (_v) - угол между направлением вектора направленной скорости влетающих молекул и продольной осью МИП (Х);

ugol_у (v) - угол между проекцией вектора направленной скорости влетающих молекул на плоскость YZ, перпендикулярную продольной оси X, и осью Y [1];

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

Таким образом, строка исходных данных имеет следующий примерный вид:

10000 0.3 1.0 180 0.95 8000 300 0 0 90 1

1.3 Описание программы MODMD79

1.3.1 Разработка модернизированной математической модели

Основной идеей, положенной в основу модернизации программы MODMD24, является гипотеза о пропорциональности относительной концентрации в объеме сумме отношений времен пролета молекул в исследуемом объеме и свободном пространстве.

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

Разобьем все влетающие частицы на классы: в класс j входят все частицы с определенной точкой влета, определенной скоростью влета, одинаково рассеиваемые в объеме 1. Пусть при зеркальном отражении в рассматриваемом объеме таких частиц будет N_j, тогда при диффузном отражении число частиц класса j в объеме 1 будет равно:

,

где t_реалj - среднее время нахождения частиц класса j в объеме 1;

t_идеалj - время нахождения частицы класса j в объеме 2.

Очевидно, что

Тогда, сложив частицы всех классов, получим, что в объеме 1 количество частиц равно:

,

где N - количество рассматриваемых частиц (число частиц в объеме 2).

При зеркальном отражении плотность внутри объема равна внешней плотности n = N/V, где V - величина рассматриваемого объема.

Рассмотрим относительную плотность в объеме 1.

(58)

Таким образом (58) - расчетная формула для всего объема.

Если объем 1 разбит на несколько частей, то отношение количества частиц класса j в частном объеме к количеству частиц класса j в объеме 1 равно t_частj/t_реалj, где t_{частj -} время нахождения в частном объеме частиц класса j. Если N_j - количество частиц класса j в объеме 1, то количество частиц класса j в частном объеме

Очевидно, что

Просуммируем по всем классам j, тогда

Плотность в частном объеме естественно задавать как

,

тогда относительная плотность в частном объеме

. (59)

Таким образом, (59) - расчетная формула для частного объема.

Покажем, что (59) соответствует формуле (58). Пусть все частные объемы одинаковы по величине, тогда

По определению:

При равных частных объемах формула (60) очевидна.

Достоинством данного подхода является то, что для его реализации достаточно внести в текст программы MODMD24 лишь небольшие изменения.

1.3.2 Разработка программы моделирования

В соответствии с изложенными теоретическими идеями была создана программа MODMD79. Модификации в текст MODMD24 вносились таким образом, чтобы основное расчетное ядро оставалось без ощутимых изменений:

а) добавлены новые переменные:

- tid - время пролета молекулы в зеркально отражающей трубе;

- q1 - новое значение вспомогательного коэффициента q;

- Tsss - старое значение переменной Tss, теперь она вычисляется по другой формуле;

- NNN - значение относительной концентрации во всем объеме, рассчитанное по новому методу;

- nk1 [] - массив, хранящий значения относительной концентрации в частных объемах, рассчитанные по новому методу, nk [] теперь вспомогательный массив, содержит периоды пролета текущей молекулы в частных объемах;

б) изменения в расчетных процедурах: в конце процедуры counttrack в Tss теперь накапливается величина taukon/tid, в Tsss сохранено старое значение Tss;

в) изменения в теле основной программы:

- значение q1 вычисляется по формуле: (r₁² - r₂²) (1 - r_A²) / (2 ks ps a);

- для каждой конкретной молекулы после моделирования влета определяется величина tid = |L/vx|, обнуляется массив nk [];

- после вылета текущей молекулы значения nk [] делятся на tid, полученный nk [] суммируется с nk1 [];

- после отработки всех n траекторий вычисляется NNN по сумме всех элементов nk1 [] и вспомогательным коэффициентам (a, q1);

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

2. Доработка математической модели торцевых стенок

Рис.12 Распределение плотности при отклонении от оси Х на 5 градусов по программе MODMD82

Рис.13 Распределение плотности при отклонении от оси Х на 10 градусов по программе MODMD82

Рис.14 Распределение плотности при отклонении от оси Х на 15 градусов по программе MODMD82

Рис.15 Распределение плотности при нахождении потока перед датчиком по программе MODMD82krug

Рис.16 Распределение плотности при нахождении потока перед датчиком по программе MODMD82krug

Рис.17 Распределение плотности при нахождении потока перед датчиком по программе MODMD82krug

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

Интерпретация выходных данных

В результате анализа графиков можно прийти к выводу, что самым оптимальным случаем является тот, при котором углы отклонения от оси Х стремятся к нулю, тогда угол раскрытия торцов можно изменять в диапазоне от 80 до 190 градусов. В этом случаи минимальное значение концентрации N1 отличается на 10 процентов от максимального значения, а N2 остается практически неизменной, достигая своего максимум.

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

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

Отклонение на 10 или 15 градусов ведет к существенным изменениям параметров. При дальнейших отклонениях следует обратиться к задаче обнаружения потока при углах раскрытия торцов 90 или 180 градусов.

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

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

Действительно, напряжение на выходе измерителей тока:

(64)

где и - коэффициенты преобразования верхней и нижний частей катода и подключенных к ним измерителей тока;

N1 и N2 - концентрация газа в нижней (Ndw) и верхней (Nup) частях ионизационных камеры, образованной частями катода.

При этом:

(65)

где Nc - концентрация равновесного газового потока;

Nv - концентрация направленного по оси динамического газового потока;

k1,l2,k2,l2 - коэффициенты пропорциональности концентрации в первой и второй частей ионизационный камеры для каждого вида газового потока соответственно.

Подставляем значения из (65) в (64), получаем:

(66)

Решение системы уравнений (3) имеет вид:

(67)

При симметричной конструкции ионизационной камеры и одинаковых характеристиках измерителей тока можно принять ==, k1=k2=k, тогда:

(68)

где A,B,C - коэффициенты, определяемые конструтивными характеристиками устройства;

Uraz= U2 - U1 - разностный сигнал двуз измерителей тока;

Usum= U1 + U2 - суммарный сигнал двуз измерителей тока;

Коэффициенты A,B,C могут быть получены при градуировке устройства и расчетным путем. Для этого при отсутствии динамического газового потока Nv = 0, Uraz = 0 при известной концентрации равновесного газового потока Nc определяется коэффициент A. Затем при известной концентрации равновесного газового потока Nc и динамического газового потока Nv с заданными динамическими характеристиками определяются коэффициенты B и C.

Заключение

{MODMD82. PAS - реализация метода АБ}

{Новые переменные:

q1 - новое значение вспомог. к-та q;

tid - (по идее) время пролета в зеркальноотражающем объеме;

Tsss - старое значение Tsss, теперь Tss выч. по др. формуле;

NNN - относительная концентрация в объеме по методу АБ;

nk1 [] - распределение отн. конц. по методу АБ. }

{Изменения в процедурах: только одна строчка в конце counttrack}

{18.04.2008: Мелкие несущественные изменения (интерфейс) }

program dat (input,output);

label ll; const ks=12; ps=3; ProgID=82; r11=0.600; r12=0.240; s=3;

var a,q,vx,vy,vz,vyz,v,vT,vTm,vTmm,v0,w0,vv,y,y0,teta,xkpr,ksi,vn,ww0,q1,su1,sum1,su2,sum2: real;

x2,y2,z2,vx00,vyz00,vz00,vy00,Tm,Tmm,r0,r1,r2,ugol,ac,ugol_x,ugol_y,Tsr,Tss,Vsr,VVsr: real;

vxTmm,vxTmmv,vyzTmm,vyTmm,vyTmmv,vzTmm,vzTmmv,tid,Tsss,NNN, R,T,L,rA,r01,r02,kpr: real;

j,p,m1,k,kR,kL,f3,f4,sk: integer; n, i,n1,n2: longint; isp_No,m3,tekp: integer; defl: boolean;

inp1,out2,out3,out4: text; Fdir,conc,Nup,Ndn: real;

nk: array [1.2*ks,1. ps] of real;

nk1: array [1.2*ks,1. ps] of real; nc: array [1. ps,1.2] of real; d: array [1. ps] of real;

procedure Maxwell;

begin

repeat

v: = (random (5) +random);

y: =v*v*exp (1-v*v);

y0: =random;

until y0<=y;

if (v<=0.05) then v: =0.05;

end;

procedure scattering (Mx,My: real);

var vx0,vyz0,vy0,vz0,gamma,beta,alfa: real;

begin

repeat

beta: =pi/2*random;

y: =2*sin (beta) *cos (beta);

y0: =random;

until y0<=y;

Maxwell;

vv: = (1-ac) *vv*vv+ac*v*v+ac* (1-ac) *vv*vv* ( (v*v) /1.5-1);

if (vv<=0.0025) then vv: =0.05

else vv: =sqrt (vv);

alfa: =2*pi*random;

vx0: =vv*cos (beta); vyz0: =vv*sin (beta);

vy0: =vyz0*cos (alfa); vz0: =vyz0*sin (alfa);

if Mx<2 then begin vx: =Mx*vx0; vy: =vy0; vz: =Mx*vz0; end

else if not (y2*z2=0) then begin

if (y2*My) >0 then gamma: =arctan (z2/y2) +Pi-arctan (vy0/vx0)

else gamma: =arctan (z2/y2) - arctan (vy0/vx0);

vx: =-vz0; vy: =sqrt (vx0*vx0+vy0*vy0) *cos (gamma);

vz: =sqrt (vx0*vx0+vy0*vy0) *sin (gamma); end

else if y2=0 then begin

vx: =-vz0; vy: =-vy0*My*z2/abs (z2); vz: =-vx0*My*z2/abs (z2); end

else begin vx: =-vz0; vy: =-vx0*My*y2/abs (y2); vz: =vy0*My*y2/abs (y2); end;

end;

procedure counttrack;

label 2,3;

var x1,y1,z1,tmin1,tmax1,tmin2,tmax2,tmin,tmax,taumin,taumax,taukon,tAnod: real;

f0,f1,f2,fr1,fr2,fa2: real;

procedure countmolec;

var w: 0. ps; b1,b2,b3,b4,a1,a2: boolean;

d1,d2,d3,d4,tau1,tau2,tau3,tau4,tkon1,tkon2,t0,t00: real; label ii;

begin

Tsr: =Tsr+abs (tmax-tmin);

for w: =0 to ps-1 do

begin

tkon1: = ( (L*w/ps) - x1) /vx; tkon2: = ( (L* (w+1) /ps) - x1) /vx;

p: =w+1;

a1: = ( (tkon1>=tmin) and (tkon1<=tmax));

a2: = ( (tkon2>=tmin) and (tkon2<=tmax));

if a1 then if a2 then if vx>0 then begin t0: =tkon1; t00: =tkon2; end

else begin t0: =tkon2; t00: =tkon1; end

else if vx>0 then begin t0: =tkon1; t00: =tmax; end

else begin t0: =tmin; t00: =tkon1; end

else if a2 then if vx>0 then begin t0: =tmin; t00: =tkon2; end

else begin t0: =tkon2; t00: =tmax; end

else if (tkon1<tmin) and (tkon2>tmax) or (tkon1>tmax) and (tkon2<tmin)

then begin t0: =tmin; t00: =tmax; end

else goto ii;

d1: =y1-z1*vy/vz; d2: =d1;

tau1: = (d1-y1) /vy; tau2: = (d2-y1) /vy;

for kR: =1 to ks do

begin

d3: =vz*y1-vy*z1;

d4: =d3/ (vz*cos (kR*Pi/ks) +vy*sin (kR*Pi/ks));

d3: =d3/ (vz*cos (kR*Pi/ks) - vy*sin (kR*Pi/ks));

tau3: = (d3*cos (kR*Pi/ks) - y1) /vy; tau4: = (d4*cos (kR*Pi/ks) - y1) /vy;

b1: = ( (d1>0) and (tau1<=t00) and (tau1>=t0));

b2: = ( (d2>0) and (tau2<=t00) and (tau2>=t0));

b3: = ( (d3>0) and (tau3<=t00) and (tau3>=t0));

b4: = ( (d4>0) and (tau4<=t00) and (tau4>=t0));

if b1 then if b3 then nk [kR,p]: =nk [kR,p] +abs (tau1-tau3)

else if ( (tau3-tau1) *d3) >0 then nk [kR,p]: =nk [kR,p] +abs (tau1-t00)

else nk [kR,p]: =nk [kR,p] +abs (tau1-t0)

else if b3 then if ( (tau1-tau3) *d1) >0 then nk [kR,p]: =nk [kR,p] +abs (tau3-t00)

else nk [kR,p]: =nk [kR,p] +abs (tau3-t0)

else if ( ( (d1>0) and (d3>0) and ( ( (tau1>t00) and (tau3<t0)) or ( (tau1<t0) and (tau3>t00))))

or ( (d1*d3<0) and ( ( (d3* (tau1-tau3) >0) and (tau1>t00) and (tau3>t00))

or ( (d1* (tau3-tau1) >0) and (tau1<t0) and (tau3<t0)))))

then nk [kR,p]: =nk [kR,p] +abs (t00-t0);

kL: = (2*ks+1) - kR;

if b2 then if b4 then nk [kL,p]: =nk [kL,p] +abs (tau2-tau4)

else if ( (tau4-tau2) *d4) >0 then nk [kL,p]: =nk [kL,p] +abs (tau2-t00)

else nk [kL,p]: =nk [kL,p] +abs (tau2-t0)

else if b4 then if ( (tau2-tau4) *d2) >0 then nk [kL,p]: =nk [kL,p] +abs (tau4-t00)

else nk [kL,p]: =nk [kL,p] +abs (tau4-t0)

else if ( ( (d2>0) and (d4>0) and ( ( (tau2>t00) and (tau4<t0)) or ( (tau2<t0) and (tau4>t00))))

or ( (d2*d4<0) and ( ( (d4* (tau2-tau4) >0) and (tau2>t00) and (tau4>t00))

or ( (d2* (tau4-tau2) >0) and (tau2<t0) and (tau4<t0)))))

then nk [kL,p]: =nk [kL,p] +abs (t00-t0);

tau1: =tau3; tau2: =tau4; d1: =d3; d2: =d4;

end;

ii: end;

end;

begin

tmax: =0; tmin: =0;

x1: =x2; y1: =y2; z1: =z2;

f0: =vy*vy+vz*vz; f1: = (vy*y1+vz*z1) /f0;

f2: =f0-sqr (vz*y1-vy*z1);

fa2: =rA*rA*f0-sqr (vz*y1-vy*z1);

fr1: =r1*r1*f0-sqr (vz*y1-vy*z1);

fr2: =r2*r2*f0-sqr (vz*y1-vy*z1);

if f2>0 then

begin

f2: =sqrt (f2) /f0; taumin: =-f1-f2; taumax: =-f1+f2;

if not (vx=0) then if vx>0 then begin taukon: = (L-x1) /vx; f3: =-1; end

else begin taukon: = (0-x1) /vx; f3: =1; end

else taukon: = (2/vyz) +0.1;

tAnod: =-1;

if fr2>0 then

begin

fr2: =sqrt (fr2) /f0; tmin2: =-f1-fr2; tmax2: =-f1+fr2;

if fr1>0 then

begin

fr1: =sqrt (fr1) /f0; tmin1: =-f1-fr1; tmax1: =-f1+fr1;

if fa2>0 then begin fa2: =sqrt (fa2) /f0; tAnod: =-f1-fa2; end;

if (tmin2>0) and (tmin2<taukon) then Tmin: =tmin2

else if (tmin2<=0) and (tmin1>0) then Tmin: =0

else goto 2;

if (tmin1<taukon) then Tmax: =tmin1

else Tmax: =taukon;

countmolec;

2: if (tmax1>0) and (tmax1<taukon) then Tmin: =tmax1

else if (tmax1<=0) and (tmax2>0) then Tmin: =0

else goto 3;

if (tmax2<taukon) then Tmax: =tmax2

else Tmax: =taukon;

if (tAnod<0) then countmolec;

goto 3;

end;

if (tmin2>0) and (tmin2<taukon) then Tmin: =tmin2

else if (tmin2<=0) and (tmax2>0) then Tmin: =0

else goto 3;

if (tmax2>taukon) then Tmax: =taukon

else Tmax: =tmax2;

countmolec;

end;

3: end;

if taumax<taukon then begin taukon: =taumax; f3: =2; f4: =1; end;

if (tAnod>0) and (tAnod<taukon) then begin taukon: =tAnod; f3: =2; f4: =-1; end;

x2: =x1+taukon*vx; y2: =y1+taukon*vy; z2: =z1+taukon*vz;

Tss: =Tss+taukon/tid; Tsss: =Tsss+taukon; {в Tss накапливаем taukon/tid; Tsss - старое значение Tss}

if x2>L then x2: =L; if x2<0 then x2: =0;

end;

procedure byben;

var y0,y,sum: real;

begin

m3: =0;

if (r0>r11-r12) and (r0<r11+r12) then

begin

sum: = (r11+r12+r0) /2;

y0: =1/3+4/pi*2*arctan (sqrt ( (sum-r11) * (sum-r12) * (sum-r0) /sum) / (sum-r12));

y: =random;

if (y<=y0) then m3: =1;

end;

end;

{procedure vvod (var Nsu: mas);

begin

assign (out4,'c: \80Nsum. txt');

reset (out4);

for i1: =1 to s do

for j1: =1 to m do

read (out4,Nsu [i1,j1]);

{writeln (a [1,2]: 7: 4);

writeln (a [2,2]: 7: 4);

writeln (a [3,2]: 7: 4);

close (out4);

end; }

function summa: real;

var su1,sum1,su2,sum2,Koef: real;

begin

su1: =0;

for p: =1 to s do

su1: =su1+nc [p,1];

sum1: =su1/3;

su2: =0;

for p: =1 to s do

su2: =su2+nc [p,2];

sum2: =su2/3;

writeln (out2,'Nup=': 8,'Ndn=': 8);

writeln (out2,sum1: 8: 4,sum2: 8: 4);

Koef: = (sum2-sum1) / (sum1+sum2);

writeln (out2,'Koef=': 8);

writeln (out2,Koef: 8: 4);

end;

begin

Randomize;

Assign (inp1,'c: \80inp. txt'); Reset (inp1);

writeln; writeln ('Start, program ID=',ProgID: 2); isp_No: =0;

repeat

Assign (out2,'c: \80all. txt'); Append (out2);

isp_No: =isp_No+1;

Read (inp1,n,L,rA,r01,r02,ugol,kpr,r1,r2,v0,ugol_x,ugol_y,R,Tm,Tmm,T,ac,m1);

writeln ('---------------------------Исходные-данные------------------№-в-серии: ', isp_No: 2);

writeln ('n=': 8,'L=': 8,'rA=': 8,'r01=': 8,'r02=': 8,'ugol=': 8,'kpr=': 8,'r1=': 8,'r2=': 8);

writeln (n: 8,L: 8: 2,rA: 8: 3,r01: 8: 3,r02: 8: 3,ugol: 8: 1,kpr: 8: 3,r1: 8: 3,r2: 8: 3);

writeln ('Vo=': 8,'ug_x=': 8,'ug_y=': 8,'Rg=': 8,'Tm=': 8,'Tmm=': 8,'T=': 8,'ac=': 8);

writeln (v0: 8: 0,ugol_x: 8: 1,ugol_y: 8: 1,R: 8: 1,Tm: 8: 1,Tmm: 8: 1,T: 8: 1,ac: 8: 3);

w0: =v0/sqrt (R*T); vTmm: =sqrt (Tmm/T); Fdir: =n;

vTm: =sqrt (Tm/T); vn: =0; Tsr: =0; Tss: =0; Vsr: =0; VVsr: =0; n1: =0; n2: =0; NNN: =0;

vx00: =w0*cos (ugol_x/180*pi); vyz00: =w0*sin (ugol_x/180*pi);

vy00: =vyz00*cos (ugol_y/180*pi); vz00: =vyz00*sin (ugol_y/180*pi);

vxTmm: =vTmm*cos (ugol_x/180*pi); vyzTmm: =vTmm*sin (ugol_x/180*pi);

vyTmm: =vyzTmm*cos (ugol_y/180*pi); vzTmm: =vyzTmm*sin (ugol_y/180*pi);

a: = (4*r12*r12 +2*r11*r12) *ugol/360;

q: = (r2*r2-r1*r1) *L/ps/ (2*ks) /a;

q1: = (r2*r2-r1*r1) /ps/ (2*ks) /a* (1-rA*rA); tekp: =0; {q1= (q/L) * (1-rA*rA), (1-rA^2) - площадь межэлектродной части торцев}

for k: =1 to 2*ks do for p: =1 to ps do nk1 [k,p]: =0;

for i: =1 to n do

begin

if ( (i/n*100) > (tekp+1)) then begin write (chr (8),chr (8),chr (8),chr (8),'%',tekp: 3); tekp: =tekp+1 end;

Maxwell; defl: =false;

vxTmmv: =v*vxTmm; vyTmmv: =v*vyTmm; vzTmmv: =v*vzTmm;

Maxwell; v: =v*vTm+0.0000000001;

{teta: =pi*random; } teta: =2*random-1;

{vx: =v*cos (teta) +vx00+vxTmmv; } vx: =v*teta+vx00+vxTmmv;

if not (vx=0) then begin

if (vx>0) then begin sk: =-1; x2: =0; n1: =n1+1 end

else begin sk: =1; x2: =L; n2: =n2+1 end;

repeat y2: =2*Random-1; z2: =2*Random-1;

r0: =sqrt (y2*y2+z2*z2);

byben;

{xkpr: =Random; }

until (y2*y2+z2*z2<r02*r02) and (y2*y2+z2*z2>r01*r01) and (m3=1) and ( ( (sk=1) and (-cos (pi*ugol/360) +y2/sqrt (y2*y2+z2*z2) >0))

or ( (sk=-1) and (cos (pi-pi*ugol/360) - y2/sqrt (y2*y2+z2*z2) >0)));

{vyz: =v*sin (teta); } vyz: =v*sqrt (1-teta*teta);

ksi: =2*pi*random;

vy: =vyz*cos (ksi) +vy00+vyTmmv; vz: =vyz*sin (ksi) +vz00+vzTmmv;

vv: =sqrt (vx*vx+vy*vy+vz*vz); VVsr: =VVsr+vv;

vn: =vn+abs (vx); j: =0; tid: =abs (L/vx); {< - (!) для полн. откр. трубы, оператор присв. для tid только здесь! }

for k: =1 to 2*ks do for p: =1 to ps do nk [k,p]: =0; Tsss: =0;

repeat counttrack;

r0: =sqrt (y2*y2+z2*z2);

xkpr: =Random;

if (xkpr<kpr) and (abs (f3) =1) and (y2*y2+z2*z2>r01*r01) and (y2*y2+z2*z2<r02*r02) and (m3=1) and ( ( (f3=-1)

and (-cos (pi*ugol/360) +y2/sqrt (y2*y2+z2*z2) >0))

or ( (f3=1) and (cos (pi-pi*ugol/360) - y2/sqrt (y2*y2+z2*z2) >0)))

then goto ll {если попала в зону прозрачности, то вылет}

else begin scattering (f3,f4); defl: =true end {иначе, отражение}

until j=1;

ll: Vsr: =Vsr+vv; if (defl) then Fdir: =Fdir-1;

for k: =1 to 2*ks do for p: =1 to ps do nk [k,p]: =nk [k,p] /tid; { отношение времен реал. (для МИП) и идеал. (для открытой трубы) }

for k: =1 to 2*ks do for p: =1 to ps do nk1 [k,p]: =nk1 [k,p] +nk [k,p];

end;

end;

for k: =1 to 2*ks do for p: =1 to ps do NNN: =NNN+nk1 [k,p]; { cуммируем все отношения времен реал. и идеал. по ч. объемам}

ww0: =vn/n; Tss: =Tss*a/n* (1-rA*rA); NNN: =NNN/n/q1/ps/ (2*ks); conc: =Tsr/n*ww0/q/ps/ (2*ks); Fdir: =Fdir/n;

{conc - старое знач. отн. конц. в МИП; в Tss до сих пор накапливалось отношение taukon/tid, теперь Tss отнормир.; нормируем NNN}

write (chr (8),chr (8),chr (8),chr (8));

writeln ('------------------------------Результаты--------------------------------');

writeln ('N_old=': 8,'N_new=': 8,'Tsr=': 8,'n1=': 8,'n2=': 8,'|Vx|=': 8,'|Vвх|=': 8,'|Vвых|=': 8,'Tss=': 8,'Fdir=': 7);

writeln (conc: 8: 4,NNN: 8: 4,Tsr/n: 8: 4,n1: 8,n2: 8,ww0: 8: 3,vvsr/n: 8: 3,vsr/n: 8: 3,Tss: 8: 4,Fdir: 7: 3);

writeln ('------------------------------------------------------------------------');

writeln (out2,'ProgID=': 8,'n=': 8,'L=': 8,'rA=': 8,'r01=': 8,'r02=': 8,'ugol=': 8,'kpr=': 8,'r1=': 8,'r2=': 8);

writeln (out2,ProgID: 8,n: 8,L: 8: 2,rA: 8: 3,r01: 8: 3,r02: 8: 3,ugol: 8: 1,kpr: 8: 3,r1: 8: 3,r2: 8: 3);

writeln (out2,'V0=': 8,'ug_x=': 8,'ug_y=': 8,'Rg=': 8,'Tm=': 8,'Tmm=': 8,'T=': 8,'ac=': 8);

writeln (out2,v0: 8: 0,ugol_x: 8: 1,ugol_y: 8: 1,R: 8: 1,Tm: 8: 1,Tmm: 8: 1,T: 8: 1,ac: 8: 3);

writeln (out2,'N_old=': 8,'N_new=': 8,'Tsr=': 8,'n1=': 8,'n2=': 8,'|Vx|=': 8,'|Vin|=': 8,'|Vout|=': 8,'Tss=': 8,'Fdir=': 8);

writeln (out2,conc: 8: 4,NNN: 8: 4,Tsr/n: 8: 4,n1: 8,n2: 8,ww0: 8: 3,vvsr/n: 8: 3,vsr/n: 8: 3,Tss: 8: 4,Fdir: 8: 3);

writeln (out2,'Nup=': 8,'Ndn=': 8,'D=': 8);

Nup: =0; Ndn: =0;

for p: =1 to ps do

begin

nc [p,1]: =0; nc [p,2]: =0; d [p]: =0;

for kR: =1 to 2*ks do

begin

if kR<=ks/2 then nc [p,1]: =nc [p,1] +nk1 [kR,p];

if (kR>ks/2) and (kR<=1.5*ks) then nc [p,2]: =nc [p,2] +nk1 [kR,p];

if kR>1.5*ks then nc [p,1]: =nc [p,1] +nk1 [kR,p];

if kR<=ks then d [p]: =d [p] +nk1 [kR,p];

if kR>ks then d [p]: =d [p] - nk1 [kR,p];

end;

d [p]: =d [p] /n/q1/ks;

nc [p,1]: =nc [p,1] /n/q1/ks; nc [p,2]: =nc [p,2] /n/q1/ks;

write (out2,nc [p,1]: 8: 4,nc [p,2]: 8: 4,d [p]: 8: 4); writeln (out2);

Nup: =Nup+nc [p,1]; Ndn: =Ndn+nc [p,2];

end;

{ vvod (Nsu); }

Nup: =Nup/ps; Ndn: =Ndn/ps;

writeln (out2,'Conc': 8,'priv': 8,'volms': 8);

for p: =1 to ps do

begin

write (out2, (nk1 [1,p] +nk1 [2*ks,p]) /2/n/q1: 8: 3);

for kR: =2 to 2*ks do write (out2, (nk1 [kR,p] +nk1 [kR-1,p]) /2/n/q1: 8: 3);

writeln (out2);

end;

summa;

writeln (out2);

close (out2);

Assign (out3,'c: \80int. txt'); Append (out3);

write (out3,ww0: 9: 4,vvsr/n: 9: 4,vsr/n: 9: 4,Tsr/n: 9: 4,Tss: 9: 4,Fdir: 9: 4,conc: 9: 4,NNN: 9: 4,Ndn: 9: 4,Nup: 9: 4,Ndn-Nup: 9: 4);

write (out3,nc [2,2] - nc [2,1]: 9: 4); {changeble addon}

for p: =1 to 2 do

for kL: =1 to ps do write (out3,nc [kL,p]: 9: 4);

writeln (out3);

close (out3);

until m1=0;

writeln ('Ready'); close (inp1);

end.

Приложение 2

Текст программы modmd82krug

{MODMD82. PAS - реализация метода АБ}

{Новые переменные:

q1 - новое значение вспомог. к-та q;

tid - (по идее) время пролета в зеркальноотражающем объеме;

Tsss - старое значение Tsss, теперь Tss выч. по др. формуле;

NNN - относительная концентрация в объеме по методу АБ;

nk1 [] - распределение отн. конц. по методу АБ. }

{Изменения в процедурах: только одна строчка в конце counttrack}

{16.04.2008: Мелкие несущественные изменения (интерфейс) }

program dat (input,output);

label ll; const ks=12; ps=3; ProgID=82; r21=0.600; r22=0.240;

var a,q,vx,vy,vz,vyz,v,vT,vTm,vTmm,v0,w0,vv,y,y0,teta,xkpr,ksi,vn,ww0,q1: real;

x2,y2,z2,vx00,vyz00,vz00,vy00,Tm,Tmm,r0,r1,r2,ugol,ac,ugol_x,ugol_y,Tsr,Tss,Vsr,VVsr: real;

vxTmm,vxTmmv,vyzTmm,vyTmm,vyTmmv,vzTmm,vzTmmv,tid,Tsss,NNN, R,T,L,rA,r01,r02,kpr: real;

j,p,m1,k,kR,kL,f3,f4,sk: integer; n, i,n1,n2: longint; isp_No,m3,tekp: integer; defl: boolean;

inp1,out2,out3: text; Fdir,conc,Nup,Ndn: real;

nk: array [1.2*ks,1. ps] of real;

nk1: array [1.2*ks,1. ps] of real; nc: array [1. ps,1.2] of real; d: array [1. ps] of real;

procedure Maxwell;

begin

repeat

v: = (random (5) +random);

y: =v*v*exp (1-v*v);

y0: =random;

until y0<=y;

if (v<=0.05) then v: =0.05;

end;

procedure scattering (Mx,My: real);

var vx0,vyz0,vy0,vz0,gamma,beta,alfa: real;

begin

repeat

beta: =pi/2*random;

y: =2*sin (beta) *cos (beta);

y0: =random;

until y0<=y;

Maxwell;

vv: = (1-ac) *vv*vv+ac*v*v+ac* (1-ac) *vv*vv* ( (v*v) /1.5-1);

if (vv<=0.0025) then vv: =0.05

else vv: =sqrt (vv);

alfa: =2*pi*random;

vx0: =vv*cos (beta); vyz0: =vv*sin (beta);

vy0: =vyz0*cos (alfa); vz0: =vyz0*sin (alfa);

if Mx<2 then begin vx: =Mx*vx0; vy: =vy0; vz: =Mx*vz0; end

else if not (y2*z2=0) then begin

if (y2*My) >0 then gamma: =arctan (z2/y2) +Pi-arctan (vy0/vx0)

else gamma: =arctan (z2/y2) - arctan (vy0/vx0);

vx: =-vz0; vy: =sqrt (vx0*vx0+vy0*vy0) *cos (gamma);

vz: =sqrt (vx0*vx0+vy0*vy0) *sin (gamma); end

else if y2=0 then begin

vx: =-vz0; vy: =-vy0*My*z2/abs (z2); vz: =-vx0*My*z2/abs (z2); end

else begin vx: =-vz0; vy: =-vx0*My*y2/abs (y2); vz: =vy0*My*y2/abs (y2); end;

end;

procedure counttrack;

label 2,3;

var x1,y1,z1,tmin1,tmax1,tmin2,tmax2,tmin,tmax,taumin,taumax,taukon,tAnod: real;

f0,f1,f2,fr1,fr2,fa2: real;

procedure countmolec;

var w: 0. ps; b1,b2,b3,b4,a1,a2: boolean;

d1,d2,d3,d4,tau1,tau2,tau3,tau4,tkon1,tkon2,t0,t00: real; label ii;

begin

Tsr: =Tsr+abs (tmax-tmin);

for w: =0 to ps-1 do

begin

tkon1: = ( (L*w/ps) - x1) /vx; tkon2: = ( (L* (w+1) /ps) - x1) /vx;

p: =w+1;

a1: = ( (tkon1>=tmin) and (tkon1<=tmax));

a2: = ( (tkon2>=tmin) and (tkon2<=tmax));

if a1 then if a2 then if vx>0 then begin t0: =tkon1; t00: =tkon2; end

else begin t0: =tkon2; t00: =tkon1; end

else if vx>0 then begin t0: =tkon1; t00: =tmax; end

else begin t0: =tmin; t00: =tkon1; end

else if a2 then if vx>0 then begin t0: =tmin; t00: =tkon2; end

else begin t0: =tkon2; t00: =tmax; end

else if (tkon1<tmin) and (tkon2>tmax) or (tkon1>tmax) and (tkon2<tmin)

then begin t0: =tmin; t00: =tmax; end

else goto ii;

d1: =y1-z1*vy/vz; d2: =d1;

tau1: = (d1-y1) /vy; tau2: = (d2-y1) /vy;

for kR: =1 to ks do

begin

d3: =vz*y1-vy*z1;

d4: =d3/ (vz*cos (kR*Pi/ks) +vy*sin (kR*Pi/ks));

d3: =d3/ (vz*cos (kR*Pi/ks) - vy*sin (kR*Pi/ks));

tau3: = (d3*cos (kR*Pi/ks) - y1) /vy; tau4: = (d4*cos (kR*Pi/ks) - y1) /vy;

b1: = ( (d1>0) and (tau1<=t00) and (tau1>=t0));

b2: = ( (d2>0) and (tau2<=t00) and (tau2>=t0));

b3: = ( (d3>0) and (tau3<=t00) and (tau3>=t0));

b4: = ( (d4>0) and (tau4<=t00) and (tau4>=t0));

if b1 then if b3 then nk [kR,p]: =nk [kR,p] +abs (tau1-tau3)

else if ( (tau3-tau1) *d3) >0 then nk [kR,p]: =nk [kR,p] +abs (tau1-t00)

else nk [kR,p]: =nk [kR,p] +abs (tau1-t0)

else if b3 then if ( (tau1-tau3) *d1) >0 then nk [kR,p]: =nk [kR,p] +abs (tau3-t00)

else nk [kR,p]: =nk [kR,p] +abs (tau3-t0)