Směs dvou ideálních plynů

Quit  [];

kB=1.380649 (10^(-23));(*J/K*)
h=6.62607015 (10^(-34)); (*J/Hz*)
Da=1.66053906660 (10^(-27));(*kg*)

m2=10 Da;

(*Celkové stavové proměnné*)
ρ=1; (*kg/(m^3)*)
u=0 (10^3); (*kg m/s*) 
e= (3/2) kB 300((ρ/m2))+(1/2)((u^2)/ρ); (*J/(m^3)*)
(*T=300;(*K*) *)

(*Parametry výpočtu*)
κnum=100;κmin=0.5;κmax=2;Ratios4κ=Table[{},{i,1,κnum+1}];

For[j=1,j≤κnum+1,j++,
κ=κmin+(j-1) (κmax-κmin)/κnum; 
rel26=((κ )^(5/2 κ))(m2^(4 (1-κ))) ((ρ1)^κ)(((4 (π/3) ((e-((u^2)/(2 ρ)))/((h^2) (κ ρ1+ρ-ρ1)))))^(3/2 (1-κ)));
rel24=((2 e (ρ^2)+(u^2) ((κ^(-1))-1)(ρ-ρ1))/(2 ρ (ρ1+(κ^(-1))(ρ-ρ1))))((ρ1)/ρ);
rel23=((ρ1)/ρ) u;
ρ1sol=FindRoot[ρ-ρ1==rel26,{ρ1,0.5 ρ}]//Flatten;
e1sol=FindRoot[e1==rel24/.ρ1sol,{e1,0.5 e}]//Flatten;
u1sol=FindRoot[u1==((ρ1)/ρ)u/.ρ1sol,{u1,0.5 u}]//Flatten;
ρ2sol={ρ2->((ρ-ρ1)/.ρ1sol)};
e2sol={e2->(e-e1/.e1sol)};
u2sol={u2->(((ρ2)/ρ) u/.ρ2sol)};
s1sol={s1->(kB ((ρ1)/m1)((5/2)+Log[(m1/(ρ1))(((((4 π m1)/(3 (h^2))) ((e1-((u1^2)/(2 ρ1)))/(ρ1/m1))))^(3/2))])/.m1->m2 (κ^(-1))/.ρ1sol/.u1sol/.e1sol)};
s2sol={s2->(kB ((ρ2)/m2)((5/2)+Log[(m2/(ρ2))(((((4 π m2)/(3 (h^2))) ((e2-((u2^2)/(2 ρ2)))/(ρ2/m2))))^(3/2))])/.ρ2sol/.u2sol/.e2sol)};
ρstar=-(kB/m1)Log[((ρ1)/m1)]+(3/2)(kB/m1)Log[((4 π m1)/(3 (h^2))) (m1/(ρ1))(e1-((u1^2)/(2 ρ1)))]+(3/2)(kB/m1) (1/(e1-((u1^2)/(2 ρ1))))((u1^2)/(2 ρ1))/.m1->m2 (κ^(-1))/.ρ1sol/.u1sol/.e1sol;
estar=D[(kB ((ρ1)/m1)((5/2)+Log[(m1/(ρ1))(((((4 π m1)/(3 (h^2))) ((e1-((u1^2)/(2 ρ1)))/(ρ1/m1))))^(3/2))])),e1]/.m1->m2 (κ^(-1))/.ρ1sol/.u1sol/.e1sol;
p=(-e+ρ μ+T s)/.μ->(-ρstar T)/.T->1/estar/.s->s1+s2/.ρ1sol/.e1sol/.s1sol/.s2sol;
p1=-e1+ρ1 μ+T s1/.μ->(-ρstar T)/.T->1/estar/.e1sol/.ρ1sol/.s1sol;
p2=-e2+ρ2 μ+T s2/.μ->(-ρstar T)/.T->1/estar/.e2sol/.ρ2sol/.s2sol;
Ratios4κ[[j]]={κ,Chop[ρ1/ρ/.ρ1sol],Chop[(((ρ1)/ρ)e/e1)/.ρ1sol/.e1sol],
Chop[(((ρ1)/ρ)s/s1)/.s->s1+s2/.ρ1sol/.s1sol/.s2sol],
Chop[((((ρ1)/ρ)e-e1)-T(((ρ1)/ρ)s-s1))/(((ρ1)/ρ)p)/.s->s1+s2/.ρ1sol/.e1sol/.s1sol/.s2sol],Chop[(ρ1)/.ρ1sol],Chop[(ρ2)/.ρ2sol],Chop[(e1)/.e1sol],Chop[(e2)/.e2sol],Chop[(s1)/.s1sol],Chop[(s2)/.s2sol],Chop[p],Chop[p1],Chop[p2],Chop[estar]}];

Celkový tlak a parciální tlaky

Clear[κ,j,p,p1,p2,estar];
ListPlot[{Table[{Ratios4κ[[j,1]],Chop[p1/.p1->Ratios4κ[[j,13]],(10^(-8))]},{j,1,κnum+1}],Table[{Ratios4κ[[j,1]],Chop[p2/.p2->Ratios4κ[[j,14]],(10^(-8))]},{j,1,κnum+1}],Table[{Ratios4κ[[j,1]],Chop[p/.p->Ratios4κ[[j,12]],(10^(-8))]},{j,1,κnum+1}]},AxesStyle -> Thickness[0.0025],PlotTheme->{"Detailed"},LabelStyle->{FontSize->12,Black},PlotLegends-> Placed [{(Subscript[p, 1]),(Subscript[p, 2]),p},{Right,Center}],FrameLabel->{Style["κ"
, FontFamily->"Times",FontSize -> 20],"p [Pa]"}]

Srovnání s Daltonovým zákonem

Clear[κ,j,p,p1,p2,estar];
ListPlot[{Table[{Ratios4κ[[j,1]],Chop[(p1-((ρ1)/(ρ1+(κ^(-1))ρ2)) p)/(((ρ1)/(ρ1+(κ^(-1))ρ2))p)/.κ->Ratios4κ[[j,1]]/.ρ1->Ratios4κ[[j,6]]/.ρ2->Ratios4κ[[j,7]]/.e1->Ratios4κ[[j,8]]/.e2->Ratios4κ[[j,9]]/.s1->Ratios4κ[[j,10]]/.s2->Ratios4κ[[j,11]]/.p->Ratios4κ[[j,12]]/.p1->Ratios4κ[[j,13]]/.p2->Ratios4κ[[j,14]]/.estar->Ratios4κ[[j,15]],(10^(-8))]},{j,1,κnum+1}],Table[{Ratios4κ[[j,1]],Chop[100*(p2-(((κ^(-1))ρ2)/(ρ1+(κ^(-1))ρ2)) p)/((((κ^(-1))ρ2)/(ρ1+(κ^(-1))ρ2))p)/.κ->Ratios4κ[[j,1]]/.ρ1->Ratios4κ[[j,6]]/.ρ2->Ratios4κ[[j,7]]/.e1->Ratios4κ[[j,8]]/.e2->Ratios4κ[[j,9]]/.s1->Ratios4κ[[j,10]]/.s2->Ratios4κ[[j,11]]/.p->Ratios4κ[[j,12]]/.p1->Ratios4κ[[j,13]]/.p2->Ratios4κ[[j,14]]/.estar->Ratios4κ[[j,15]],(10^(-8))]},{j,1,κnum+1}]},AxesStyle -> Thickness[0.0025],PlotTheme->{"Detailed"},LabelStyle->{FontSize->12,Black},PlotLegends-> Placed [{((Subscript[p, 1])-(Subscript[x, 1])p)/((Subscript[x, 1])p),((Subscript[p, 2])-(Subscript[x, 2])p)/((Subscript[x, 2])p)},{Right,Top}],PlotMarkers->{{"▼",11},{"●",6}},FrameLabel->{Style["κ"
, FontFamily->"Times",FontSize -> 20],"[%]"}]

Srovnání s rozkladem tlaku dle hustotního zlomku

Clear[κ,j,p,p1,p2,estar];
ListPlot[{Table[{Ratios4κ[[j,1]],Chop[((((ρ1)/ρ)(e1+e2)-e1)-(1/estar)(((ρ1)/ρ)(s1+s2)-s1))/(((ρ1)/ρ)p)/.ρ1->Ratios4κ[[j,6]]/.ρ2->Ratios4κ[[j,7]]/.e1->Ratios4κ[[j,8]]/.e2->Ratios4κ[[j,9]]/.s1->Ratios4κ[[j,10]]/.s2->Ratios4κ[[j,11]]/.p->Ratios4κ[[j,12]]/.p1->Ratios4κ[[j,13]]/.p2->Ratios4κ[[j,14]]/.estar->Ratios4κ[[j,15]],(10^(-8))]},{j,1,κnum+1}],Table[{Ratios4κ[[j]][[1]],Chop[((((ρ2)/ρ)(e1+e2)-e2)-(1/estar)(((ρ2)/ρ)(s1+s2)-s2))/(((ρ2)/ρ)p)/.ρ1->Ratios4κ[[j,6]]/.ρ2->Ratios4κ[[j,7]]/.e1->Ratios4κ[[j,8]]/.e2->Ratios4κ[[j,9]]/.s1->Ratios4κ[[j,10]]/.s2->Ratios4κ[[j,11]]/.p->Ratios4κ[[j,12]]/.p1->Ratios4κ[[j,13]]/.p2->Ratios4κ[[j,14]]/.estar->Ratios4κ[[j,15]],(10^(-8))]},{j,1,κnum+1}]},AxesStyle -> Thickness[0.0025],PlotTheme->{"Detailed"},LabelStyle->{FontSize->12,Black},PlotLegends-> Placed [{((Subscript[p, 1])-(Subscript[δ, 1]) p)/((Subscript[δ, 1])p),((Subscript[p, 2])-(Subscript[δ, 2])p)/((Subscript[δ, 2])p)},{Right,Center}],FrameLabel->{Style["κ"
, FontFamily->"Times",FontSize -> 20]," "}]

Směs dvou van der Waalsových pynů

Quit  [];

s1=(3/2)(kB/m1) ρ1 Log[(T1/Tref)]-(kB/m1) ρ1 Log[C1 ((ρ1/m1)/(1-ν1 ρ1/m1))]/.T1->(e1-(1/2)((u1^2)/(ρ1))+a1 ((ρ1)^2))((((3/2)(kB/m1)ρ1))^(-1));
s2=(3/2)(kB/m2) ρ2 Log[(T2/Tref)]-(kB/m2) ρ2 Log[C2 ((ρ2/m2)/(1-ν2 ρ2/m2))]/.T2->(e2-(1/2)((u2^2)/(ρ2))+a2 ((ρ2)^2))((((3/2)(kB/m2)ρ2))^(-1));
s=s1+s2;

FirstMEstep={D[s,ρ1]==ρstar,D[s,ρ2]==ρstar,D[s,u1]==ustar,D[s,u2]==ustar,D[s,e1]==estar,D[s,e2]==estar}//Simplify;

Projections={ρ==ρ1+ρ2,u==u1+u2,e==e1+e2};

kB=1.380649 (10^(-23));(*J/K*)
h=6.62607015 (10^(-34)); (*J/Hz*)
Da=1.66053906660 (10^(-27));(*kg*)

(*Celkové stavové proměnné*)
ρ=1; (*kg/(m^3)*)
u=0 (10^3); (*kg m/s*) 
T=300;  (*K*)

(*Ostatní veličiny*)
m2=10 Da;
Tref=300; (*K*)
C2=(1/10)(m2/ρ);
ν1=(10^(-28));
ν2=0.5 (10^(-28));
a2=3 (10^(-49))/(m2^2);
estar=1/T;

(*Parametry výpočtu*)
κnum=100;κmin=0.1;κmax=5;Ratios4κ=Table[{},{i,1,κnum+1}];

Clear[κ,j,m1,e,p,p1,p2];
For[j=1,j≤κnum+1,j++,
κ=κmin+(j-1) (κmax-κmin)/κnum;
m1=m2/κ; (*κ=m2/m1*) 
C1=(1/10)(m1/ρ);
a1=3 (10^(-49))/(m1^2);
eestimate=(3/2)kB/estar((ρ/((m1+m2)/2)))+(1/2)((u^2)/ρ);
numEstimate=Chop[FindRoot[FirstMEstep⋃Projections,{{ρ1,ρ(0.8-(0.6/4)κ)},{ρ2,ρ(1-(0.8-(0.6/4)κ))},{u1,u/2},{u2,u/2},{e1,eestimate/2},{e2,eestimate/2},{ρstar,(10^3)},{ustar,1},{e,eestimate}}]];
p=(-e+ρ μ+T s)/.μ->(-ρstar T)/.numEstimate;
p1=-e1+ρ1 μ+T s1/.μ->(-ρstar T)/.numEstimate;
p2=-e2+ρ2 μ+T s2/.μ->(-ρstar T)/.numEstimate;
Ratios4κ[[j]]={κ,ρ1/.numEstimate,ρ2/.numEstimate,u1/.numEstimate,u2/.numEstimate,e1/.numEstimate,e2/.numEstimate,p,p1,p2,ρstar/.numEstimate}];

Parciální hustoty

Clear[κ,j,m1,p,p1,p2];
ListPlot[{Table[{Ratios4κ[[j,1]],ρ1/.ρ1->Ratios4κ[[j,2]]},{j,1,κnum+1}],Table[{Ratios4κ[[j,1]],ρ2/.ρ2->Ratios4κ[[j,3]]},{j,1,κnum+1}],Table[{Ratios4κ[[j,1]],ρ1+ρ2/.ρ1->Ratios4κ[[j,2]]/.ρ2->Ratios4κ[[j,3]]},{j,1,κnum+1}]},AxesStyle -> Thickness[0.0025],PlotTheme->{"Detailed"},LabelStyle->{FontSize->12,Black},PlotLegends-> Placed [{"(Subscript[ρ, 1])","(Subscript[ρ, 2])","ρ"  },{Right,Center}],FrameLabel->{Style["κ"
, FontFamily->"Times",FontSize -> 20],"ρ [kg/(m^3)]"}]

Celková hustota celkové energie a parciální hustoty celkových energiií

Clear[κ,j,m1,p,p1,p2];
ListPlot[{Table[{Ratios4κ[[j,1]],e1/.e1->Ratios4κ[[j,6]]},{j,1,κnum+1}],Table[{Ratios4κ[[j,1]],e2/.e2->Ratios4κ[[j,7]]},{j,1,κnum+1}],Table[{Ratios4κ[[j,1]],e1+e2/.e1->Ratios4κ[[j,6]]/.e2->Ratios4κ[[j,7]]},{j,1,κnum+1}]},AxesStyle -> Thickness[0.0025],PlotTheme->{"Detailed"},LabelStyle->{FontSize->12,Black},PlotLegends-> Placed [{"(Subscript[e, 1])","(Subscript[e, 2])",e  },{Right,Center}],FrameLabel->{Style["κ"
, FontFamily->"Times",FontSize -> 20],"e [J/(m^3)]"}]

Celkový tlak a parciální tlaky

Clear[κ,j,m1,p,p1,p2];
ListPlot[{Table[{Ratios4κ[[j,1]],p1/.p1->Ratios4κ[[j,9]]},{j,1,κnum+1}],Table[{Ratios4κ[[j,1]],p2/.p2->Ratios4κ[[j,10]]},{j,1,κnum+1}],Table[{Ratios4κ[[j,1]],p/.p->Ratios4κ[[j,8]]},{j,1,κnum+1}]},AxesStyle -> Thickness[0.0025],PlotTheme->{"Detailed"},LabelStyle->{FontSize->12,Black},PlotLegends-> Placed [{(Subscript[p, 1]),(Subscript[p, 2]),p},{Right,Center}],FrameLabel->{Style["κ"
, FontFamily->"Times",FontSize -> 20],"p [Pa]"}]

Srovnání s Daltonovým zákonem

Clear[κ,j,m1,p,p1,p2];
ListPlot[{Table[{Ratios4κ[[j,1]], 100(p1-((ρ1/m1)/(ρ1/m1+ρ2/m2))p)/(((ρ1/m1)/(ρ1/m1+ρ2/m2)) p)/.m1->m2/Ratios4κ[[j,1]]/.ρ1->Ratios4κ[[j,2]]/.ρ2->Ratios4κ[[j,3]]/.p->Ratios4κ[[j,8]]/.p1->Ratios4κ[[j,9]]},{j,1,κnum+1}],Table[{Ratios4κ[[j,1]],100(p2-((ρ2/m2)/(ρ1/m1+ρ2/m2))p)/(((ρ2/m2)/(ρ1/m1+ρ2/m2)) p)/.m1->m2/Ratios4κ[[j,1]]/.ρ1->Ratios4κ[[j,2]]/.ρ2->Ratios4κ[[j,3]]/.p->Ratios4κ[[j,8]]/.p2->Ratios4κ[[j,10]]},{j,1,κnum+1}]},AxesStyle -> Thickness[0.0025],PlotTheme->{"Detailed"},LabelStyle->{FontSize->12,Black},PlotLegends-> Placed [{((Subscript[p, 1])-(Subscript[x, 1])p)/((Subscript[x, 1])p),((Subscript[p, 2])-(Subscript[x, 2])p)/((Subscript[x, 2])p)},{Right,Center}],FrameLabel->{Style["κ"
, FontFamily->"Times",FontSize -> 20],"[%]"}]

Srovnání s rozkladem tlaků dle hustotního zlomku

Clear[κ,j,m1,p,p1,p2];
ListPlot[{Table[{Ratios4κ[[j,1]],((((ρ1)/ρ)(e1+e2)-e1)-T(((ρ1)/ρ)s-s1))/(((ρ1)/ρ) p)/.m1->m2/Ratios4κ[[j,1]]/.ρ1->Ratios4κ[[j,2]]/.ρ2->Ratios4κ[[j,3]]/.u1->Ratios4κ[[j,4]]/.u2->Ratios4κ[[j,5]]/.e1->Ratios4κ[[j,6]]/.e2->Ratios4κ[[j,7]]/.p->Ratios4κ[[j,8]] },{j,1,κnum+1}],Table[{Ratios4κ[[j,1]],((((ρ2)/ρ)(e1+e2)-e2)-T(((ρ2)/ρ)s-s2))/(((ρ2)/ρ) p)/.m1->m2/Ratios4κ[[j,1]]/.ρ1->Ratios4κ[[j,2]]/.ρ2->Ratios4κ[[j,3]]/.u1->Ratios4κ[[j,4]]/.u2->Ratios4κ[[j,5]]/.e1->Ratios4κ[[j,6]]/.e2->Ratios4κ[[j,7]]/.p->Ratios4κ[[j,8]] },{j,1,κnum+1}]},AxesStyle -> Thickness[0.0025],PlotTheme->{"Detailed"},LabelStyle->{FontSize->12,Black},PlotLegends-> Placed [{((Subscript[p, 1])-(Subscript[δ, 1]) p)/((Subscript[δ, 1])p),((Subscript[p, 2])-(Subscript[δ, 2])p)/((Subscript[δ, 2])p)},{Right,Center}],FrameLabel->{Style["κ"
, FontFamily->"Times",FontSize -> 20]," "}]