From 4ac7fa02c5a416217d0f626a4dd806f82d555ff8 Mon Sep 17 00:00:00 2001
From: acrovato <a.crovato@uliege.be>
Date: Mon, 12 Sep 2022 13:51:11 +0200
Subject: [PATCH] =?UTF-8?q?Change=20speed=20of=20sound=20clamp=20from=20Pa?=
=?UTF-8?q?d=C3=A9=20to=20constant.?=
MIME-Version: 1.0
Content-Type: text/plain; charset=UTF-8
Content-Transfer-Encoding: 8bit
---
dart/src/dart.h | 2 +-
dart/src/wF0El.cpp | 80 +++++++++++++++++++++++-----------------------
dart/src/wF0El.h | 31 ++++++++++--------
3 files changed, 58 insertions(+), 55 deletions(-)
diff --git a/dart/src/dart.h b/dart/src/dart.h
index ad5628c..4770c32 100644
--- a/dart/src/dart.h
+++ b/dart/src/dart.h
@@ -69,7 +69,7 @@ class Adjoint;
class F0Ps;
class F0El;
class F0ElC;
-class F0ElRhoBase;
+class F0ElSpeedSound;
class F0ElRho;
class F0ElRhoClamp;
class F0ElRhoLin;
diff --git a/dart/src/wF0El.cpp b/dart/src/wF0El.cpp
index bea2834..4cfe6cf 100644
--- a/dart/src/wF0El.cpp
+++ b/dart/src/wF0El.cpp
@@ -46,40 +46,40 @@ double F0ElC::evalGrad(Element const &e, std::vector<double> const &phi, size_t
* The input _mC2 is the square of the critical Mach number (used to compute
* the critical velocity).
*/
-F0ElRhoBase::F0ElRhoBase(double _mInf, double _mC2) : F0El(), gamma(1.4), mInf(_mInf)
+F0ElSpeedSound::F0ElSpeedSound(double _mInf, double _mC2) : F0El(), gamma(1.4), mInf(_mInf)
{
uC = sqrt(_mC2 * (1 / (mInf * mInf) + (gamma - 1) / 2) / (1 + (gamma - 1) / 2 * _mC2)); // critical velocity
- rC = 1 + 0.5 * (gamma - 1) * mInf * mInf * (1 - uC * uC); // r(uC)
- beta = ((gamma - 1) * mInf * mInf * uC) / rC; // dr/du(uC) / -r(uC)
+ aC = 1 + 0.5 * (gamma - 1) * mInf * mInf * (1 - uC * uC); // a^2(uC)
+ beta = 0 * ((gamma - 1) * mInf * mInf * uC) / aC; // da^2/du(uC) / -a^2(uC)
}
/**
* @brief Evaluate the base of the isentropic density
*
- * r = [1 + (gamma-1)/2*Mi^2*(1-u^2)], u <= uC
- * r = r(uC) / (1-(dr/du(uC)/f(uC))*(u/uC-1)), u > uC
+ * a^2 = [1 + (gamma-1)/2*Mi^2*(1-u^2)], u <= uC
+ * a^2 = a^2(uC) / (1-(da^2/du(uC)/a^2(uC))*(u/uC-1)), u > uC
*/
-double F0ElRhoBase::eval(Element const &e, std::vector<double> const &phi, size_t k) const
+double F0ElSpeedSound::eval(Element const &e, std::vector<double> const &phi, size_t k) const
{
double u = e.computeGradient(phi, k).norm(); // norm of velocity
if (u <= uC)
return 1 + 0.5 * (gamma - 1) * mInf * mInf * (1 - u * u);
else
- return rC / (1 + beta * (u / uC - 1));
+ return aC / (1 + beta * (u / uC - 1));
}
/**
* @brief Evaluate the gradient of the base of the isentropic density
*
- * dr = -(gamma-1)*Mi^2 * u, u <= uC
- * dr = -r(uC) / (1-(dr/du(uC)/f(uC))*(u/uC-1))^2 * (dr/du(uC)/f(uC)) / uC, u > uC
+ * da^2 = -(gamma-1)*Mi^2 * u, u <= uC
+ * da^2 = -a^2(uC) / (1-(da^2/du(uC)/a^2(uC))*(u/uC-1))^2 * (da^2/du(uC)/a^2(uC)) / uC, u > uC
* Computed value has to be multiplied by u to recover actual gradient
*/
-double F0ElRhoBase::evalGrad(Element const &e, std::vector<double> const &phi, size_t k) const
+double F0ElSpeedSound::evalGrad(Element const &e, std::vector<double> const &phi, size_t k) const
{
double u = e.computeGradient(phi, k).norm(); // norm of velocity
if (u <= uC)
return -(gamma - 1) * mInf * mInf;
else
- return -rC / ((1 + beta * (u / uC - 1)) * (1 + beta * (u / uC - 1))) * beta / uC / u;
+ return -aC / ((1 + beta * (u / uC - 1)) * (1 + beta * (u / uC - 1))) * beta / uC / u;
}
/**
@@ -90,8 +90,8 @@ double F0ElRhoBase::evalGrad(Element const &e, std::vector<double> const &phi, s
double F0ElRho::eval(Element const &e, std::vector<double> const &phi, size_t k) const
{
double u = e.computeGradient(phi, k).norm(); // norm of velocity
- double r = 1 + 0.5 * (gamma - 1) * mInf * mInf * (1 - u * u);
- return sqrt(r * r * r * r * r); // particularized to gamma=1.4
+ double a2 = 1 + 0.5 * (gamma - 1) * mInf * mInf * (1 - u * u);
+ return sqrt(a2 * a2 * a2 * a2 * a2); // particularized to gamma=1.4
}
/**
* @brief Evaluate the isentropic density gradient
@@ -102,31 +102,31 @@ double F0ElRho::eval(Element const &e, std::vector<double> const &phi, size_t k)
double F0ElRho::evalGrad(Element const &e, std::vector<double> const &phi, size_t k) const
{
double u = e.computeGradient(phi, k).norm(); // norm of velocity
- double r = 1 + 0.5 * (gamma - 1) * mInf * mInf * (1 - u * u);
- return -mInf * mInf * sqrt(r * r * r); // particularized to gamma=1.4
+ double a2 = 1 + 0.5 * (gamma - 1) * mInf * mInf * (1 - u * u);
+ return -mInf * mInf * sqrt(a2 * a2 * a2); // particularized to gamma=1.4
}
/**
* @brief Evaluate the (clamped) isentropic density
*
- * rho = r^(1/(gamma-1))
+ * rho = a^(2/(gamma-1))
*/
double F0ElRhoClamp::eval(Element const &e, std::vector<double> const &phi, size_t k) const
{
- double r = F0ElRhoBase::eval(e, phi, k);
- return sqrt(r * r * r * r * r); // particularized to gamma=1.4
+ double a2 = F0ElSpeedSound::eval(e, phi, k);
+ return sqrt(a2 * a2 * a2 * a2 * a2); // particularized to gamma=1.4
}
/**
* @brief Evaluate the (clamped) isentropic density gradient
*
- * drho = 1/(gamma-1) * r^(2-gamma) * dr
+ * drho = 1/(gamma-1) * a^(2*(2-gamma)) * da^2
* Computed value has to be multiplied by u to recover actual gradient
*/
double F0ElRhoClamp::evalGrad(Element const &e, std::vector<double> const &phi, size_t k) const
{
- double r = F0ElRhoBase::eval(e, phi, k);
- double dr = F0ElRhoBase::evalGrad(e, phi, k);
- return 1 / (gamma - 1) * sqrt(r * r * r) * dr; // particularized to gamma=1.4
+ double a2 = F0ElSpeedSound::eval(e, phi, k);
+ double da2 = F0ElSpeedSound::evalGrad(e, phi, k);
+ return 1 / (gamma - 1) * sqrt(a2 * a2 * a2) * da2; // particularized to gamma=1.4
}
/**
@@ -171,26 +171,26 @@ double F0ElMach::evalGrad(Element const &e, std::vector<double> const &phi, size
/**
* @brief Evaluate the (clamped) isentropic Mach number
*
- * M = u / sqrt(r/Mi^2)
+ * M = u / sqrt(a^2/Mi^2)
*/
double F0ElMachClamp::eval(Element const &e, std::vector<double> const &phi, size_t k) const
{
double u = e.computeGradient(phi, k).norm(); // norm of velocity
- double r = F0ElRhoBase::eval(e, phi, k);
- return mInf * u / sqrt(r);
+ double a2 = F0ElSpeedSound::eval(e, phi, k);
+ return mInf * u / sqrt(a2);
}
/**
* @brief Evaluate the (clamped) isentropic Mach number gradient
*
- * dM = 1 / sqrt(r/Mi^2) + u / sqrt(r^3/Mi^2) * dr
+ * dM = 1 / sqrt(a^2/Mi^2) + u / sqrt(a^6/Mi^2) * da^2
* Computed value has to be multiplied by u to recover actual gradient
*/
double F0ElMachClamp::evalGrad(Element const &e, std::vector<double> const &phi, size_t k) const
{
double u = e.computeGradient(phi, k).norm(); // norm of velocity
- double r = F0ElRhoBase::eval(e, phi, k);
- double dr = F0ElRhoBase::evalGrad(e, phi, k);
- return mInf * (1 / sqrt(r) / u - 0.5 * u / sqrt(r * r * r) * dr);
+ double a2 = F0ElSpeedSound::eval(e, phi, k);
+ double da2 = F0ElSpeedSound::evalGrad(e, phi, k);
+ return mInf * (1 / sqrt(a2) / u - 0.5 * u / sqrt(a2 * a2 * a2) * da2);
}
/**
@@ -216,8 +216,8 @@ double F0ElMachLin::evalGrad(Element const &e, std::vector<double> const &phi, s
double F0ElCp::eval(Element const &e, std::vector<double> const &phi, size_t k) const
{
double u2 = e.computeGradient(phi, k).squaredNorm(); // velocity squared
- double r = 1 + 0.5 * (gamma - 1) * mInf * mInf * (1 - u2);
- return 2 / (gamma * mInf * mInf) * (sqrt(r * r * r * r * r * r * r) - 1); // particularized to gamma=1.4
+ double a2 = 1 + 0.5 * (gamma - 1) * mInf * mInf * (1 - u2);
+ return 2 / (gamma * mInf * mInf) * (sqrt(a2 * a2 * a2 * a2 * a2 * a2 * a2) - 1); // particularized to gamma=1.4
}
/**
* @brief Evaluate the isentropic pressure coefficient gradient
@@ -228,31 +228,31 @@ double F0ElCp::eval(Element const &e, std::vector<double> const &phi, size_t k)
double F0ElCp::evalGrad(Element const &e, std::vector<double> const &phi, size_t k) const
{
double u2 = e.computeGradient(phi, k).squaredNorm(); // velocity squared
- double r = 1 + 0.5 * (gamma - 1) * mInf * mInf * (1 - u2);
- return -2 * sqrt(r * r * r * r * r); // particularized to gamma=1.4
+ double a2 = 1 + 0.5 * (gamma - 1) * mInf * mInf * (1 - u2);
+ return -2 * sqrt(a2 * a2 * a2 * a2 * a2); // particularized to gamma=1.4
}
/**
* @brief Evaluate the (clamped) isentropic pressure coefficient
*
- * Cp = 2/(gamma*Mi^2) * (r^(gamma/(gamma-1)) - 1)
+ * Cp = 2/(gamma*Mi^2) * (a^(2*gamma/(gamma-1)) - 1)
*/
double F0ElCpClamp::eval(Element const &e, std::vector<double> const &phi, size_t k) const
{
- double r = F0ElRhoBase::eval(e, phi, k);
- return 2 / (gamma * mInf * mInf) * (sqrt(r * r * r * r * r * r * r) - 1); // particularized to gamma=1.4
+ double a2 = F0ElSpeedSound::eval(e, phi, k);
+ return 2 / (gamma * mInf * mInf) * (sqrt(a2 * a2 * a2 * a2 * a2 * a2 * a2) - 1); // particularized to gamma=1.4
}
/**
* @brief Evaluate the (clamped) isentropic pressure coefficient gradient
*
- * dCp = 2/((gamma-1)*Mi^2) * r^(1/(gamma-1)) * dr
+ * dCp = 2/((gamma-1)*Mi^2) * a^(2/(gamma-1)) * da^2
* Computed value has to multiplied by u to recover actual gradient
*/
double F0ElCpClamp::evalGrad(Element const &e, std::vector<double> const &phi, size_t k) const
{
- double r = F0ElRhoBase::eval(e, phi, k);
- double dr = F0ElRhoBase::evalGrad(e, phi, k);
- return 2 / ((gamma - 1) * mInf * mInf) * sqrt(r * r * r * r * r) * dr; // particularized to gamma=1.4
+ double a2 = F0ElSpeedSound::eval(e, phi, k);
+ double da2 = F0ElSpeedSound::evalGrad(e, phi, k);
+ return 2 / ((gamma - 1) * mInf * mInf) * sqrt(a2 * a2 * a2 * a2 * a2) * da2; // particularized to gamma=1.4
}
/**
diff --git a/dart/src/wF0El.h b/dart/src/wF0El.h
index 4f4a774..e2ff691 100644
--- a/dart/src/wF0El.h
+++ b/dart/src/wF0El.h
@@ -58,25 +58,26 @@ public:
};
/**
- * @brief Base of isentropic density
+ * @brief Speed of sound squared
*
- * This function computes the base of the isentropic density, and is then used
+ * This function computes the square of the speed of sound, and is then used
* to compute the actual density, Mach number and pressure coefficient. This
* function is clamped for large velocities using a Padé approximation and
- * particularized to diatomic gas.
+ * particularized to diatomic gas. The function is normalized by its
+ * freestream value.
*/
-class DART_API F0ElRhoBase : public F0El
+class DART_API F0ElSpeedSound : public F0El
{
protected:
double gamma; ///< heat capacity ratio (diatomic gas only)
double mInf; ///< freestream Mach number
private:
double uC; ///< critical velocity
- double rC; ///< critical density base
+ double aC; ///< critical speed of sound squared
double beta; ///< ratio for Padé
public:
- F0ElRhoBase(double _mInf, double _mC2 = 5);
+ F0ElSpeedSound(double _mInf, double _mC2 = 5);
virtual double eval(tbox::Element const &e, std::vector<double> const &phi,
size_t k) const override;
@@ -87,7 +88,8 @@ public:
/**
* @brief Isentropic density
*
- * The density is particularized to diatomic gas.
+ * The density is particularized to diatomic gas. The density is normalized by
+ * its freestream value.
*/
class DART_API F0ElRho : public F0El
{
@@ -107,12 +109,13 @@ public:
* @brief Isentropic clamped density
*
* The density is clamped (Padé approximation) for large values of velocity and
- * particularized to diatomic gas.
+ * particularized to diatomic gas. The density is normalized by its
+ * freestream value.
*/
-class DART_API F0ElRhoClamp : public F0ElRhoBase
+class DART_API F0ElRhoClamp : public F0ElSpeedSound
{
public:
- F0ElRhoClamp(double _mInf, double _mC2 = 5) : F0ElRhoBase(_mInf, _mC2) {}
+ F0ElRhoClamp(double _mInf, double _mC2 = 5) : F0ElSpeedSound(_mInf, _mC2) {}
virtual double eval(tbox::Element const &e, std::vector<double> const &phi,
size_t k) const override;
@@ -159,10 +162,10 @@ public:
* The Mach number is clamped (Padé approximation) for large values of velocity
* and particularized to diatomic gas.
*/
-class DART_API F0ElMachClamp : public F0ElRhoBase
+class DART_API F0ElMachClamp : public F0ElSpeedSound
{
public:
- F0ElMachClamp(double _mInf, double _mC2 = 5) : F0ElRhoBase(_mInf, _mC2) {}
+ F0ElMachClamp(double _mInf, double _mC2 = 5) : F0ElSpeedSound(_mInf, _mC2) {}
virtual double eval(tbox::Element const &e, std::vector<double> const &phi,
size_t k) const override;
@@ -209,10 +212,10 @@ public:
* The pressure coefficient is clamped (Padé approximation) for large values of
* velocity and particularized to diatomic gas.
*/
-class DART_API F0ElCpClamp : public F0ElRhoBase
+class DART_API F0ElCpClamp : public F0ElSpeedSound
{
public:
- F0ElCpClamp(double _mInf, double _mC2 = 5) : F0ElRhoBase(_mInf, _mC2) {}
+ F0ElCpClamp(double _mInf, double _mC2 = 5) : F0ElSpeedSound(_mInf, _mC2) {}
virtual double eval(tbox::Element const &e, std::vector<double> const &phi,
size_t k) const override;
--
GitLab