PR状态方程程序C++

// N2-CH4-CO2.cpp : 定义控制台应用程序的入口点。
//

#include "stdafx.h"
#include <cmath>
#include <iostream>
#include <fstream>
#include <iomanip>
#include <cstdlib>
using namespace std;



double PressZ(double T, double P, double y1, double y2, double y3)    //压缩因子计算，1为甲烷，2为二氧化碳，3为氮气
{

double Tc[5] = { 0, 190.6, 304.2, 126.2 };
double Pc[5] = { 0, 4.6103e+6, 7.3765e+6, 3.3944e+6 };
double w[5] = { 0, 0.008, 0.225, 0.040 };
double R = 8.3144625;

double ac[5] = { 0 };
double bc[5] = { 0 };
double k[5] = { 0 };
double alpha[5] = { 0 };
double a[5] = { 0 };
double b[5] = { 0 };

int i = 1;
for (i = 1; i<4; i++)

{
ac[i] = 0.45724*R*R*Tc[i] * Tc[i] / Pc[i];

bc[i] = 0.07780*R*Tc[i] / Pc[i];

k[i] = 0.37464 + 1.54226*w[i] - 0.26992*w[i] * w[i];

alpha[i] = (1 + k[i] * (1 - sqrt((float)T / Tc[i])))*(1 + k[i] * (1 - sqrt((float)T / Tc[i])));

a[i] = ac[i] * alpha[i];

b[i] = bc[i];
}

double aa, bb, ax[9] = { 0 };
double kij[9] = { 0, 0.120, 0.0311, 0.0170 };   //ax[i]为气体交互作用系数  aij[1]为甲烷-二氧化碳交互系数，aij[2]为甲烷-氮气交互系数，aij[3]为二氧化碳-氮气交互系数

ax[1] = (1 - kij[1])*sqrt(a[1] * a[2]);
ax[2] = (1 - kij[2])*sqrt(a[1] * a[3]);
ax[3] = (1 - kij[3])*sqrt(a[2] * a[3]);

aa = y1*y1*a[1] + y1*y2*ax[1] + y1*y3*ax[2] + y2*y1*ax[1] + y2*y2*a[2] + y2*y3*ax[3] + y3*y1*ax[2] + y3*y2*ax[3] + y3*y3*a[3];
bb = y1*b[1] + y2*b[2] + y3*b[3];


double A, B, Z1, Z2;
A = aa*P / (R*R*T*T);
B = bb*P / (R*T);

Z2 = 0.9;