Matlab uppgift: van der Waals ekvation.

The van der Waals equation of state, (p+a/v²)(v-b)=RT
relates the pressure p, the specific volume v, and the temperature T of a gas,
where R is a universal constant and a and b are constants that depend on the
particular gas. In appropriate units, R = 0,082054, and for carbon dioxide,
a = 3,592 and b = 0,04267. Use a zero finder to compute the specific
volume v for a temperature of T = 300 and for pressures p of 1, 10, and
100. Compare your results with those of the ideal gas law, pv = RT. The
latter can be used as a starting guess for an iterative method to solve the van
der Waals equation (if one is needed).

Så jag har satt upp två m-filer med lite hjälp av en kamrat:

m-fil 1:
function f = we(v)

a = 3.592;
b = 0.04267;
r = 0.082054;
T = 300;
p = 100;

f = (p+(a/v^2))*(v-b)-r*T;

m-fil 2:
fsolve (@we,24)

r = 0.082054;
T = 300;
p = 100;
v = r*T/p
error = abs(v-ans)

Hur ser det ut i era ögon?

1)
Jag får en "error message" på första m-filen när jag försöker köra den: variable "v" undefined. Det kankse bara är mening att den andra m-filen ska fungera genom den första?

2)
Skulle behöva förklaring till utrycket fsolve (@we,24). Som jag förstår det är "fsolve" en lösare till ickelinjära ekvationer. Men var är ideen bakom (@we,24). Och varför 24?

mvh

Först ett alternativ som inte kräver några m-filer.

>> R= 0.082054;
>> p = [1 10 100];
>> T = 300;
>> v_ideal = R*T./p;
>> a = 3.592;
>> b = 0.04267;
>> v_vw(1)=fzero(@(v) (p(1)+a./v^2).*(v-b)-R*T,v_ideal(1));
>> v_vw(2)=fzero(@(v) (p(2)+a./v^2).*(v-b)-R*T,v_ideal(2));
>> v_vw(3)=fzero(@(v) (p(3)+a./v^2).*(v-b)-R*T,v_ideal(3));

Ideala gaslagens volymer finns då i vektorn v_ideal och van der Waals volymer i v_vw.

1) Ska du köra den första måste du ange v genom att kalla den som ex. we(1).
2) @we anger vilken ekvation som sätts lika med 0 och löses. Man kan också skriva 'we'. Siffran 24 är en startgissning för v som är passande för p = 1 men inte för p = 100 vilket är angivet i filerna.

Jag vet inte om "roots" inte får användas, annars kan ju ekvationen:

(p+a/v²)(v-b)=RT

Förvandlas till en tredjegradare i v och med roots få lösningarna (där två kommer bli komplexa och en bli reell). På så vis slipper man definera funktionen och kan jämföra med allmäna gaslagen lätt.

Sitter på samma uppgift: men med följande:
function f = we(v)
v = w(1)
a = 3.592;
b = 0.04267;
R = 0.082054;
T = 300;
p = 100;
f = (p+(a/v^2))*(v-b)-R*T;

får man en error message.

Ok. Hur skulle det se ut i matlab i så fall?

Observera att koden jag skrivit tar att man kan ange a,b,T också. Om de är fixa (som bara för ett problem) så kan man modifiera koden efter tycke och smak. Man kan också modifiera så att det går att ta par av (P,T) i vektorer om man känner för det, dvs givna tryck och temperaturer medan repulsion- och attraktionskonstanterna är fixa.


Nu ligger de från VDW-ekvation i Vw och de från ideal i VI. Observera att P är en tryckvektor nu, men kan även vara ett rent tal. Exempel:

vanw(3.592,0.04267,[1 10 100],300) ger:

VI =

24.6162 2.4616 0.2462


ans =

24.5126 2.3545 0.0795

Och ans är här alltså VW. Man kan naturligtvis trycka bort så inte VI syns eller modifiera på annat vis.

För härledning av tredjegradspolynomet så se:
(p + a/v^2)(v - b) = RT

Multiplicera med v^2 ger:
(p*v^2 + a)(v - b) = RT*v^2
p*v^3 - b*p*v^2 + a*v - a*b = RT*v^2
v^3*(p) + v^2*(-b*p-RT) + a*v + (-a*b) = 0

roots tar ett polynoms koefficienter och löser ekvationen. Med raden Vw(i) = sol(sol == real(sol)) ; så tar vi bort de komplexa lösningarna. Observera att om VDW ekvation ger flera reella lösningar (jag vet inte om det scenariot kan uppkomma) blir det problem.

Ok, ska gå igenom detta, tack.

Det första alternativet förstår jag väldigt lite utav. Men om jag håller mig till m-filerna som jag förstår behöver jag bara definiera "v" i m-fil 1 så är det klart. om jag definierar v = we(1) får jag följande feedback:
Maximum recursion limit of 500 reached...

Varför ska du göra det? Den andra filen kallar den första. Sen ska du välja p = 1 om du har startgissning 24 som sagt. Försök förstå vad det är du gör.

Ok men det var det jag trodde från början. Jag nämnde att m-fil 1 inte fungerar när man kör den. men att det kanske berodde på att den inte ska fungera självständigt utan bara genom m-fil 2.
Stämmer detta?

Titta på vad första filen gör. Vad har du för värde av den separat?

Första m-filen ställer up van der Waals ekvationen så att den = 0 och sedan ska den väl lösa ut v för att kunna skicka den till m-fil 2?