olb::MultiComponentPengRobinson Class Reference

#include <interactionPotential.h>

Collaboration diagram for olb::MultiComponentPengRobinson:

Public Member Functions
	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)
double	G_Excess (std::vector< double > x_i, std::vector< double > g_jk, std::vector< double > beta_jk)
double	b_mix (std::vector< double > x_i)
double	a_mix (std::vector< double > x_i, std::vector< double > a_i, double b, double G_E)
double	p_PR (double v, double a_m, double b_m)
double	gamma_i (std::vector< double > x_i, double a_i, std::vector< double > g_jk, std::vector< double > beta_jk, double _T, int i)
double	lnfugacity_i (std::vector< double > x_i, double v, double a_m, double b_m, double gamma_i, double _T, int i)
std::vector< double >	iterate_VLE (double residuum, double beta0)
std::vector< double >	getChis (int n)

Detailed Description

Definition at line 38 of file interactionPotential.h.

Constructor & Destructor Documentation

MultiComponentPengRobinson()

olb::MultiComponentPengRobinson::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 )

Definition at line 191 of file interactionPotential.hh.

                                                                                                                                  :
p(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)
{ }

Member Function Documentation

a_mix()

double olb::MultiComponentPengRobinson::a_mix	(	std::vector< double >	x_i,
		std::vector< double >	a_i,
		double	b,
		double	G_E )

Definition at line 216 of file interactionPotential.hh.

                                                                                                            {
  double a_sum = 0;
  for(int i = 0; i < num_comp; i++){
    a_sum = a_sum + x_i[i]*a_i[i]/b[i];
  }
  /*double a_m = 0;
  for(int i = 0; i < num_comp; i++){
    a_i[i] = a_c[i] * util::pow((1 + m[i] * (1 - util::pow((_T/T_c[i]),0.5))),2);
    for(int j = 0; j < num_comp; j++){
      a_m += x_i[i]*x_i[j]*util::pow(a_i[i]*a_i[j],0.5);
    }
  }
  return a_m;*/
  return b_m * (a_sum - G_E / lambda);
}

Here is the caller graph for this function:

b_mix()

double olb::MultiComponentPengRobinson::b_mix

(

std::vector< double >

x_i

)

Definition at line 209 of file interactionPotential.hh.

                                                             {
  double b_m = 0;
  for(int i = 0; i < num_comp; i++){
    b_m = b_m + b[i]*x_i[i];
  }
  return b_m;
}

Here is the caller graph for this function:

G_Excess()

double olb::MultiComponentPengRobinson::G_Excess	(	std::vector< double >	x_i,
		std::vector< double >	g_jk,
		std::vector< double >	beta_jk )

Definition at line 196 of file interactionPotential.hh.

                                                                                                                   {
  double G_E = 0;
  for(int i = 0; i < num_comp; i++){
    double sum1 = 0;
    double sum2 = 0;
    for(int j = 0; j < num_comp; j++){
      sum1 = sum1 + x_i[j]*g_jk[i+num_comp*j]*beta_jk[i+num_comp*j];
      sum2 = sum2 + x_i[j]*beta_jk[i+num_comp*j];
    }
    G_E = G_E + x_i[i] * sum1/sum2;
  }
  return G_E;
}

Here is the caller graph for this function:

gamma_i()

double olb::MultiComponentPengRobinson::gamma_i	(	std::vector< double >	x_i,
		double	a_i,
		std::vector< double >	g_jk,
		std::vector< double >	beta_jk,
		double	_T,
		int	i )

Definition at line 234 of file interactionPotential.hh.

                                                                                                                                                {
  double sum1 = 0;
  double sum2 = 0;
  double sum3 = 0;
  double sum4 = 0;
  for(int j = 0; j < num_comp; j++){
    sum1 = sum1 + x_i[j]*g_jk[num_comp*j+i]*beta_jk[num_comp*j+i];
    sum2 = sum2 + x_i[j]*beta_jk[num_comp*j+i];
    double sub31 = 0;
    double sub41 = 0;
    double sub42 = 0;
    for(int k = 0; k < num_comp; k++){
      sub31 = sub31 + x_i[k]*beta_jk[j+num_comp*k];
      sub41 = sub41 + x_i[k]*g_jk[j+num_comp*k]*beta_jk[j+num_comp*k];
      sub42 = sub42 + x_i[k]*beta_jk[j+num_comp*k];
    }
    sum3 = sum3 + x_i[j]*g_jk[num_comp*i+j]*beta_jk[num_comp*i+j] / sub31;
    sum4 = sum4 + x_i[j]*beta_jk[num_comp*i+j]*sub41/util::pow(sub42,2);
  }
  return 1/(R*_T) * (a_i / b[i] - 1/lambda * (sum1/sum2 + sum3 - sum4));
}

References olb::util::pow().

Here is the call graph for this function:

Here is the caller graph for this function:

getChis()

std::vector< double > olb::MultiComponentPengRobinson::getChis

(

int

n

)

Definition at line 369 of file interactionPotential.hh.

                                                          {
  std::vector<double> n_chis(n);
  for(int i = 0; i<n; i++){
    n_chis[i] = chis[i];
  }
  return n_chis;
}

iterate_VLE()

std::vector< double > olb::MultiComponentPengRobinson::iterate_VLE	(	double	residuum,
		double	beta0 )

Definition at line 275 of file interactionPotential.hh.

                                                                                      {
  std::vector<double> a_i(num_comp), g_jk(num_comp*num_comp), beta_jk(num_comp*num_comp);
  for(int i = 0; i < num_comp; i++){
    a_i[i] = a_c[i] * util::pow((1 + m[i] * (1 - util::pow((T/T_c[i]),0.5))),2);
    for(int j = 0; j < num_comp; j++){
      g_jk[num_comp*i+j] = g_I[num_comp*i+j] + g_II[num_comp*i+j] * T;
      beta_jk[num_comp*i+j] = b[j] * util::exp(-alpha[num_comp*i+j]*g_jk[num_comp*i+j]/(R*T));
    }
  }
  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);
  std::vector<double> deltaLogFugacity(num_comp), errorFugacity(num_comp);
  double totalError = 0;
  for(int i = 0; i < num_comp; i++){
    //adjust prefactor or initial values for partition coefficients, 0.9 for hydrocarbons
    K[i] = 1.0*(p_c[i]/(p))*util::exp(5.373*(1+omega[i])*(1-T_c[i]/T));
    errorFugacity[i] = 1;
  }
  for(int i = 0; i < 2; i++){
    vx[i] = init_v[i];
  }
  double beta = beta0;
  double f = 0;
  do{
    totalError = 0;
    f = 0;
    for(int i = 0; i < num_comp; i++){
      f += (K[i]-1)*z[i]/(1+(K[i]-1)*beta);
    }
    while(abs(f) > residuum){
      double dfdbeta = 0;
      for(int i = 0; i < num_comp; i++){
        dfdbeta += -1*(K[i]-1)*(K[i]-1)*z[i]/((1+(K[i]-1)*beta)*(1+(K[i]-1)*beta));
      }
      beta = beta - f / dfdbeta;
      f = 0;
      for(int i = 0; i < num_comp; i++){
        f += (K[i]-1)*z[i]/(1+(K[i]-1)*beta);
      }
    }
    double x_sum = 0;
    double y_sum = 0;
    for(int i = 0; i < num_comp; i++){
      x[i] = z[i]/(1+(K[i]-1)*beta);
      y[i] = K[i]*x[i];
      x_sum += x[i];
      y_sum += y[i];
    }
    if(x_sum != 1 || y_sum != 1){
      for(int i = 0; i < num_comp; i++){
        vx[i+2] = x[i]/x_sum;
        x[i] = vx[i+2];
        vx[i+2+num_comp] = y[i]/y_sum;
        y[i] = vx[i+2+num_comp];
      }
    }
 
    G_E[0] = G_Excess(x, g_jk, beta_jk);
    G_E[1] = G_Excess(y, g_jk, beta_jk);
    b_m[0] = b_mix(x);
    b_m[1] = b_mix(y);
    a_m[0] = a_mix(x, a_i, b_m[0], G_E[0]);
    a_m[1] = a_mix(y, a_i, b_m[1], G_E[1]);
    for(int i = 0; i < num_comp; i++){
      gammas[i] = gamma_i(x, a_i[i], g_jk, beta_jk, T, i);
      gammas[i+num_comp] = gamma_i(y, a_i[i], g_jk, beta_jk, T, i);
    }
 
    for(int i = 0; i < 2; i++){
      while(abs(p-p_PR(vx[i], a_m[i], b_m[i])) > residuum){
        double dpdv = -R*T/util::pow(vx[i]-b_m[i],2) +
                2*a_m[i]*(vx[i]+b_m[i])/util::pow(util::pow(vx[i],2) + 2*vx[i]*b_m[i] - util::pow(b_m[i],2),2);
        vx[i] = vx[i] - (p - p_PR(vx[i], a_m[i], b_m[i])) / (-dpdv);
      }
    }
    for(int i = 0; i < num_comp; i++){
      deltaLogFugacity[i] = lnfugacity_i(x, vx[0], a_m[0], b_m[0], gammas[i], T, i) -
              lnfugacity_i(y, vx[1], a_m[1], b_m[1], gammas[i+num_comp], T, i);
      errorFugacity[i] = util::exp(deltaLogFugacity[i]) - 1;
      K[i] = util::exp(deltaLogFugacity[i]) / (x[i]/y[i]);
      totalError += util::pow(errorFugacity[i],2);
    }
  }while(util::pow(totalError,0.5) > residuum);
  for(int i = 1; i<num_comp; i++){
    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]));
  }
  double vol_sum = 0;
  for(int i = 0; i<num_comp; i++){
    vol_sum += volatilities[i];
  }
  for(int i = 0; i<num_comp; i++){
    chis[i] = volatilities[i] / vol_sum;
  }
  return vx;
}

References a_mix(), olb::abs(), b_mix(), olb::util::exp(), G_Excess(), gamma_i(), lnfugacity_i(), p_PR(), and olb::util::pow().

Here is the call graph for this function:

lnfugacity_i()

double olb::MultiComponentPengRobinson::lnfugacity_i	(	std::vector< double >	x_i,
		double	v,
		double	a_m,
		double	b_m,
		double	gamma_i,
		double	_T,
		int	i )

Definition at line 255 of file interactionPotential.hh.

                                                                                                                                        {
  //TODO change inputs for non-isothermal cases!!!
  double term1 = b[i]/(v-b_m) + util::log(x_i[i]*R*_T/(v-b_m));
  double term2 = -gamma_i / util::pow(2,1.5) * util::log((v + (util::pow(2,0.5)+1)*b_m)/(v-(util::pow(2,0.5)-1)*b_m));
  double term3 = -b[i]*a_m/(util::pow(2,1.5)*b_m*R*_T)*
          ((util::pow(2,0.5)+1)/(v+(util::pow(2,0.5)+1)*b_m)+(util::pow(2,0.5)-1)/(v-(util::pow(2,0.5)-1)*b_m));
  /*std::vector<double> a_i(num_comp);
  a_m = 0;
  for(int i = 0; i < num_comp; i++){
    a_i[i] = a_c[i] * util::pow((1 + m[i] * (1 - util::pow((_T/T_c[i]),0.5))),2);
    for(int j = 0; j < num_comp; j++){
      a_m += x_i[i]*x_i[j]*util::pow(a_i[i]*a_i[j],0.5);
    }
  }
  gamma_i = -a_m/(b_m*R*_T)*(b[i]/b_m-2/a_m*(x_i[0]*util::pow(a_i[0]*a_i[i],0.5)+x_i[1]*util::pow(a_i[1]*a_i[i],0.5)));
  term1 = b[i]/(v-b_m) + util::log(x_i[i]*R*_T/(v-b_m));
  term2 = -gamma_i / util::pow(2,1.5) * util::log((v + (util::pow(2,0.5)+1)*b_m)/(v-(util::pow(2,0.5)-1)*b_m));
  term3 = -b[i]*a_m/(util::pow(2,1.5)*b_m*R*_T)*((util::pow(2,0.5)+1)/(v+(util::pow(2,0.5)+1)*b_m)+(util::pow(2,0.5)-1)/(v-(util::pow(2,0.5)-1)*b_m));*/
  return term1 + term2 + term3;
}

References gamma_i(), olb::util::log(), and olb::util::pow().

Here is the call graph for this function:

Here is the caller graph for this function:

p_PR()

double olb::MultiComponentPengRobinson::p_PR	(	double	v,
		double	a_m,
		double	b_m )

Definition at line 231 of file interactionPotential.hh.

                                                                       {
  return R*T/(v - b_m) - a_m/(util::pow(v,2) + 2*v*b_m - util::pow(b_m,2));
}

References olb::util::pow().

Here is the call graph for this function:

Here is the caller graph for this function:

---

The documentation for this class was generated from the following files:

• src/dynamics/interactionPotential.h
• src/dynamics/interactionPotential.hh