| 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 FreeflowModels | ||
| 10 | * \brief This file contains different functions for estimating turbulence properties. | ||
| 11 | */ | ||
| 12 | |||
| 13 | #ifndef DUMUX_TURBULENCE_PROPERTIES_HH | ||
| 14 | #define DUMUX_TURBULENCE_PROPERTIES_HH | ||
| 15 | |||
| 16 | #include <iostream> | ||
| 17 | |||
| 18 | #include<dune/common/fvector.hh> | ||
| 19 | |||
| 20 | namespace Dumux { | ||
| 21 | |||
| 22 | /*! | ||
| 23 | * \brief This class contains different functions for estimating turbulence properties. | ||
| 24 | */ | ||
| 25 | template<class Scalar, unsigned dim, bool verbose = false> | ||
| 26 | class TurbulenceProperties | ||
| 27 | { | ||
| 28 | public: | ||
| 29 | /*! | ||
| 30 | * \brief Estimates dimensionless wall distance \f$ y^+ \f$ based on a formula given in | ||
| 31 | * http://www.cfd-online.com/Wiki/Y_plus_wall_distance_estimation | ||
| 32 | */ | ||
| 33 | Scalar yPlusEstimation(const Scalar velocity, | ||
| 34 | const Dune::FieldVector<Scalar, dim> position, | ||
| 35 | const Scalar kinematicViscosity, | ||
| 36 | const Scalar density, | ||
| 37 | int yCoordDim=dim-1, | ||
| 38 | bool print=verbose) const | ||
| 39 | { | ||
| 40 | using std::pow; | ||
| 41 | using std::log10; | ||
| 42 | using std::sqrt; | ||
| 43 | const Scalar re_x = reynoldsNumber(velocity, position[0], kinematicViscosity, false); | ||
| 44 | const Scalar c_f = pow((2.0 * log10(re_x) - 0.65), -2.3); // for re_x < 10^9 | ||
| 45 | const Scalar wallShearStress = 0.5 * c_f * density * velocity * velocity; | ||
| 46 | const Scalar frictionVelocity = sqrt(wallShearStress / density); | ||
| 47 | const Scalar yPlus = position[yCoordDim] * frictionVelocity / kinematicViscosity; | ||
| 48 | if (print) | ||
| 49 | { | ||
| 50 | std::cout << "turbulence properties at ("; | ||
| 51 | for (unsigned int dimIdx = 0; dimIdx < dim; ++dimIdx) | ||
| 52 | std::cout << position[dimIdx] << ","; | ||
| 53 | std::cout << ")" << std::endl; | ||
| 54 | std::cout << "estimated Re_x : " << re_x << " [-]" << std::endl; | ||
| 55 | std::cout << "estimated c_f : " << c_f << " [-]" << std::endl; | ||
| 56 | std::cout << "estimated tau_w : " << wallShearStress << " [kg/(m*s^2)]" << std::endl; | ||
| 57 | std::cout << "estimated UStar : " << frictionVelocity << " [m/s]" << std::endl; | ||
| 58 | std::cout << "estimated yPlus : " << yPlus << " [-]" << std::endl; | ||
| 59 | std::cout << std::endl; | ||
| 60 | } | ||
| 61 | return yPlus; | ||
| 62 | } | ||
| 63 | |||
| 64 | /*! | ||
| 65 | * \brief Estimates the entrance length for this pipe | ||
| 66 | */ | ||
| 67 | Scalar entranceLength(const Scalar velocity, | ||
| 68 | const Scalar diameter, | ||
| 69 | const Scalar kinematicViscosity, | ||
| 70 | bool print=verbose) const | ||
| 71 | { | ||
| 72 | using std::pow; | ||
| 73 | const Scalar re_d = reynoldsNumber(velocity, diameter, kinematicViscosity, false); | ||
| 74 | const Scalar entranceLength = 4.4 * pow(re_d, 1.0/6.0) * diameter; | ||
| 75 | if (print) | ||
| 76 | { | ||
| 77 | std::cout << "estimated Re_d : " << re_d << " [-]" << std::endl; | ||
| 78 | std::cout << "estimated l_ent : " << entranceLength << " [m]"<< std::endl; | ||
| 79 | std::cout << std::endl; | ||
| 80 | } | ||
| 81 | return entranceLength; | ||
| 82 | } | ||
| 83 | |||
| 84 | /*! | ||
| 85 | * \brief Calculates the Reynolds number | ||
| 86 | */ | ||
| 87 | Scalar reynoldsNumber(const Scalar velocity, | ||
| 88 | const Scalar charLengthScale/*e.g. diameter*/, | ||
| 89 | const Scalar kinematicViscosity, | ||
| 90 | bool print=verbose) const | ||
| 91 | { | ||
| 92 | using std::abs; | ||
| 93 | 180 | return abs(velocity * charLengthScale / kinematicViscosity); | |
| 94 | } | ||
| 95 | |||
| 96 | /*! | ||
| 97 | * \brief Estimates the turbulence intensity based on a formula given | ||
| 98 | * in the ANSYS Fluent user guide \cite ANSYSUserGuide12 | ||
| 99 | */ | ||
| 100 | 150 | Scalar turbulenceIntensity(const Scalar reynoldsNumber, | |
| 101 | bool print=verbose) const | ||
| 102 | { | ||
| 103 | using std::pow; | ||
| 104 | 150 | const Scalar turbulenceIntensity = 0.16 * pow(reynoldsNumber, -0.125); | |
| 105 |
2/2✓ Branch 0 taken 30 times.
✓ Branch 1 taken 120 times.
|
150 | if (print) |
| 106 | { | ||
| 107 | 90 | std::cout << "estimated I : " << turbulenceIntensity << " [-]" << std::endl; | |
| 108 | } | ||
| 109 | 150 | return turbulenceIntensity; | |
| 110 | } | ||
| 111 | |||
| 112 | /*! | ||
| 113 | * \brief Estimates the turbulence length scale based on a formula given | ||
| 114 | * in the ANSYS Fluent user guide \cite ANSYSUserGuide12 | ||
| 115 | */ | ||
| 116 | 60 | Scalar turbulenceLengthScale(const Scalar charLengthScale/*e.g. diameter*/, | |
| 117 | bool print=verbose) const | ||
| 118 | { | ||
| 119 | 60 | const Scalar turbulenceLengthScale = 0.07 * charLengthScale; | |
| 120 |
2/2✓ Branch 0 taken 30 times.
✓ Branch 1 taken 30 times.
|
60 | if (print) |
| 121 | { | ||
| 122 | 90 | std::cout << "estimated l_turb: " << turbulenceLengthScale << " [m]" << std::endl; | |
| 123 | } | ||
| 124 | 60 | return turbulenceLengthScale; | |
| 125 | } | ||
| 126 | |||
| 127 | /*! | ||
| 128 | * \brief Estimates the turbulent kinetic energy based on a formula given | ||
| 129 | * in the ANSYS Fluent user guide \cite ANSYSUserGuide12 | ||
| 130 | */ | ||
| 131 | 60 | Scalar turbulentKineticEnergy(const Scalar velocity, | |
| 132 | const Scalar diameter, | ||
| 133 | const Scalar kinematicViscosity, | ||
| 134 | bool print=verbose) const | ||
| 135 | { | ||
| 136 | 60 | const Scalar re_d = reynoldsNumber(velocity, diameter, kinematicViscosity, false); | |
| 137 | 120 | const Scalar k = 1.5 * velocity * velocity | |
| 138 | 60 | * turbulenceIntensity(re_d, false) * turbulenceIntensity(re_d, false); | |
| 139 |
2/2✓ Branch 0 taken 30 times.
✓ Branch 1 taken 30 times.
|
60 | if (print) |
| 140 | { | ||
| 141 | 90 | std::cout << "estimated k : " << k << " [m^2/s^2]" << std::endl; | |
| 142 | } | ||
| 143 | 60 | return k; | |
| 144 | } | ||
| 145 | |||
| 146 | /*! | ||
| 147 | * \brief Estimates the dissipation based on a formula given | ||
| 148 | * in the ANSYS Fluent user guide \cite ANSYSUserGuide12 | ||
| 149 | */ | ||
| 150 | 20 | Scalar dissipation(const Scalar velocity, | |
| 151 | const Scalar diameter, | ||
| 152 | const Scalar kinematicViscosity, | ||
| 153 | bool print=verbose) const | ||
| 154 | { | ||
| 155 | using std::pow; | ||
| 156 | 20 | const Scalar k = turbulentKineticEnergy(velocity, diameter, kinematicViscosity, false); | |
| 157 | 20 | const Scalar factor = 0.1643; // = cMu^(3/4) = 0.09^(3/4) | |
| 158 | 20 | const Scalar epsilon = factor * pow(k, 1.5) / turbulenceLengthScale(diameter, false); | |
| 159 |
1/2✓ Branch 0 taken 20 times.
✗ Branch 1 not taken.
|
20 | if (print) |
| 160 | { | ||
| 161 | 60 | std::cout << "estimated eps. : " << epsilon << " [m^2/s^3]" << std::endl; | |
| 162 | } | ||
| 163 | 20 | return epsilon; | |
| 164 | } | ||
| 165 | |||
| 166 | /*! | ||
| 167 | * \brief Estimates the dissipation rate based on a formula given | ||
| 168 | * in the ANSYS Fluent user guide \cite ANSYSUserGuide12 | ||
| 169 | * \f[ \omega = \frac{k^{1/2}}{C_{\mu}^{1/4}L} \f] | ||
| 170 | */ | ||
| 171 | 10 | Scalar dissipationRate(const Scalar velocity, | |
| 172 | const Scalar diameter, | ||
| 173 | const Scalar kinematicViscosity, | ||
| 174 | bool print=verbose) const | ||
| 175 | { | ||
| 176 | using std::pow; | ||
| 177 | 10 | const Scalar k = turbulentKineticEnergy(velocity, diameter, kinematicViscosity, false); | |
| 178 | 10 | const Scalar factor = 0.54772; // = cMu^(1/4) = 0.09^(1/4) | |
| 179 | 10 | const Scalar L = turbulenceLengthScale(diameter , false); | |
| 180 | 10 | const Scalar omega = pow(k, 0.5) / (L*factor); | |
| 181 |
1/2✓ Branch 0 taken 10 times.
✗ Branch 1 not taken.
|
10 | if (print) |
| 182 | { | ||
| 183 | 30 | std::cout << "estimated omega : " << omega << " [1/s]" << std::endl; | |
| 184 | |||
| 185 | } | ||
| 186 | 10 | return omega; | |
| 187 | } | ||
| 188 | |||
| 189 | /*! | ||
| 190 | * \brief Estimates the viscosity tilde based on a formula given in | ||
| 191 | * in the ANSYS Fluent user guide \cite ANSYSUserGuide12 | ||
| 192 | */ | ||
| 193 | 30 | Scalar viscosityTilde(const Scalar velocity, | |
| 194 | const Scalar diameter, | ||
| 195 | const Scalar kinematicViscosity, | ||
| 196 | bool print=verbose) const | ||
| 197 | { | ||
| 198 | using std::abs; | ||
| 199 | using std::sqrt; | ||
| 200 | 30 | const Scalar re_d = reynoldsNumber(velocity, diameter, kinematicViscosity, false); | |
| 201 | 30 | const Scalar viscosityTilde = sqrt(1.5) * abs(velocity) | |
| 202 | 30 | * turbulenceIntensity(re_d) | |
| 203 | 30 | * turbulenceLengthScale(diameter); | |
| 204 |
1/2✓ Branch 0 taken 30 times.
✗ Branch 1 not taken.
|
30 | if (print) |
| 205 | { | ||
| 206 | 90 | std::cout << "estimated nu~ : " << viscosityTilde << " [m^2/s]" << std::endl; | |
| 207 | } | ||
| 208 | 30 | return viscosityTilde; | |
| 209 | } | ||
| 210 | }; | ||
| 211 | } // end namespace Dumux | ||
| 212 | |||
| 213 | #endif | ||
| 214 |