4. Математические методы решения

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

Пусть требуется найти решение дифференциального уравнения первого порядка

,

удовлетворяющее начальному условию

.

Численный метод Эйлера дает возможность в некоторых точках x₁, x₂, … , x_n найти приближения y₁,y₂, …, y_nдля значений точного решенияy(x_k), гдеk=1,2, …,n. Обычно полагаютx_k=x₀+kh, гдеk=1,2, … ,n;h0 – шаг. Вычислительная формула метода Эйлера имеет вид:

,k=0,1, … , n.

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

Первая схема (Усовершенствованный метод Эйлера)

Вычисления на каждом шаге kвыполняются в следующем порядке:

, , .

Вторая схема (усовершенствованный метод Эйлера – Коши)

Сначала, по основной вычислительногй формуле метода Эйлера вычисляется грубое значение y⁽⁰⁾_k+1:

Затем вычисленное значение уточняется по формуле:

.

Погрешность этого метода имеет порядок h³.

Метод Эйлера и его модификации фактически являются элементарными схемами метода Рунге-Кутта.

5. Программирование

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

6. Оформление домашнего задания

Домашнее задание должно быть оформлено по единой схеме, включающей следующие материалы:

1. титульный лист;

2. введение;

3. расчетная схема;

4. математическая модель;

5. программа расчета на ЭВМ (исходный текст);

6. краткая инструкция по эксплуатации программы;

7. графическое изображение результатов расчета;

8. анализ полученных результатов.

По всем трем этапам выполнения домашнего задания выпускается единый отчет. Зачет по домашнему заданию принимается индивидуально с демонстрацией работы программы на ЭВМ. В том случае, если программа была написана на персональном компьютере вне кафедры, то к зачету, в том числе предъявляется дискета с загрузочным модулем.

7. Пример

В качестве примера выполнения домашнего задания составим математическую модель и программу расчета на ЭВМ для пневматической системы, расчетная схема которой показана на рис.2.

Расчет выполним в три этапа:

0. без учета теплообмена между рабочим телом и стенками емкостей пневмосистемы;

1. с учетом теплообмена между рабочим телом и стенками емкостей пневмосистемы;

2. с учетом заправки емкостей пневмосистемы раздельно воздухом и гелием, с последующим образованием смеси газов.

Этап 1Рис.2

МАТЕМАТИЧЕСКАЯ МОДЕЛЬ

Исходные данные:

V₁; V₂; V₃; V₄; d₁₃; d₂₃; d₃₄;

₁₃;₂₃;₃₄

R = 287.2 Дж/(кг×К); k=1.41;

_кр=0.528;π=3.14;=6 Вт/(м²×К)

Начальные условия:

t=0;

Р₁=Р_1Н; Р₂=Р_2Н; Р₃=Р_3Н; Р₄=Р_4Н;

Т₁=300 К; Т₂=300 К; Т₃=300 К; Т₄=300 К.

Предварительные расчеты:

F₁₃=₁₃d₁₃²/4

F₂₃=₂₃d₂₃²/4

F₃₄=₃₄d₃₄²/4

S₁= (6V₁/)^2/3

S₂= (6V₂/)^2/3

S₃=(6V₃/)^2/3

S₄=(6V₄/)^2/3

Уравнения расхода:

G₁₃=G(F₁₃,P₁,P₃,T_Pmax1)

G₂₃=G(F₂₃,P₂,P₃,T_Pmax2)

G₃₄=G(F₃₄,P₃,P₄,T_Pmax3),

где T_Pmaxi – температура газа в той из двух рассматриваемых емкостей, где давление больше.

Уравнения расчета изменения параметров в емкостях ПС:

1-я емкость:

2-я емкость

3-я емкость

4-я емкость

Конечное условие: P_{min }0.95*P_max

ПРОГРАММА РАСЧЕТА НА ЭВМ

Program Raschet_PS1;

uses CRT;

const

R=287.2; K=1.41; BK=0.528; PI=3.14; TST=300;

Var

ALF,M1,M2,M3,D1,D2,D3,V1,V2,V3,V4,F1,F2,F3,PA1,PA2,PA3,PB1,PB2,PB3, TA1,TA2,TA3,AA,BB,CC,DD,EE,FF,G1,G2,G3,DT,PM,PX,H,tt:real;

i,l,m,n:integer;

Y,V:array [1..8] of real;

BEGIN

clrscr;

writeln('Vvod ishodnih dannih');

writeln;

writeln('Zadaite ALF,M1,M2,M3,D1,D2,D3,V1,V2,V3,V4,DT'); readln(ALF,M1,M2,M3,D1,D2,D3,V1,V2,V3,V4,DT);

writeln;

writeln('Zadaite P1,T1,P2,T2,P3,T3,P4,T4');

for i:=1 to 8 do

read(Y[i]);

writeln;

writeln('Kontrol vvedennih dannih');

writeln('V1=',V1:4:4,' ','V2=',V2:4:4,' ','V3=',V3:4:4,' ','V4=',V4:4:4);

writeln('P1=',Y[1]:3:2,' ','P2=',Y[3]:3:2,' ','P3=',Y[5]:3:2,' ','P4=',Y[7]:3:2);

writeln('T1=',Y[2]:3:2,' ','T2=',Y[4]:3:2,' ','T3=',Y[6]:3:2,' ','T4=',Y[8]:3:2);

writeln;

writeln;

writeln('Raschet');

writeln;

writeln('Secheniya truboprovodov');

writeln;

F1:=M1*PI*sqr(D1)/4;

F2:=M2*PI*sqr(D2)/4;

F3:=M3*PI*sqr(D3)/4;

S1:=PI*exp((2/3)*ln(6*V1/PI));

S2:=PI*exp((2/3)*ln(6*V2/PI));

S3:=PI*exp((2/3)*ln(6*V3/PI));

S4:=PI*exp((2/3)*ln(6*V4/PI));

writeln('F1=',F1:4:4,' ','F2=',F2:4:4,' ','F3=',F3:4:4);

writeln('S1=',S1:4:4,' ','S2=',S2:4:4,' ','S3=',S3:4:4,' ','S4=',S4:4:4);

writeln;

tt:=0;

n:=1;

REPEAT

writeln('Opredelenie napravleniya rashoda');

writeln;

if Y[1]>Y[5] then

begin

PA1:=Y[1]; PB1:=Y[5]; TA1:=Y[2]; AA:=-1; BB:=1

end

else begin

PA1:=Y[5]; PB1:=Y[1]; TA1:=Y[6]; AA:=1; BB:=-1

end;

writeln('PA1=',PA1:3:2,' ','PB1=',PB1:3:2,' ','TA1=',TA1:3:2,' ','AA=',AA:3:2,' ','BB=',BB:3:2);

writeln;

if Y[3]>Y[5] then

begin

PA2:=Y[3]; PB2:=Y[5]; TA2:=Y[4]; CC:=-1; DD:=1

end

else begin

PA2:=Y[5]; PB2:=Y[3]; TA2:=Y[6]; CC:=1; DD:=-1

end;

writeln('PA2=',PA2:3:2,' ','PB2=',PB2:3:2,' ','TA2=',TA2:3:2,' ','CC=',CC:3:2,' ','DD=',DD:3:2);

writeln;

if Y[5]>Y[7] then

begin

PA3:=Y[5]; PB3:=Y[7]; TA3:=Y[6]; EE:=-1; FF:=1

end

else begin

PA3:=Y[7]; PB3:=Y[5]; TA3:=Y[8]; EE:=1; FF:=-1

end;

writeln('PA3=',PA3:3:2,' ','PB3=',PB3:3:2,' ','TA3=',TA3:3:2,' ','EE=',EE:3:2,' ','FF=',FF:3:2);

writeln;

writeln('Raschet rashodov');

writeln;

if (PB1/PA1)>BK then

begin

H:= 2*K*(exp((2/K)*ln(PB1/PA1))-exp(((K+1)/K)*ln(PB1/PA1)))/(K-1)/R/TA1;

if H<0 then G1:=0

else

G1:=F1*PA1*sqrt(H);

end

else begin

H:= 2*K*(exp((2/K)*ln(BK))-exp(((K+1)/K)*ln(BK)))/(K-1)/R/TA1;

if H<0 then G1:=0

else

G1:=F1*PA1*sqrt(H);

end;

if (PB2/PA2)>BK then

begin

H:= 2*K*(exp((2/K)*ln(PB2/PA2))-exp(((K+1)/K)*ln(PB2/PA2)))/(K-1)/R/TA2;

if H<0 then G2:=0

else

G2:=F2*PA2*sqrt(H)

end

else begin

H:= 2*K*(exp((2/K)*ln(BK))-exp(((K+1)/K)*ln(BK)))/(K-1)/R/TA2;

if H<0 then G2:=0

else

G2:=F2*PA2*sqrt(H);

end;

if (PB3/PA3)>BK then

begin

H:= 2*K*(exp((2/K)*ln(PB3/PA3))-exp(((K+1)/K)*ln(PB3/PA3)))/(K-1)/R/TA3;

if H<0 then G3:=0

else

G3:=F3*PA3*sqrt(H)

end

else begin

H:= 2*K*(exp((2/K)*ln(BK))-exp(((K+1)/K)*ln(BK)))/(K-1)/R/TA3;

if H<0 then G3:=0

else

G3:=F3*PA3*sqrt(H);

end;

writeln(‘Shag’, n:4,' ', 't=',tt:3:3);

writeln;

writeln('P1=',Y[1]:3:2,' ','P2=',Y[3]:3:2,' ','P3=',Y[5]:3:2,' ','P4=',Y[7]:3:2);

writeln('T1=',Y[2]:3:2,' ','T2=',Y[4]:3:2,' ','T3=',Y[6]:3:2,' ','T4=',Y[8]:3:2);

writeln('G1=',G1:6:6,' ','G2=',G2:6:6,' ','G3=',G3:6:6);

writeln;

n:=n+1;

tt:=tt+DT;

V[1]:=K*R*TA1*G1*AA/V1;

V[2]:=Y[2]*(V1*V[1]-R*Y[2]*G1*AA)/(Y[1]*V1);

V[3]:=K*R*TA2*G2*CC/V2;

V[4]:=Y[4]*(V2*V[3]-R*Y[4]*G2*CC)/(Y[3]*V2);

V[5]:=K*(R*TA1*G1*BB+R*TA2*G2*DD+R*TA3*G3*EE)/V3;

V[6]:=Y[6]*(V3*V[5]-R*Y[6]*(G1*BB+G2*DD+G3*EE))/(Y[5]*V3);

V[7]:=K*R*TA3*G3*FF/V4;

V[8]:=Y[8]*(V4*V[7]-R*Y[8]*G3*FF)/(Y[7]*V4);

for i:=1 to 8 do

Y[i]:=Y[i]+V[i]*DT;

PX=MAX (Y(1),Y(3),Y(5), Y(7))

PM=MIN(Y(1),Y(3),Y(5),Y(7))

UNTILPM>=0.95*PX

END.

Результаты расчета ПС по данной программе в графическом виде приведены на рис.3, 4,5.

Расчеты проведены для следующих начальных условий: Р₁=2·10⁵Па; Р₂=1·10⁵Па; Р₃=6·10⁵Па; Р₄=12·10⁵Па;

Т₁=300 К; Т₂=300 К; Т₃=300 К; Т₄=300 К.

Рис.3

Рис.4

Рис.5

Этап 2

МАТЕМАТИЧЕСКАЯ МОДЕЛЬ

Исходные данные:

V₁; V₂; V₃; V₄; d₁₃; d₂₃; d₃₄;

₁₃;₂₃;₃₄

R = 287.2 Дж/(кг×К); k=1.41;

_кр=0.528;π=3.14;=6 Вт/(м²×К)

Начальные условия:

t=0;

Р₁=Р_1Н; Р₂=Р_2Н; Р₃=Р_3Н; Р₄=Р_4Н;

Т₁=300 К; Т₂=300 К; Т₃=300 К; Т₄=300 К.

Предварительные расчеты:

F₁₃=₁₃d₁₃²/4

F₂₃=₂₃d₂₃²/4

F₃₄=₃₄d₃₄²/4

S₁= (6V₁/)^2/3

S₂= (6V₂/)^2/3

S₃=(6V₃/)^2/3

S₄=(6V₄/)^2/3

Уравнения расхода:

G₁₃=G(F₁₃,P₁,P₃,T_Pmax1)

G₂₃=G(F₂₃,P₂,P₃,T_Pmax2)

G₃₄=G(F₃₄,P₃,P₄,T_Pmax3),

где T_Pmaxi – температура газа в той из двух рассматриваемых емкостей, где давление больше.

Уравнения расчета изменения параметров в емкостях ПС:

1-я емкость:

2-я емкость

3-я емкость

4-я емкость

Конечное условие: P_{min }0.95*P_max.

ПРОГРАММА РАСЧЕТА НА ЭВМ

Program Raschet_PS2;

uses CRT; const

R=287.2; K=1.41; BK=0.528; PI=3.14; TST=300;

var

ALF,M1,M2,M3,D1,D2,D3,V1,V2,V3,V4,F1,F2,F3,S1,S2,S3,S4,PA1,PA2,PA3,PB1,PB2,PB3, TA1,TA2,TA3,AA,BB,CC,DD,EE,FF,G1,G2,G3,DT,DQ1,DQ2,DQ3,DQ4,PM,PX,H,tt:real;

i,l,m,n:integer;

Y,V:array [1..8] of real;

BEGIN

clrscr;

writeln('Vvod ishodnih dannih');

writeln;

writeln('Zadaite ALF,M1,M2,M3,D1,D2,D3,V1,V2,V3,V4,DT'); readln(ALF,M1,M2,M3,D1,D2,D3,V1,V2,V3,V4,DT);

writeln;

writeln('Zadaite P1,T1,P2,T2,P3,T3,P4,T4');

for i:=1 to 8 do

read(Y[i]);

writeln;

writeln('Kontrol vvedennih dannih');

writeln('V1=',V1:4:4,' ','V2=',V2:4:4,' ','V3=',V3:4:4,' ','V4=',V4:4:4);

writeln('P1=',Y[1]:3:2,' ','P2=',Y[3]:3:2,' ','P3=',Y[5]:3:2,' ','P4=',Y[7]:3:2);

writeln('T1=',Y[2]:3:2,' ','T2=',Y[4]:3:2,' ','T3=',Y[6]:3:2,' ','T4=',Y[8]:3:2);

writeln;

writeln;

writeln('Raschet');

writeln;

writeln('Secheniya i poverhnosti');

writeln;

F1:=M1*PI*sqr(D1)/4;

F2:=M2*PI*sqr(D2)/4;

F3:=M3*PI*sqr(D3)/4;

S1:=PI*exp((2/3)*ln(6*V1/PI));

S2:=PI*exp((2/3)*ln(6*V2/PI));

S3:=PI*exp((2/3)*ln(6*V3/PI));

S4:=PI*exp((2/3)*ln(6*V4/PI));

writeln('F1=',F1:4:4,' ','F2=',F2:4:4,' ','F3=',F3:4:4);

writeln('S1=',S1:4:4,' ','S2=',S2:4:4,' ','S3=',S3:4:4,' ','S4=',S4:4:4);

writeln;

tt:=0;

n:=1;

REPEAT

writeln('Opredelenie napravleniya rashoda');

writeln;

if Y[1]>Y[5] then

begin

PA1:=Y[1]; PB1:=Y[5]; TA1:=Y[2]; AA:=-1; BB:=1

end

else begin

PA1:=Y[5]; PB1:=Y[1]; TA1:=Y[6]; AA:=1; BB:=-1

end;

writeln('PA1=',PA1:3:2,' ','PB1=',PB1:3:2,' ','TA1=',TA1:3:2,' ','AA=',AA:3:2,' ','BB=',BB:3:2);

writeln;

if Y[3]>Y[5] then

begin

PA2:=Y[3]; PB2:=Y[5]; TA2:=Y[4]; CC:=-1; DD:=1

end

else begin

PA2:=Y[5]; PB2:=Y[3]; TA2:=Y[6]; CC:=1; DD:=-1

end;

writeln('PA2=',PA2:3:2,' ','PB2=',PB2:3:2,' ','TA2=',TA2:3:2,' ','CC=',CC:3:2,' ','DD=',DD:3:2);

writeln;

if Y[5]>Y[7] then

begin

PA3:=Y[5]; PB3:=Y[7]; TA3:=Y[6]; EE:=-1; FF:=1

end

else begin

PA3:=Y[7]; PB3:=Y[5]; TA3:=Y[8]; EE:=1; FF:=-1

end;

writeln('PA3=',PA3:3:2,' ','PB3=',PB3:3:2,' ','TA3=',TA3:3:2,' ','EE=',EE:3:2,' ','FF=',FF:3:2);

writeln;

writeln('Raschet rashodov');

writeln;

if (PB1/PA1)>BK then

begin

H:= 2*K*(exp((2/K)*ln(PB1/PA1))-exp(((K+1)/K)*ln(PB1/PA1)))/(K-1)/R/TA1;

if H<0 then G1:=0

else

G1:=F1*PA1*sqrt(H);

end

else begin

H:= 2*K*(exp((2/K)*ln(BK))-exp(((K+1)/K)*ln(BK)))/(K-1)/R/TA1;

if H<0 then G1:=0

else

G1:=F1*PA1*sqrt(H);

end;

if (PB2/PA2)>BK then

begin

H:= 2*K*(exp((2/K)*ln(PB2/PA2))-exp(((K+1)/K)*ln(PB2/PA2)))/(K-1)/R/TA2;

if H<0 then G2:=0

else

G2:=F2*PA2*sqrt(H)

end

else begin

H:= 2*K*(exp((2/K)*ln(BK))-exp(((K+1)/K)*ln(BK)))/(K-1)/R/TA2;

if H<0 then G2:=0

else

G2:=F2*PA2*sqrt(H);

end;

if (PB3/PA3)>BK then

begin

H:= 2*K*(exp((2/K)*ln(PB3/PA3))-exp(((K+1)/K)*ln(PB3/PA3)))/(K-1)/R/TA3;

if H<0 then G3:=0

else

G3:=F3*PA3*sqrt(H)

end

else begin

H:= 2*K*(exp((2/K)*ln(BK))-exp(((K+1)/K)*ln(BK)))/(K-1)/R/TA3;

if H<0 then G3:=0

else

G3:=F3*PA3*sqrt(H);

end;

writeln(‘Shag’, n:4,' ', 't=',tt:3:3);

writeln;

writeln('P1=',Y[1]:3:2,' ','P2=',Y[3]:3:2,' ','P3=',Y[5]:3:2,' ','P4=',Y[7]:3:2);

writeln('T1=',Y[2]:3:2,' ','T2=',Y[4]:3:2,' ','T3=',Y[6]:3:2,' ','T4=',Y[8]:3:2);

writeln('G1=',G1:6:6,' ','G2=',G2:6:6,' ','G3=',G3:6:6);

writeln;

n:=n+1;

tt:=tt+DT;

l:=1; m:=1; i:=1;

repeat

i:=i+2;

if Y[i]>Y[l] then l:=i;

if Y[i]<Y[m] then m:=i;

until i=7;

PX:=Y[l];

PM:=Y[m];

writeln('Maksimalnoe davlenie PX=',PX:3:2)

writeln('Minimalnoe davlenie PM=',PM:3:2);

writeln('PM/PX=',PM/PX:3:2);

writeln;

DQ1:=ALF*S1*(TST-Y[2]);

V[1]:=K*(R*TA1*G1*AA+DQ1*(K-1))/V1;

V[2]:=Y[2]*(V1*V[1]-R*Y[2]*G1*AA+DQ1*(K-1))/(Y[1]*V1);

DQ2:=ALF*S2*(TST-Y[4]);

V[3]:=K*(R*TA2*G2*CC+DQ2*(K-1))/V2;

V[4]:=Y[4]*(V2*V[3]-R*Y[4]*G2*CC+DQ2*(K-1))/(Y[3]*V2);

DQ3:=ALF*S3*(TST-Y[6]);

V[5]:=K*(R*TA1*G1*BB+R*TA2*G2*DD+R*TA3*G3*EE+DQ3*(K-1))/V3;

V[6]:=Y[6]*(V3*V[5]-R*Y[6]*(G1*BB+G2*DD+G3*EE)+DQ3*(K-1))/(Y[5]*V3);

DQ4:=ALF*S4*(TST-Y[8]);

V[7]:=K*(R*TA3*G3*FF+DQ4*(K-1))/V4;

V[8]:=Y[8]*(V4*V[7]-R*Y[8]*G3*FF+DQ4*(K-1))/(Y[7]*V4);

for i:=1 to 8 do

Y[i]:=Y[i]+V[i]*DT;

UNTIL PM>=0.95*PX

END.

Результаты расчета ПС по данной программе в графическом виде приведены на рис.6, 7,8.

Расчеты проведены для следующих начальных условий: Р₁=2·10⁵Па; Р₂=1·10⁵Па; Р₃=6·10⁵Па; Р₄=12·10⁵Па;

Т₁=300 К; Т₂=300 К; Т₃=300 К; Т₄=300 К.

Рис.6

Рис.7

Рис.8

Этап 3

Расчетная схема представлена на рис.9. Емкость2 заправлена гелием, а емкости 1 и 3 -воздухом.

МАТЕМАТИЧЕСКАЯ МОДЕЛЬ

Исходные данные:

V₁; V₂; V₃; d₁₃; d₂₃;

₁₃;₂₃;

R_В= 287.2 Дж/(кг×К);k_В=1.41;_КРВ=0.528;

R_He= 2077.4 Дж/(кг×К);k_He=1.66;

π=3.14;=6 Вт/(м²×К)

Начальные условия:

t=0;

Р₁=Р_1Н; Р₂=Р_2Н; Р₃=Р_3Н;

Т₁=300 К; Т₂=300 К; Т₃=300 К. Рис.9

Предварительные расчеты:

F₁₃=₁₃d₁₃²/4

F₂₃=₂₃d₂₃²/4

S₁=(6V₁/)^2/3

S₂=(6V₂/)^2/3

S₃=(6V₃/)^2/3

Уравнения расхода:

G₁₃=G(F₁₃,P₁,P₃,T_Pmax1)

G₂₃=G(F₂₃,P₂,P₃,T_Pmax2),

где T_Pmaxi – температура газа в той из двух рассматриваемых емкостей, где давление больше.

Расчет концентрации, газовой постоянной и показателя изоэнтропы смеси:

концентрации гелия и воздуха:

газовая постоянная смеси:

показатель изоэнтропы смеси:

Уравнения расчета изменения параметров в емкостях ПС:

1-я емкость:

2-я емкость

3-я емкость

Конечное условие: P_{min }0.95*P_max

ПРОГРАММА РАСЧЕТА НА ЭВМ

Program Raschet_PS3; uses CRT;

const

R=287,2; K=1.41; R1=2077,4; K1=1.66; BK=0.528; PI=3.14; TST=300;

var

ALF,M1,M2,D1,D2,V1,V2,V3,F1,F2,S1,S2,S3,BK1,PA1,PA2,PB1,PB2,TA1,TA2,

AA,BB,CC,DD,G1,G2,DT,DQ1,DQ2,DQ3,PM,PX,H,GH,GV,gHe,KM,RM,tt:real;

i,l,m,n:integer;

Y,V:array [1..6] of real;

BEGIN

clrscr;

writeln('Vvod ishodnih dannih');

writeln;

writeln('Zadaite ALF,M1,M2,D1,D2,V1,V2,V3,DT');

readln(ALF,M1,M2,D1,D2,V1,V2,V3,DT);

writeln;

writeln('Zadaite P1,T1,P2,T2,P3,T3');

for i:=1 to 6 do

read(Y[i]);

writeln;

writeln('Kontrol vvedennih dannih');

writeln('V1=',V1:4:4,' ','V2=',V2:4:4,' ','V3=',V3:4:4);

writeln('P1=',Y[1]:3:2,' ','P2=',Y[3]:3:2,' ','P3=',Y[5]:3:2);

writeln('T1=',Y[2]:3:2,' ','T2=',Y[4]:3:2,' ','T3=',Y[6]:3:2);

writeln;

writeln;

writeln('Raschet');

writeln;

writeln('Secheniya i poverhnosti');

writeln;

F1:=M1*PI*sqr(D1)/4;

F2:=M2*PI*sqr(D2)/4;

S1:=PI*exp((2/3)*ln(6*V1/PI));

S2:=PI*exp((2/3)*ln(6*V2/PI));

S3:=PI*exp((2/3)*ln(6*V3/PI));

writeln('F1=',F1:4:4,' ','F2=',F2:4:4);

writeln('S1=',S1:4:4,' ','S2=',S2:4:4,' ','S3=',S3:4:4);

writeln;

writeln('Kriticheskii perepad davleniya dlya geliya');

writeln;

BK1:=exp((K1/(K1-1))*ln(2/(K1+1)));

writeln('BK1=',BK1:3:2);

writeln;

tt:=0;

n:=1;

GH:=0;

GV:=0;

REPEAT

writeln('Opredelenie napravleniya rashoda');

writeln;

if Y[1]>Y[5] then

begin

PA1:=Y[1]; PB1:=Y[5]; TA1:=Y[2]; AA:=-1; BB:=1

end

else begin

PA1:=Y[5]; PB1:=Y[1]; TA1:=Y[6]; AA:=1; BB:=- 1

end;

writeln('PA1=',PA1:3:2,'','PB1=',PB1:3:2,'','TA1=',TA1:3:2,'','AA=',AA:3:2,' ','BB=',BB:3:2);

writeln;

if Y[3]>Y[5] then

begin

PA2:=Y[3]; PB2:=Y[5]; TA2:=Y[4]; CC:=-1; DD:=1;

end

else begin

PA2:=Y[5]; PB2:=Y[3]; TA2:=Y[6]; CC:=1; DD:=-1;

end;

writeln('PA2=',PA2:3:2,'','PB2=',PB2:3:2,'','TA2=',TA2:3:2,' ','CC=',CC:3:2,' ','DD=',DD:3:2);

writeln;

writeln('Raschet rashodov');

writeln;

if (PB1/PA1)>BK then

begin

H:= 2*K*(exp((2/K)*ln(PB1/PA1))- exp(((K+1)/K)*ln(PB1/PA1))) /(K-1)/R/TA1;

if H<0 then G1:=0

else

G1:=F1*PA1*sqrt(H);

end

else begin

H:= 2*K*(exp((2/K)*ln(BK))-exp(((K+1)/K)*ln(BK)))/(K-1)/R/TA1;

if H<0 then G1:=0

else

G1:=F1*PA1*sqrt(H);

end;

if (PB2/PA2)>BK1 then

begin

H:= 2*K1*(exp((2/K1)*ln(PB2/PA2))-exp(((K1+1)/K1)*ln(PB2/PA2)))/(K1-1)/R1/TA2;

if H<0 then G2:=0

else

G2:=F2*PA2*sqrt(H)

end

else begin

H:= 2*K1*(exp((2/K1)*ln(BK1))-exp(((K1+1)/K1)*ln(BK1)))/(K1-1)/R1/TA2;

if H<0 then G2:=0

else

G2:=F2*PA2*sqrt(H);

end;

GV:=GV+G1*DT;

GH:=GH+G2*DT;

gHe:=GH/(GH+GV);

gV:=GV/(GH+GV);

KM:=K1*gHe+K*(1-gHe);

RM:=R1*gHe+R*(1-gHe);

writeln('Shag',n:4,' ', 't=',tt:3:3);

writeln;

writeln('P1=',Y[1]:3:2,'','P2=',Y[3]:3:2,' ','P3=',Y[5]:3:2);

writeln('T1=',Y[2]:3:2,'','T2=',Y[4]:3:2,' ','T3=',Y[6]:3:2);

writeln('G1=',G1:6:6,' ','G2=',G2:6:6);

writeln('Koncentraciya geliya v kamere smesheniya gHe=',gHe:5:5);

writeln('Koncentraciya vozduha v kamere smesheniya gV=',gV:5:5);

writeln('Pokazatel adiabati smesi gazov KM=',KM:3:2);

writeln('Universalnaya gazovaya postoyannaya smesi RM=',RM:3:2);

writeln;

n:=n+1;

tt:=tt+DT;

l:=1; m:=1; i:=1;

repeat

i:=i+2;

if Y[i]>Y[l] then l:=i;

if Y[i]<Y[m] then m:=i;

until i=5;

PX:=Y[l];

PM:=Y[m];

writeln('Maksimalnoe davlenie PX=',PX:3:2);

writeln('Minimalnoe davlenie PM=',PM:3:2);

writeln;

writeln('PM/PX=',PM/PX:3:2);

writeln;

writeln;

DQ1:=ALF*S1*(TST-Y[2]);

V[1]:=K*(R*TA1*G1*AA+DQ1*(K-1))/V1;

V[2]:=Y[2]*(V1*V[1]-R*Y[2]*G1*AA+DQ1*(K-1))/(Y[1]*V1);

DQ2:=ALF*S2*(TST-Y[4]);

V[3]:=K1*(R1*TA2*G2*CC+DQ2*(K1-1))/V2;

V[4]:=Y[4]*(V2*V[3]-R1*Y[4]*G2*CC+DQ2*(K1-1))/(Y[3]*V2);

DQ3:=ALF*S3*(TST-Y[6]);

V[5]:=KM*(RM*TA1*G1*BB+RM*TA2*G2*DD+DQ3*(KM-1))/V3;

V[6]:=Y[6]*(V3*V[5]-RM*Y[6]*(G1*BB+G2*DD)+DQ3*(KM-1))/(Y[5]*V3);

for i:=1 to 6 do

Y[i]:=Y[i]+V[i]*DT;

UNTIL PM>=0.95*PX

END.

Результаты расчета ПС по данной программе в графическом виде приведены на рис.10,11,12,13,14

Расчеты проведены для следующих начальных условий: Р₁=12·10⁵Па; Р₂=11,5·10⁵Па; Р₃=0,01·10⁵Па;

Т₁=300 К; Т₂=300 К; Т₃=300 К.

Рис.10

Рис.11

Рис.12

Рис.13

Рис.14

анализ полученных результатов

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

СОДЕРЖАНИЕ

Введение…………………………………………………………………………………………………...1