I'm writing code to simulate a chemical reaction. This requires the calculation of the fugacity of a mixture of gases. The underlying model is the SRK equation of state.
First the physical properties are declared:
*fugacity calculations according to SRK equation*)
(*Critical properties for selected gases*)
r = 8.314;
presc = {35.0, 73.8, 20.5, 220.5, 81.0, 33.9}*100000;
tempc = {132.9, 304.2, 43.6, 647.3, 512.6, 126.2};
volc = {93.1, 94.0, 51.5, 56.0, 118.0}*10^-6;
omegac = {0.049, 0.255, 0.0, 0.344, 0.572, 0.040};
The equation of state is set up:
(*SRK equation*)
a = 0.42748*(r*tempc)^2/presc;
b = 0.08664*(r*tempc)/presc;
alpha[temp_] = (1 + (0.480 + 1.574 omegac - 0.176 omegac^2) (1-
Sqrt[temp/tempc]))^2;
A[temp_, pres_] = Table[(a[[i]]*alpha[temp][[i]]*pres)/(r*temp)^2, {i, 1,
6}];
B[temp_, pres_] = Table[(b[[i]]*pres)/(r*temp), {i, 1, 6}];
Amix[temp_, pres_, yfrac_] := Sum[yfrac[[i]]*yfrac[[j]]*Sqrt[A[temp,
pres][[i]]*A[temp, pres][[j]]], {i, 1, 6}, {j, 1, 6}];
Bmix[temp_, pres_, yfrac_] := Sum[yfrac[[i]]*B[temp, pres][[i]], {i, 1, 6}];
Finally, the fugacity is calculated:
(Fugacity calculation)
fugacitymix[temp_, pres_, i_, yfrac_] := With[{Amix2 = Amix[temp, pres,
yfrac], Bmix2 = Bmix[temp, pres, yfrac]},
Exp[(z - 1) B[temp, pres][[i]]/Bmix2 - Log[z -
Bmix2] - Amix2/Bmix2 (2 (A[temp, pres][[i]]/Amix2)^0.5 -
B[temp, pres][[i]]/Bmix2) Log[1 + Bmix2/z] /.
z -> Last[
Cases[z /. {ToRules@NRoots[z^3 -
z^2 + (Amix2 - Bmix2 - (Bmix2)^2) z - (Amix2) (Bmix2) ==
0, z]}, _Real]]
]]
The issue is that the code is somewhat slow. It makes running the model costly. For example:
yfra = N@{3, 20, 75, 1, 1, 1}/100;
Plot[Table[fugacitymix[273 + 200, x*100000, i, yfra], {i, 1, 6}], {x,
1, 300}] // AbsoluteTiming
Takes approximately 40s on my machine.
I tried compiling the function with:
cfugacity =
Compile[{{temp, _Real}, {pres, _Real}, {i, _Integer}, {yfrac, _Real,
6}},
With[{Amix2 = Amix[temp, pres, yfrac],
Bmix2 = Bmix[temp, pres, yfrac]},
Exp[
(z - 1) B[temp, pres][[i]]/Bmix2 - Log[z - Bmix2] -
Amix2/Bmix2 (2 (A[temp, pres][[i]]/Amix2)^0.5 -
B[temp, pres][[i]]/Bmix2) Log[1 + Bmix2/z] /.
z -> Last[
Cases[z /. {ToRules@
NRoots[z^3 -
z^2 + (Amix2 - Bmix2 - (Bmix2)^2) z - (Amix2) (Bmix2) ==
0, z]}, _Real]]
]]
]
Plot[Table[cfugacity[273 + 200, x*100000, i, yfra], {i, 1, 6}], {x,
1, 300}] // AbsoluteTiming
This one takes in fact longer than the uncompiled version, about 47s. I'm looking for suggestions to improve the speed.