teqp 0.22.0
Loading...
Searching...
No Matches
vdW.hpp
Go to the documentation of this file.
1#pragma once
2
3#include "teqp/constants.hpp"
4#include "teqp/types.hpp"
5#include "teqp/exceptions.hpp"
6
7namespace teqp {
8
10class vdWEOS1 {
11private:
12 double a, b;
13public:
15 vdWEOS1(double a, double b) : a(a), b(b) {};
16
18 double get_a() const{ return a; }
19 double get_b() const{ return b; }
20
21 const double Ru = 1.380649e-23 * 6.02214076e23;
22
25 template<class VecType>
26 auto R(const VecType& /*molefrac*/) const { return Ru; }
27
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);
35 }
36
38 double p(double T, double v) {
39 return Ru*T/(v - b) - a/(v*v);
40 }
41};
42
44template <typename NumType>
45class vdWEOS {
46protected:
47 std::valarray<NumType> ai, bi;
48 std::valarray<std::valarray<NumType>> k;
49
50 template<typename TType, typename IndexType>
51 auto get_ai(TType /*T*/, IndexType i) const { return ai[i]; }
52
53 template<typename TType, typename IndexType>
54 auto get_bi(TType /*T*/, IndexType i) const { return bi[i]; }
55
56public:
62 vdWEOS(const std::valarray<NumType>& Tc_K, const std::valarray<NumType>& pc_Pa)
63 {
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");
66 }
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];
72 }
73 k = std::valarray<std::valarray<NumType>>(std::valarray<NumType>(0.0, Tc_K.size()), Tc_K.size());
74 };
75
78 template<typename TType, typename CompType>
79 auto a(TType /*T*/, 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;
85 }
86 }
87 return forceeval(a_);
88 }
89
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];
97 }
98 return forceeval(b_);
99 }
100
101 const NumType Ru = get_R_gas<double>();
102
106 template<class VecType>
107 auto R(const VecType& /*molefrac*/) const {
108 return Ru;
109 }
110
115 template<typename TType, typename RhoType, typename MoleFracType>
116 auto alphar(const TType &T,
117 const RhoType& rho,
118 const MoleFracType &molefrac) const
119 {
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()));
122 }
123 auto Psiminus = -log(1.0 - b(molefrac) * rho);
124 auto Psiplus = rho;
125 auto val = Psiminus - a(T, molefrac) / (Ru * T) * Psiplus;
126 return forceeval(val);
127 }
128};
129
130}; // namespace teqp
A (very) simple implementation of the van der Waals EOS.
Definition vdW.hpp:10
double get_b() const
Definition vdW.hpp:19
double p(double T, double v)
For testing, provide the pressure explicit form of the EOS.
Definition vdW.hpp:38
auto R(const VecType &) const
Get the universal gas constant.
Definition vdW.hpp:26
vdWEOS1(double a, double b)
Intializer, taking the a and b constants directly.
Definition vdW.hpp:15
double get_a() const
Accessor functions.
Definition vdW.hpp:18
const double Ru
Exact value, given by k_B*N_A.
Definition vdW.hpp:21
auto alphar(const TType &T, const RhoType &rhotot, const VecType &molefrac) const
Definition vdW.hpp:33
A slightly more involved implementation of van der Waals, this time with mixture properties.
Definition vdW.hpp:45
auto alphar(const TType &T, const RhoType &rho, const MoleFracType &molefrac) const
The evaluation of .
Definition vdW.hpp:116
vdWEOS(const std::valarray< NumType > &Tc_K, const std::valarray< NumType > &pc_Pa)
Initializer, taking the arrays of critical temperatures and pressures.
Definition vdW.hpp:62
auto get_ai(TType, IndexType i) const
Definition vdW.hpp:51
std::valarray< std::valarray< NumType > > k
Definition vdW.hpp:48
auto get_bi(TType, IndexType i) const
Definition vdW.hpp:54
std::valarray< NumType > bi
Definition vdW.hpp:47
std::valarray< NumType > ai
Definition vdW.hpp:47
const NumType Ru
Universal gas constant, exact number.
Definition vdW.hpp:101
auto a(TType, const CompType &molefracs) const
Calculate the a parameter, based on quadratic mixing rules.
Definition vdW.hpp:79
auto R(const VecType &) const
Get the universal gas constant Here the real universal gas constant, with no composition dependence.
Definition vdW.hpp:107
auto b(const CompType &molefracs) const
Calculate the b parameter, based on linear mixing rules.
Definition vdW.hpp:93
auto pow(const double &x, const double &e)
Definition types.hpp:195
auto forceeval(T &&expr)
Definition types.hpp:52
auto get_R_gas()
< Gas constant, according to CODATA 2019, in the given number type
Definition constants.hpp:22