| Line | Branch | Exec | Source |
|---|---|---|---|
| 1 | // -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*- | ||
| 2 | // vi: set et ts=4 sw=4 sts=4: | ||
| 3 | // | ||
| 4 | // SPDX-FileCopyrightInfo: Copyright © DuMux Project contributors, see AUTHORS.md in root folder | ||
| 5 | // SPDX-License-Identifier: GPL-3.0-or-later | ||
| 6 | // | ||
| 7 | /*! | ||
| 8 | * \file | ||
| 9 | * \ingroup Components | ||
| 10 | * \brief Shomate equations for enthalpy and heat capacity. | ||
| 11 | * | ||
| 12 | * The Shomate equations were inroduced in \cite Shomate1954. | ||
| 13 | * They offer analytical equations for computing the heat capacity as well as enthalpy. | ||
| 14 | * Both equations use a set of component-specific parameters \f$A,B,C,D,E,F,G,H\f$. | ||
| 15 | * For the heat capacity \f$C_p^0\f$, one obtains | ||
| 16 | \f[ | ||
| 17 | C_p^0 = A + Bt + Ct^2 + Dt^3 + E/t^2, | ||
| 18 | \f] | ||
| 19 | * while for the enthalpy with a reference state at \f$T=298.15\text{K}\f$, one uses | ||
| 20 | \f[ | ||
| 21 | H^0-H_{298.15}^0 = At + \frac{Bt^2}{2} + \frac{Ct^3}{3} + \frac{Dt^4}{4} -\frac{E}{t} + F - H | ||
| 22 | \f] | ||
| 23 | *where: | ||
| 24 | * * \f$ C_p \f$ is the heat capacity in J/(mol*K), | ||
| 25 | * * \f$ H^0 \f$ represents the standard enthalpy in kJ/mol, | ||
| 26 | * * \f$ t \f$ is the temperature in K divided by 1000. | ||
| 27 | * | ||
| 28 | */ | ||
| 29 | #ifndef DUMUX_MATERIAL_COMPONENTS_SHOMATE_HH | ||
| 30 | #define DUMUX_MATERIAL_COMPONENTS_SHOMATE_HH | ||
| 31 | |||
| 32 | #include <algorithm> | ||
| 33 | #include <array> | ||
| 34 | #include <iostream> | ||
| 35 | #include <type_traits> | ||
| 36 | #include <vector> | ||
| 37 | |||
| 38 | #include <dune/common/exceptions.hh> | ||
| 39 | |||
| 40 | namespace Dumux { | ||
| 41 | |||
| 42 | /*! | ||
| 43 | * \brief The Shomate method to compute enthalpy and heat capacity | ||
| 44 | * \tparam Scalar the type of the scalar values | ||
| 45 | * \tparam intervals the static number of intervals for the temperature range, each interval has its own set of coefficients. | ||
| 46 | * If -1 the number of intervals is dynamic. | ||
| 47 | */ | ||
| 48 | template<class Scalar, int intervals = -1> | ||
| 49 | class ShomateMethod | ||
| 50 | { | ||
| 51 | static_assert(intervals == -1 || intervals > 0, "Number of intervals must be -1 (dynamic) or greater than 0 (static)."); | ||
| 52 | |||
| 53 | public: | ||
| 54 | using CoefficientSet = struct { Scalar A, B, C, D, E, F, G, H; }; | ||
| 55 | using Temperatures = std::conditional_t<intervals == -1, std::vector<Scalar>, std::array<Scalar, std::size_t(intervals+1)>>; | ||
| 56 | using Coefficients = std::conditional_t<intervals == -1, std::vector<CoefficientSet>, std::array<CoefficientSet, std::size_t(intervals)>>; | ||
| 57 | |||
| 58 | /*! | ||
| 59 | * \brief Constructor for Shomate method class | ||
| 60 | * \param temperatures lower bound of the temperature intervals plus the upper bound of the last interval | ||
| 61 | * \param coeffs contains the sets of coefficients for each temperature interval bounded by the entries of temperatures | ||
| 62 | */ | ||
| 63 | constexpr ShomateMethod(const Temperatures& temperatures, const Coefficients& coeffs) | ||
| 64 | : temperatures_(temperatures) | ||
| 65 | , coeffs_(coeffs) | ||
| 66 | { | ||
| 67 | checkInput_(); | ||
| 68 | } | ||
| 69 | |||
| 70 | /*! | ||
| 71 | * \brief Return enthalpy in kJ/mol | ||
| 72 | * \param temperature in K | ||
| 73 | * \return Scalar | ||
| 74 | */ | ||
| 75 | 21186516 | Scalar enthalpy(const Scalar temperature) const | |
| 76 | { | ||
| 77 | 21186516 | const auto& p = paramsAtTemperature_(temperature); | |
| 78 | 21186516 | const Scalar t = temperature/1000.0; | |
| 79 | 21186516 | const Scalar standardEnthalpyDiff = p.A*t + p.B*t*t/2.0 + p.C*t*t*t/3.0 + p.D*t*t*t*t/4.0 - p.E/t + p.F - p.H; | |
| 80 | 21186516 | return standardEnthalpyDiff; | |
| 81 | } | ||
| 82 | |||
| 83 | /*! | ||
| 84 | * \brief Return heat capacity in J/(mol*K) | ||
| 85 | * \param temperature in K | ||
| 86 | * \return Scalar | ||
| 87 | */ | ||
| 88 | 1047 | Scalar heatCapacity(const Scalar temperature) const | |
| 89 | { | ||
| 90 | 1047 | const auto& p = paramsAtTemperature_(temperature); | |
| 91 | 1047 | const Scalar t = temperature/1000.0; | |
| 92 | 1047 | const Scalar heatCapacity = p.A + p.B*t + p.C*t*t + p.D*t*t*t + p.E/(t*t); | |
| 93 | 1047 | return heatCapacity; | |
| 94 | } | ||
| 95 | |||
| 96 | private: | ||
| 97 | 21187563 | const CoefficientSet& paramsAtTemperature_(const Scalar T) const | |
| 98 | { | ||
| 99 | // check if T is smaller or higher than allowed min/max T | ||
| 100 |
8/8✓ Branch 0 taken 21186235 times.
✓ Branch 1 taken 318 times.
✓ Branch 2 taken 21186235 times.
✓ Branch 3 taken 318 times.
✓ Branch 4 taken 102 times.
✓ Branch 5 taken 21186133 times.
✓ Branch 6 taken 102 times.
✓ Branch 7 taken 21186133 times.
|
42375126 | if (T < temperatures_.front() || T > temperatures_.back()) |
| 101 | { | ||
| 102 |
2/2✓ Branch 0 taken 4 times.
✓ Branch 1 taken 416 times.
|
822 | if (!warningPrinted_) |
| 103 | { | ||
| 104 | 21 | std::cout << "Temperature "<< T << " [K] is out of range. Enthalpy values are extrapolated." << std::endl; | |
| 105 | 7 | warningPrinted_ = true; | |
| 106 | } | ||
| 107 | } | ||
| 108 | |||
| 109 | // find the interval for the given temperature | ||
| 110 | 42375126 | const auto index = std::min<std::size_t>( | |
| 111 | 42375126 | coeffs_.size()-1, | |
| 112 |
6/6✓ Branch 0 taken 21161619 times.
✓ Branch 1 taken 24934 times.
✓ Branch 2 taken 21161619 times.
✓ Branch 3 taken 24934 times.
✓ Branch 4 taken 21161619 times.
✓ Branch 5 taken 24934 times.
|
127125378 | std::distance( |
| 113 | temperatures_.begin(), | ||
| 114 | std::lower_bound(temperatures_.begin(), temperatures_.end(), T) | ||
| 115 | ) | ||
| 116 | ); | ||
| 117 | |||
| 118 | 42375126 | return coeffs_[index]; | |
| 119 | } | ||
| 120 | |||
| 121 | ✗ | void checkInput_() const | |
| 122 | { | ||
| 123 | if constexpr (intervals == -1) | ||
| 124 | { | ||
| 125 | if (temperatures_.size() < 2) | ||
| 126 | DUNE_THROW(Dune::InvalidStateException, "Temperature range must have at least two entries."); | ||
| 127 | |||
| 128 | for (size_t i = 0; i < temperatures_.size()-1; i++) | ||
| 129 | if (temperatures_[i] >= temperatures_[i+1]) | ||
| 130 | DUNE_THROW(Dune::InvalidStateException, "Temperature range must be strictly increasing."); | ||
| 131 | |||
| 132 | if (temperatures_.size()-1 != coeffs_.size()) | ||
| 133 | DUNE_THROW(Dune::InvalidStateException, "If temperature vector is size n+1, there must be n coefficient sets."); | ||
| 134 | } | ||
| 135 | ✗ | } | |
| 136 | |||
| 137 | Temperatures temperatures_; | ||
| 138 | Coefficients coeffs_; | ||
| 139 | static bool warningPrinted_; | ||
| 140 | }; | ||
| 141 | |||
| 142 | template<class Scalar, int intervals> | ||
| 143 | bool ShomateMethod<Scalar, intervals>::warningPrinted_ = false; | ||
| 144 | |||
| 145 | } // namespace Dumux | ||
| 146 | |||
| 147 | #endif | ||
| 148 |