Newer
Older
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
#pragma once
#include "teqp/models/multifluid.hpp"
#include <Eigen/Dense>
namespace teqp {
namespace Mie{
/**
Equation of state for the Mie ($\lambda_{\mathrm{r}}$,6) fluid with a repulsive exponent from 11 to 13
Pohl, S.; Fingerhut, R., Thol, M.; Vrabec, J.; Span, R.
*/
class Mie6Pohl2023 {
private:
using EArray6 = Eigen::Array<double, 6, 1>;
using EArray4 = Eigen::Array<double, 4, 1>;
const EArray6 c1_pol = (EArray6() << -0.0192944,1.38,-2.2653,1.6291,-1.974,0.40412 ).finished();
const EArray6 c1_exp = (EArray6() << 0.1845,-0.3227,1.1351,2.232,-2.344,-0.4238 ).finished();
const EArray4 c1_gbs = (EArray4() << -4.367,0.0371,1.3895,2.835 ).finished();
const EArray6 c2_pol = (EArray6() << 0.26021,-5.525,8.329,-19.492,25.8,-3.8133 ).finished();
const EArray6 c2_exp = (EArray6() << -5.05,2.7842,-9.523,-30.383,17.902,2.2264 ).finished();
const EArray4 c2_gbs = (EArray4() << 48.445,-5.506,-11.643,-24.36 ).finished();
const EArray6 t_pol = (EArray6() << 1,0.236,0.872,0.313,0.407,0.703 ).finished();
const EArray6 t_exp = (EArray6() << 1.78,2.99,2.866,1.2,3.06,1.073 ).finished();
const EArray4 t_gbs = (EArray4() << 1.50,1.03,4.02,1.57 ).finished();
const EArray6 d_pol = (EArray6() << 4,1,1,2,2,3 ).finished();
const EArray6 d_exp = (EArray6() << 1,1,3,2,2,5 ).finished();
const EArray4 d_gbs = (EArray4() << 2,3,2,2 ).finished();
const EArray6 p = (EArray6() << 1,2,2,1,2,1 ).finished();
const EArray4 eta = (EArray4() << 0.362,0.313,1.17,0.957 ).finished();
const EArray4 beta = (EArray4() << 0.0761,0.143,0.63,1.32 ).finished();
const EArray4 gam = (EArray4() << 1.55,-0.0826,1.505,1.07 ).finished();
const EArray4 eps = (EArray4() << -1,-1,-0.195,-0.287 ).finished();
const double m_lambda_a;
const EArray6 n_pol, n_exp;
const EArray4 n_gbs;
const double Tc, rhoc; // In simulation units
public:
Mie6Pohl2023(double lambda_a) : m_lambda_a(lambda_a),
n_pol(c1_pol + c2_pol / m_lambda_a),
n_exp(c1_exp + c2_exp / m_lambda_a),
n_gbs(c1_gbs + c2_gbs / m_lambda_a),
Tc(0.668 + 6.84 / m_lambda_a + 145 / pow(m_lambda_a, 3)), // T^*
rhoc(0.2516 + 0.049 * log10(m_lambda_a)) // rho^*
{}
auto get_lambda_a() const { return m_lambda_a; }
// We are in "simulation units", so R is 1.0, and T and rho that
// go into alphar are actually T^* and rho^*
template<typename MoleFracType>
double R(const MoleFracType &) const { return 1.0; }
template<typename TTYPE, typename RHOTYPE, typename MoleFracType>
auto alphar(const TTYPE& Tstar, const RHOTYPE& rhostar, const MoleFracType& /*molefrac*/) const {
auto tau = Tc / Tstar; auto delta = rhostar / rhoc;
auto alphar_ = (
(n_pol * pow(tau, t_pol) * pow(delta, d_pol)).sum() +
(n_exp * pow(tau, t_exp) * pow(delta, d_exp) * exp(-pow(delta, p))).sum() +
(n_gbs * pow(tau, t_gbs) * pow(delta, d_gbs) * exp(-eta * pow(delta - eps, 2) - beta * pow(tau - gam, 2))).sum()
);
return forceeval(alphar_);
}
};
};
};