Institutt for fysikk, NTNU TFY4165 og FY1005 Termisk fysikk, våren 2011.

Regneøving 9.

(Veiledning: Fredag 18. mars kl. 12.15 - 14.00 og Mandag 21. mars kl. 17.15 - 19.00)

Oppgavene 1 - 4

\underline{Oppgave 5}

For temperaturer T under kritisk temperatur T_{c} gir van der Waals tilstandsligning,

p(V) = \frac{RT}{V − b} − \frac{a}{V^{2}} ,

koeksistens mellom væske og gass for en bestemt verdi av trykket, p_{co}, det såkalte koeksistenstrykket. Vi går ikke inn på begrunnelsen her, men det kan vises at p_{co} fastlegges ved hjelp av Maxwells regel om like arealer, dvs betingelsen

p_{co} · (V_{g} − V_{l}) = ∫_{V_{l}}^{V_{g}} p(V) dV,

der V_{l} og V_{g} er hhv molart væske- og gassvolum.

Bruk verdiene a = 1.368 bar (m^{3}/kmol)^{2} og b = 0.0367 (m^{3}/kmol) for luft (jf øving 1), samt gasskonstanten R = 8.314 J/mol K. Kritisk temperatur for luft er med disse parametrene og van der Waals tilstandsligning omtrent 133 K.

Oppgave: Bruk Maxwells regel om like arealer til å fastlegge koeksistenstrykket p_{co}, med en nøyaktighet på 1 mbar, ved en gitt temperatur T = 123 K. Du kan ta utgangspunkt i Matlab-programmet i øving 1, der van der Waals isotermen for (blant annet) denne temperaturen ble plottet i andre del av oppgave 4.

En enkel algoritme for å løse dette problemet numerisk kan se omtrent slik ut:

• Basert på isotermen p(V), gjett en verdi for væskevolumet V_{l}.

• Bestem tilhørende trykkverdi p(V_{l}).

• Bestem tilsvarende gassvolum V_{g}, dvs slik at p(V_{g}) = p(V_{l})

• Bestem

∫_{V_{l}}^{V_{g}} p(V) dV

numerisk og sammenlign med arealet p(V_{l}) · (V_{g} − V_{l}).

• Gjett en ny verdi for V_{l} og gjenta prosedyren inntil p_{co} er bestemt med den ønskede nøyaktighet.

En innledende del til programmet er tatt med nedenfor, inkludert noen tips for den videre programmeringen:

%%FY1005/TFY4165, Oving 9, Oppgave 5, innledende del (hentet fra oving 1).
%%
%%Velg utskrift av tallverdier med flere gjeldende siffer enn det som er standard
format long;
%%Temperatur i Kelvin
T=123;
%%V = molart volum (L/mol)
Vmin=0.05;
Vmax=0.50;
NV=500;
%%V = vektor med verdier mellom Vmin og Vmax, i alt NV verdier
V=linspace(Vmin,Vmax,NV);
%%Verdier for a og b for luft: 1.368  0.0367
%%Enheter: [a] = bar*(L/mol)^2 og [b] = L/mol
R=8.134;
a=1.368;
b=0.0367;
%%van der Waals tilstandsligning
p = (R*T./(V-b))/100 - a./(V.*V);
%%Plott av isotermen
fig = plot(V,p);
axis([0 Vmax 0 70]);
%%Tekst paa figuren
title('Isoterm: van der Waals gass','fontsize',18);
xlabel('Molart volum (L/mol)','fontsize',18);
ylabel('Trykk (bar)','fontsize',18);
%%Tegn alle kurver i en og samme figur, saa vi kan sammenlikne.
hold on;
%%HERFRA maa du selv skrive matlab-koden, men du faar noen tips:
%%Elementer som maa eller boer inngaa, inklusive eksempler:
%%
%%A Innlesing av gjetning paa vaeskevolumet V1:
%%
%%Eksempel som ber om innlesing av tallverdi for stoerrelsen x1:
%%x1 = input('Les inn tallverdi for x1:  ')
%%

%%B Bestem trykkverdi P1 som tilsvarer volumet V1:
%%
%%Eksempel som bestemmer y1 = x1 - x1^3:
%%y1 = x1 - x1^3;
%%
%%C Bestem tilsvarende gassvolum V2. Tips: Velg en startverdi V_stor som
%%ut fra kurven p(V) helt sikkert er stoerre enn V2, reduser V_stor
%%med en passende skrittlengde dV, og bruk en while (...) ... end loekke
%%til aa lokalisere V2 slik at P2 = p(V2) er (omtrent) lik P1.
%%
%%Eksempel som lokaliserer x2 slik at y2 = x2 - x2^3 er (omtrent) lik
%%y1 = x1 - x1^3 med verdien av x1 et sted mellom -1.15 og -0.58 (valgt
%%slik at vi faktisk har en loesning y2 = y1):
%%x_stor = 1.3;
%%dx = 0.001; %%skrittlengde
%%y_find = x_stor - x_stor^3;
%%while (y_find < y1)
%%      x_stor = x_stor - dx;
%%      y_find = x_stor - x_stor^3;
%%end;
%%x2 = x_stor;
%%y2 = y_find;
%%
%%D Integrer p(V) numerisk fra V1 til V2, og sammenlign med P1*(V2-V1),
%%dvs. arealet under den horisontale linjen.
%%
%%Eksempel som integrerer funksjonen
%%y = x - x^3 fra x1 til x2 med saakalt adaptiv Simpson-kvadratur, og
%%regner ut forskjellen mellom integralet I og arealet y1*(x2-x1) under
%%den horisontale linjen:
%%y = @(x) (x - x^3);
%%I = quad(y, x1, x2);
%%avvik = I - y1*(x2-x1);
%%
%%E Tegn inn linjen mellom (P1,V1) og (P2,V2). Skriv ut beregnede stoerrelser.
%%
%%Eksempel:
%%plot([x1,x2],[y1,y2]);
%%disp('x1,x2: ')
%%disp(x1,x2)
%%disp('y1,y2: ')
%%disp(y1,y2)
%%disp('Avvik mellom areal under kurve og rektangel: ')
%%disp(avvik)
%%
%%Endre x1, og eventuelt ogsaa x_stor og dx, og kjoer paa nytt,
%%for aa oppnaa et mindre avvik.

Til slutt et par kommentarer, samt tips til videre arbeid:

• Merk at skrittlengden dV setter en begrensning på hvor nøyaktig koeksistenstrykket p_{co} kan bestemmes.

• Merk at fortegnet på avviket mellom arealet under kurven p(V) og rektanglet kan brukes til å avgjøre om gjetningen på V_{l} skal justeres opp eller ned.

• Den iterative prosessen for lokalisering av V_{l}, og dermed p_{co}, kan gjøres manuelt eller automatisk. I algoritmen og eksemplet ovenfor er det skissert en manuell løsning, der programmet simpelthen kjøres på nytt og ny verdi for V 1 leses inn fra ”kommandolinjen”.