python - Thermodynamics - VLE code for a binary mixture

Question

I'm a beginner in Python programming, so I've been struggling a lot with it. I've started with a very simple thermodynamics problem but I'm stuck on it. I am going to write my code down here. If someone is available to help, I'd really appreciate it.

So, I need to calculate, system pressure, liquid and vapor compositions using Raoult's Law and compare them with some experimental data. I've done that already (certainly my code is very poor and repetitive, but consider that I'm new to this).

Then, I'd like to store the liquid (x) and vapor (y) compositions in this vectors:

x = np.array([x1, 1 - x1]) and y[i,j] = x[j] * Psat[j] / P[i]) 

I'd like to use the calculated Raoult's Law results as initial guesses so that I can apply them on a cubic equation of state (the Peng-Robinson one), but I don't know how to do that. Could someone please help?

import numpy as np
import math
import matplotlib.pyplot as plt

# EXPERIMENTAL DATA
x1_exp = np.array([0, 0.022, 0.043, 0.083, 0.166, 0.206, 0.274, 0.304, 0.369, 0.452, 0.477, 0.519, 0.584, 0.694, 0.803, 0.851, 1])
y1_exp = np.array([0, 0.139, 0.240, 0.380, 0.532, 0.590, 0.653, 0.680, 0.711, 0.755, 0.764, 0.776, 0.800, 0.840, 0.885, 0.910, 1])
P_exp = np.array([836.30, 997.70, 1128.70, 1392.20, 1922.00, 2166.60, 2573.90, 2726.70, 3004.90, 3449.30, 3562.10, 3756.60, 4035.70, 4492.80, 4966.60, 5159.80, 5723.70])  # kPa
y2_exp = np.array([1 - y1_exp])

# DATA
componentes = 2
T = 293.15
R = 8.314
Tc = np.array([304.15, 369.82])
Pc = np.array([7374.99, 4255.76])
w = np.array([0.225, 0.152])

# SATURATION PRESSURE
# WAGNER'S EQUATION
# 1983 - MCGARRY - Correlation and Prediction of the Vapor Pressure of Pure Liquids
# ln(Pr) = 1/Tr*(a*tau + b*tau^1.5 + c*tau^3 + d*tau^6)
# tau = 1 - Tr

wag_a = np.array([-6.95626, -6.67833])  # SPECIES 1, SPECIES 2
wag_b = np.array([1.19695, 1.15437])
wag_c = np.array([-3.12614, -1.64984])
wag_d = np.array([2.99448, -2.70017])

Tr = np.zeros(componentes)
tau = np.zeros(componentes)
Psat = np.zeros(componentes)

j = 0
for j in range(componentes):
    Tr[j] = T / Tc[j]
    tau[j] = 1 - Tr[j]
    Psat[j] = np.exp(1/Tr[j] * (wag_a[j]*tau[j] + wag_b[j]*tau[j]**1.5 + wag_c[j]*tau[j]**3 + wag_d[j]*tau[j]**6)) * Pc[j]
j = j + 1
print('PRESS?O DE SATURA??O - WAGNER')
print('Psat1 e Psat2 (kPa)')
print(np.around(Psat, 4))

# EQUILIBRIUM - RAOULT'S LAW (TO BE USED AS INITIAL GUESS FOR PENG-ROBINSON)
# RAOULT'S LAW
# y * phi * P = x * gamma * Psat
# P = sum(xi * Pisat)

x1_teorico = np.arange(0, 1.01, 0.05)
x = np.zeros((x1_teorico.size, componentes))
y = np.zeros((x1_teorico.size, componentes))
P = np.zeros(x1_teorico.size)

i = 0
for x1 in x1_teorico:
    x = np.array([x1, 1 - x1])
    P[i] = np.sum(np.dot(x, Psat))
    j = 0
    for j in range(componentes):
        y[i,j] = x[j] * Psat[j] / P[i]
    i = i + 1
    j = j + 1

print()
print('LEI DE RAOULT', '
y1 e y2', 'PRESS?O (kPa)')
for a, b in zip(y, P):
    np.set_printoptions(formatter={'all': lambda y: '{:0.4f}'.format(y)})
    print(a, np.around(b,4), sep='')

# EQUILIBRIUM - PENG-ROBINSON
ac = np.zeros(componentes)
b = np.zeros(componentes)
kappa = np.zeros(componentes)
alpha = np.zeros(componentes)
a = np.zeros(componentes)
a_mis = np.zeros((x1_teorico.size, componentes))
b_mis = np.zeros((x1_teorico.size, componentes))

j = 0
# PURE SPECIES
for j in range(componentes):
    ac[j] = 0.45724 * R**2 * Tc[j]**2 / Pc[j]
    b[j] = 0.07780 * R * Tc[j] / Pc[j]
    kappa[j] = 0.37464 + 1.54226 * w[j] - 0.269932 * w[j]**2
    alpha[j] = [1 + kappa * (1 - math.sqrt(Tr[j]))]**2  **# HERE I GET THE ERROR: TypeError: unsupported operand type(s) for ** or pow(): 'list' and 'int'**
    a[j] = ac[j] * alpha[j]
# MIXTURE
    a_mis[j] = x[j]**2 * a[j] + (1 - x[j])**2 * a[1] + 2 * x[j] * (1-x[j]) * math.sqrt(a[j] * a[1])
    b_mis[j] = x[j] * b[j] + (1 - x[j]) * b[1]
j = j + 1

**# FROM THIS POINT ON I DON'T KNOW HOW TO INSERT PENG-ROBINSON EQUATION OF STATE.**


# CHART - EXPERIMENTAL DATA
plt.scatter(x1_exp, P_exp, label='Fra??o molar fase líquida, $x_1$', marker='s', s=20, facecolors='none',
            edgecolors='blue')
plt.scatter(y1_exp, P_exp, label='Fra??o molar fase vapor, $y_1$', marker='^', s=20, facecolors='none',
            edgecolors='red')
plt.title('Diagrama Pxy CO2(1) + propano(2)')
plt.xlabel('$x_1$, $y_1$')
plt.ylabel('Press?o (kPa)')
plt.legend(loc='best')
plt.xlim(xmin=0, xmax=1)
plt.ylim(ymin=0)

# CHART - THEORETICAL VALUES
plt.plot(x1_teorico, P, color='black', linewidth=1)
plt.plot(y[:,0], P, color='black', linewidth=1)
plt.show()

Answer 1

等待大神答复