【发布时间】：2020-08-24 22:19:34
【问题描述】：
我必须找到以下等式的零：
这是一个状态方程，如果您不确切知道 EoS 是什么，这并不重要。根据上述等式的根，我计算（除其他外）气态物质在不同压力和温度下的压缩因子 Z。通过这些解决方案，我可以绘制以压力为横坐标、Zs 为纵坐标、温度为参数的曲线族。 Beta、delta、eta 和 phi 是常数，pr 和 Tr 也是如此。
在用 Newton-Raphson 方法（与其他几个 EoS 配合得很好）失败后，我决定尝试 Scipy 的 root() 函数。令我不满的是，我得到了这张图表：
我们很容易看出，这张锯齿状的图表是完全有缺陷的。我应该得到平滑的曲线。此外，Z 通常介于 0.25 和 2.0 之间。因此，等于或大于 3 的 Z 完全不符合标准。然而 Z
然后我尝试了 Octave 的 fzero() 求解器，得到了这个：
这正是我应该得到的，因为它们是具有正确/预期形状的曲线！
我的问题来了。显然 Scipy 的 root() 和 Octave 的 fzero() 基于来自 MINPACK 的相同算法hybrid。尽管如此，结果显然并不相同。有谁知道为什么吗？
我绘制了一条由 Octave（横坐标）获得的 Zs 与使用 Scipy 获得的 Zs 的曲线，并得到了这个：
底部暗示直线的点代表y = x，即 Octave 和 Scipy 在他们提出的解决方案中同意的点。其他观点完全不同，不幸的是，它们太多了，不能简单地忽略。
从现在开始我可能会一直使用 Octave，因为它可以工作，但我想继续使用 Python。
您对此有何看法？有什么建议吗？
PS：这是原始 Python 代码。它会生成此处显示的第一个图表。
import numpy
from scipy.optimize import root
import matplotlib.pyplot as plt
def fx(x, beta, delta, eta, phi, pr_, Tr_):
tmp = phi*x**2
etmp = numpy.exp(-tmp)
f = x*(1.0 + beta*x + delta*x**4 + eta*x**2*(1.0 + tmp)*etmp) - pr_/Tr_
return f
def zsbwr(pr_, Tr_, pc_, Tc_, zc_, w_, MW_, phase=0):
d1 = 0.4912 + 0.6478*w_
d2 = 0.3000 + 0.3619*w_
e1 = 0.0841 + 0.1318*w_ + 0.0018*w_**2
e2 = 0.075 + 0.2408*w_ - 0.014*w_**2
e3 = -0.0065 + 0.1798*w_ - 0.0078*w_**2
f = 0.770
ee = (2.0 - 5.0*zc_)*numpy.exp(f)/(1.0 + f + 3.0*f**2 - 2*f**3)
d = (1.0 - 2.0*zc_ - ee*(1.0 + f - 2.0*f**2)*numpy.exp(-f))/3.0
b = zc_ - 1.0 - d - ee*(1.0 + f)*numpy.exp(-f)
bc = b*zc_
dc = d*zc_**4
ec = ee*zc_**2
phi = f*zc_**2
beta = bc + 0.422*(1.0 - 1.0/Tr_**1.6) + 0.234*w_*(1.0- 1.0/Tr_**3)
delta = dc*(1.0+ d1*(1.0/Tr_ - 1.0) + d2*(1.0/Tr_ - 1.0)**2)
eta = ec + e1*(1.0/Tr_ - 1.0) + e2*(1.0/Tr_ - 1.0)**2 \
+ e3*(1.0/Tr_ - 1.0)**3
if Tr_ > 1:
y0 = pr_/Tr_/(1.0 + beta*pr_/Tr_)
else:
if phase == 0:
y0 = pr_/Tr_/(1.0 + beta*pr_/Tr_)
else:
y0 = 1.0/zc_**(1.0 + (1.0 - Tr_)**(2.0/7.0))
raiz = root(fx,y0,args=(beta, delta, eta, phi, pr_, Tr_),method='hybr',tol=1.0e-06)
return pr_/raiz.x[0]/Tr_
if __name__ == "__main__":
Tc = 304.13
pc = 73.773
omega = 0.22394
zc = 0.2746
MW = 44.01
Tr = numpy.array([0.8, 0.93793103])
pr = numpy.linspace(0.5, 14.5, 25)
zfactor = numpy.zeros((2, 25))
for redT in Tr:
j = numpy.where(Tr == redT)[0][0]
for redp in pr:
indp = numpy.where(pr == redp)[0][0]
zfactor[j][indp] = zsbwr(redp, redT, pc, Tc, zc, omega, MW, 0)
for key, value in enumerate(zfactor):
plt.plot(pr, value, '.-', linewidth=1, color='#ef082a')
plt.figure(1, figsize=(7, 6))
plt.xlabel('$p_R$', fontsize=16)
plt.ylabel('$Z$', fontsize=16)
plt.grid(color='#aaaaaa', linestyle='--', linewidth=1)
plt.show()
现在是 Octave 脚本：
function SoaveBenedictWebbRubin
format long;
nTr = 11;
npr = 43;
ic = 1;
nome = {"CO2"; "N2"; "H2O"; "CH4"; "C2H6"; "C3H8"};
comp = [304.13, 73.773, 0.22394, 0.2746, 44.0100; ...
126.19, 33.958, 0.03700, 0.2894, 28.0134; ...
647.14, 220.640, 0.34430, 0.2294, 18.0153; ...
190.56, 45.992, 0.01100, 0.2863, 16.0430; ...
305.33, 48.718, 0.09930, 0.2776, 30.0700; ...
369.83, 42.477, 0.15240, 0.2769, 44.0970];
Tc = comp(ic,1);
pc = comp(ic,2);
w = comp(ic,3);
zc = comp(ic,4);
MW = comp(ic,5);
Tr = linspace(0.8, 2.8, nTr);
pr = linspace(0.2, 7.2, npr);
figure(1, 'position',[300,150,600,500])
for i=1:size(Tr, 2)
icont = 1;
zval = zeros(1, npr);
for j=1:size(pr, 2)
[Z, phi, density] = SBWR(Tr(i), pr(j), Tc, pc, zc, w, MW, 0);
zval(icont) = Z;
icont = icont + 1;
endfor
plot(pr,zval,'o','markerfacecolor','white','linestyle','-','markersize',3);
hold on;
endfor
str = strcat("Soave-Benedict-Webb-Rubin para","\t",nome(ic));
xlabel("p_r",'fontsize',15);
ylabel("Z",'fontsize',15);
title(str,'fontsize',12);
end
function [Z,phi,density] = SBWR(Tr, pr, Tc, pc, Zc, w, MW, phase)
R = 8.3144E-5; % universal gas constant (bar·m3/(mol·K))
% Definition of parameters
d1 = 0.4912 + 0.6478*w;
d2 = 0.3 + 0.3619*w;
e1 = 0.0841 + 0.1318*w + 0.0018*w**2;
e2 = 0.075 + 0.2408*w - 0.014*w**2;
e3 = -0.0065 + 0.1798*w - 0.0078*w**2;
f = 0.77;
ee = (2.0 - 5.0*Zc)*exp(f)/(1.0 + f + 3.0*f**2 - 2.0*f**3);
d = (1.0 - 2.0*Zc - ee*(1.0 + f - 2.0*f**2)*exp(-f))/3.0;
b = Zc - 1.0 - d - ee*(1.0 + f)*exp(-f);
bc = b*Zc;
dc = d*Zc**4;
ec = ee*Zc**2;
ff = f*Zc**2;
beta = bc + 0.422*(1.0 - 1.0/Tr**1.6) + 0.234*w*(1.0 - 1.0/Tr**3);
delta = dc*(1.0 + d1*(1.0/Tr - 1.0) + d2*(1.0/Tr - 1.0)**2);
eta = ec + e1*(1.0/Tr - 1.0) + e2*(1.0/Tr - 1.0)**2 + e3*(1.0/Tr - 1.0)**3;
if Tr > 1
y0 = pr/Tr/(1.0 + beta*pr/Tr);
else
if phase == 0
y0 = pr/Tr/(1.0 + beta*pr/Tr);
else
y0 = 1.0/Zc**(1.0 + (1.0 - Tr)**(2.0/7.0));
end
end
fun = @(y)y*(1.0 + beta*y + delta*y**4 + eta*y**2*(1.0 + ff*y**2)*exp(-ff*y**2)) - pr/Tr;
options = optimset('TolX',1.0e-06);
yi = fzero(fun,y0,options);
Z = pr/yi/Tr;
density = yi*pc*MW/(1000.0*R*Tc);
phi = exp(Z - 1.0 - log(Z) + beta*yi + 0.25*delta*yi**4 - eta/ff*(exp(-ff*yi**2)*(1.0 + 0.5*ff*yi**2) - 1.0));
end
【问题讨论】：
-
向我们展示python代码。可能有几个根要收敛
-
@ev-br 完成。希望它对你有用。否则让我知道。谢谢。
-
@ev-br 当然有多个根。令我惊讶的是，Octave 从提供的初始猜测中找到了正确的。该方法（很可能）与
scipy.optimize.root()使用的方法相同。 -
@ev-br 另外，
zsbwr()函数返回的phi不是上面 f(x) 等式中显示的希腊字母 phi。后者应该是ff。我让我的代码有点混乱。很抱歉。 -
你介意分享八度代码吗？