15 vdWEOS1(
double a,
double b) : a(a), b(b) {};
18 double get_a()
const{
return a; }
19 double get_b()
const{
return b; }
21 const double Ru = 1.380649e-23 * 6.02214076e23;
25 template<
class VecType>
26 auto R(
const VecType& )
const {
return Ru; }
32 template<
typename TType,
typename RhoType,
typename VecType>
33 auto alphar(
const TType &T,
const RhoType& rhotot,
const VecType &molefrac)
const {
34 return forceeval(-log(1.0 - b * rhotot) - (a / (
R(molefrac) * T)) * rhotot);
38 double p(
double T,
double v) {
39 return Ru*T/(v - b) - a/(v*v);
44template <
typename NumType>
47 std::valarray<NumType>
ai,
bi;
48 std::valarray<std::valarray<NumType>>
k;
50 template<
typename TType,
typename IndexType>
51 auto get_ai(TType , IndexType i)
const {
return ai[i]; }
53 template<
typename TType,
typename IndexType>
54 auto get_bi(TType , IndexType i)
const {
return bi[i]; }
62 vdWEOS(
const std::valarray<NumType>& Tc_K,
const std::valarray<NumType>& pc_Pa)
64 if (Tc_K.size() != pc_Pa.size()){
65 throw teqp::InvalidArgument(
"Sizes of Tc_K " + std::to_string(Tc_K.size()) +
" and pc_Pa" + std::to_string(pc_Pa.size()) +
" do not agree");
67 ai.resize(Tc_K.size());
68 bi.resize(Tc_K.size());
69 for (
auto i = 0U; i < Tc_K.size(); ++i) {
70 ai[i] = 27.0 / 64.0 *
pow(
Ru * Tc_K[i], 2) / pc_Pa[i];
71 bi[i] = 1.0 / 8.0 *
Ru * Tc_K[i] / pc_Pa[i];
73 k = std::valarray<std::valarray<NumType>>(std::valarray<NumType>(0.0, Tc_K.size()), Tc_K.size());
78 template<
typename TType,
typename CompType>
79 auto a(TType ,
const CompType& molefracs)
const {
80 typename CompType::value_type a_ = 0.0;
81 for (
auto i = 0; i < molefracs.size(); ++i) {
82 for (
auto j = 0; j < molefracs.size(); ++j) {
83 auto aij = (1 -
k[i][j]) * sqrt(this->ai[i] * this->ai[j]);
84 a_ = a_ + molefracs[i] * molefracs[j] * aij;
92 template<
typename CompType>
93 auto b(
const CompType& molefracs)
const {
94 typename CompType::value_type b_ = 0.0;
95 for (
auto i = 0; i < molefracs.size(); ++i) {
96 b_ = b_ + molefracs[i] *
bi[i];
106 template<
class VecType>
107 auto R(
const VecType& )
const {
115 template<
typename TType,
typename RhoType,
typename MoleFracType>
118 const MoleFracType &molefrac)
const
120 if (
static_cast<std::size_t
>(molefrac.size()) !=
ai.size()) {
121 throw teqp::InvalidArgument(
"mole fractions must be of size " + std::to_string(
ai.size()) +
" but are of size " + std::to_string(molefrac.size()));
123 auto Psiminus = -log(1.0 -
b(molefrac) * rho);
125 auto val = Psiminus -
a(T, molefrac) / (
Ru * T) * Psiplus;
A (very) simple implementation of the van der Waals EOS.
double p(double T, double v)
For testing, provide the pressure explicit form of the EOS.
auto R(const VecType &) const
Get the universal gas constant.
vdWEOS1(double a, double b)
Intializer, taking the a and b constants directly.
double get_a() const
Accessor functions.
const double Ru
Exact value, given by k_B*N_A.
auto alphar(const TType &T, const RhoType &rhotot, const VecType &molefrac) const
A slightly more involved implementation of van der Waals, this time with mixture properties.
auto alphar(const TType &T, const RhoType &rho, const MoleFracType &molefrac) const
The evaluation of .
vdWEOS(const std::valarray< NumType > &Tc_K, const std::valarray< NumType > &pc_Pa)
Initializer, taking the arrays of critical temperatures and pressures.
auto get_ai(TType, IndexType i) const
std::valarray< std::valarray< NumType > > k
auto get_bi(TType, IndexType i) const
std::valarray< NumType > bi
std::valarray< NumType > ai
const NumType Ru
Universal gas constant, exact number.
auto a(TType, const CompType &molefracs) const
Calculate the a parameter, based on quadratic mixing rules.
auto R(const VecType &) const
Get the universal gas constant Here the real universal gas constant, with no composition dependence.
auto b(const CompType &molefracs) const
Calculate the b parameter, based on linear mixing rules.
auto pow(const double &x, const double &e)
auto get_R_gas()
< Gas constant, according to CODATA 2019, in the given number type