192 std::vector<double> M_L, std::vector<double> Tc_L, std::vector<double> pc_L, std::vector<double> omega_, std::vector<double> devi,
193 std::vector<double> alpha_, std::vector<double> gI_L, std::vector<double> gII_L) :
194p(p_), T(T_), z(z_), num_comp(z.size()), a_c(a_L), b(b_L), M(M_L), T_c(Tc_L), p_c(pc_L), omega(omega_), m(devi), alpha(alpha_), g_I(gI_L), g_II(gII_L)
239 for(
int j = 0; j < num_comp; j++){
240 sum1 = sum1 + x_i[j]*g_jk[num_comp*j+i]*beta_jk[num_comp*j+i];
241 sum2 = sum2 + x_i[j]*beta_jk[num_comp*j+i];
245 for(
int k = 0; k < num_comp; k++){
246 sub31 = sub31 + x_i[k]*beta_jk[j+num_comp*k];
247 sub41 = sub41 + x_i[k]*g_jk[j+num_comp*k]*beta_jk[j+num_comp*k];
248 sub42 = sub42 + x_i[k]*beta_jk[j+num_comp*k];
250 sum3 = sum3 + x_i[j]*g_jk[num_comp*i+j]*beta_jk[num_comp*i+j] / sub31;
251 sum4 = sum4 + x_i[j]*beta_jk[num_comp*i+j]*sub41/
util::pow(sub42,2);
253 return 1/(R*_T) * (a_i / b[i] - 1/lambda * (sum1/sum2 + sum3 - sum4));
276 std::vector<double> a_i(num_comp), g_jk(num_comp*num_comp), beta_jk(num_comp*num_comp);
277 for(
int i = 0; i < num_comp; i++){
279 for(
int j = 0; j < num_comp; j++){
280 g_jk[num_comp*i+j] = g_I[num_comp*i+j] + g_II[num_comp*i+j] * T;
281 beta_jk[num_comp*i+j] = b[j] *
util::exp(-alpha[num_comp*i+j]*g_jk[num_comp*i+j]/(R*T));
284 std::vector<double> vx(2+2*num_comp), K(num_comp), x(num_comp), y(num_comp), G_E(2), b_m(2), a_m(2), gammas(2*num_comp);
285 std::vector<double> deltaLogFugacity(num_comp), errorFugacity(num_comp);
286 double totalError = 0;
287 for(
int i = 0; i < num_comp; i++){
289 K[i] = 1.0*(p_c[i]/(p))*
util::exp(5.373*(1+omega[i])*(1-T_c[i]/T));
290 errorFugacity[i] = 1;
292 for(
int i = 0; i < 2; i++){
300 for(
int i = 0; i < num_comp; i++){
301 f += (K[i]-1)*z[i]/(1+(K[i]-1)*beta);
303 while(
abs(f) > residuum){
305 for(
int i = 0; i < num_comp; i++){
306 dfdbeta += -1*(K[i]-1)*(K[i]-1)*z[i]/((1+(K[i]-1)*beta)*(1+(K[i]-1)*beta));
308 beta = beta - f / dfdbeta;
310 for(
int i = 0; i < num_comp; i++){
311 f += (K[i]-1)*z[i]/(1+(K[i]-1)*beta);
316 for(
int i = 0; i < num_comp; i++){
317 x[i] = z[i]/(1+(K[i]-1)*beta);
322 if(x_sum != 1 || y_sum != 1){
323 for(
int i = 0; i < num_comp; i++){
324 vx[i+2] = x[i]/x_sum;
326 vx[i+2+num_comp] = y[i]/y_sum;
327 y[i] = vx[i+2+num_comp];
331 G_E[0] =
G_Excess(x, g_jk, beta_jk);
332 G_E[1] =
G_Excess(y, g_jk, beta_jk);
335 a_m[0] =
a_mix(x, a_i, b_m[0], G_E[0]);
336 a_m[1] =
a_mix(y, a_i, b_m[1], G_E[1]);
337 for(
int i = 0; i < num_comp; i++){
338 gammas[i] =
gamma_i(x, a_i[i], g_jk, beta_jk, T, i);
339 gammas[i+num_comp] =
gamma_i(y, a_i[i], g_jk, beta_jk, T, i);
342 for(
int i = 0; i < 2; i++){
343 while(
abs(p-
p_PR(vx[i], a_m[i], b_m[i])) > residuum){
344 double dpdv = -R*T/
util::pow(vx[i]-b_m[i],2) +
346 vx[i] = vx[i] - (p -
p_PR(vx[i], a_m[i], b_m[i])) / (-dpdv);
349 for(
int i = 0; i < num_comp; i++){
350 deltaLogFugacity[i] =
lnfugacity_i(x, vx[0], a_m[0], b_m[0], gammas[i], T, i) -
351 lnfugacity_i(y, vx[1], a_m[1], b_m[1], gammas[i+num_comp], T, i);
352 errorFugacity[i] =
util::exp(deltaLogFugacity[i]) - 1;
353 K[i] =
util::exp(deltaLogFugacity[i]) / (x[i]/y[i]);
354 totalError +=
util::pow(errorFugacity[i],2);
356 }
while(
util::pow(totalError,0.5) > residuum);
357 for(
int i = 1; i<num_comp; i++){
358 volatilities[i] = (M[i]*(vx[i+2]/vx[0] - vx[i+2+num_comp]/vx[1])) / (M[0]*(vx[2]/vx[0] - vx[2+num_comp]/vx[1]));
361 for(
int i = 0; i<num_comp; i++){
362 vol_sum += volatilities[i];
364 for(
int i = 0; i<num_comp; i++){
365 chis[i] = volatilities[i] / vol_sum;
MultiComponentPengRobinson(double p_, double T_, std::vector< double > z_, std::vector< double > a_L, std::vector< double > b_L, std::vector< double > M_L, std::vector< double > Tc_L, std::vector< double > pc_L, std::vector< double > omega_, std::vector< double > devi, std::vector< double > alpha_, std::vector< double > gI_L, std::vector< double > gII_L)