From 7902eb10fee1a69f6288d583897f480f8b59aa2f Mon Sep 17 00:00:00 2001
From: Ian Bell <ian.bell@nist.gov>
Date: Thu, 18 Mar 2021 13:19:13 -0400
Subject: [PATCH] Fix second virial coefficient calculation
Refactor what is implemented
---
include/teqp/derivs.hpp | 14 +++-----------
include/teqp/models/multifluid.hpp | 11 +++++++++--
2 files changed, 12 insertions(+), 13 deletions(-)
diff --git a/include/teqp/derivs.hpp b/include/teqp/derivs.hpp
index 218f669..9fb166d 100644
--- a/include/teqp/derivs.hpp
+++ b/include/teqp/derivs.hpp
@@ -97,17 +97,9 @@ typename ContainerType::value_type get_Ar01(const Model& model, const TType T, c
template <typename Model, typename TType, typename ContainerType>
typename ContainerType::value_type get_B2vir(const Model& model, const TType T, const ContainerType& molefrac) {
- auto rhovec = 0.0*molefrac;
- double max_mole_frac = *std::max_element(std::begin(molefrac), std::end(molefrac));
- if (max_mole_frac != 1.0) {
- std::cout << "The value for B_2 is incorrect when not pure!" << std::endl;
- }
- decltype(molefrac[0] * T) B2 = 0.0;
- for (auto i = 0; i < rhovec.size(); ++i) {
- auto dalphar_drhoi__constTrhoj = derivrhoi([&model](const auto& T, const auto& rhovec) { return model.alphar(T, rhovec); }, T, rhovec, i);
- B2 += molefrac[i] * dalphar_drhoi__constTrhoj;
- }
-
+ double h = 1e-100;
+ // B_2 = lim_rho\to 0 dalphar/drho|T,z
+ auto B2 = model.alphar(T, std::complex<double>(0, h), molefrac).imag()/h;
return B2;
}
diff --git a/include/teqp/models/multifluid.hpp b/include/teqp/models/multifluid.hpp
index eeafd62..ee469d5 100644
--- a/include/teqp/models/multifluid.hpp
+++ b/include/teqp/models/multifluid.hpp
@@ -84,11 +84,18 @@ public:
{
RhoType::value_type rhotot_ = (rhotot.has_value()) ? rhotot.value() : std::accumulate(std::begin(rhovec), std::end(rhovec), (decltype(rhovec[0]))0.0);
auto molefrac = rhovec / rhotot_;
+ return alphar(T, rhotot_, molefrac);
+ }
+
+ template<typename TType, typename RhoType, typename MoleFracType>
+ auto alphar(TType T,
+ const RhoType rho,
+ const MoleFracType& molefrac) const
+ {
auto Tred = redfunc.get_Tr(molefrac);
auto rhored = redfunc.get_rhor(molefrac);
- auto delta = rhotot_ / rhored;
+ auto delta = rho / rhored;
auto tau = Tred / T;
- using resulttype = decltype(T*rhovec[0]);
return corr.alphar(tau, delta, molefrac) + dep.alphar(tau, delta, molefrac);
}
};
--
GitLab