Skip to content
Snippets Groups Projects
Commit 4ac7fa02 authored by Adrien Crovato's avatar Adrien Crovato
Browse files

Change speed of sound clamp from Padé to constant.

parent 23d4022a
No related branches found
No related tags found
No related merge requests found
Pipeline #8677 passed
Hide whitespace changes
Inline Side-by-side
dart/src/dart.h
+ 1
− 1
View file @ 4ac7fa02
......@@ -69,7 +69,7 @@ class Adjoint;
class F0Ps;
class F0El;
class F0ElC;
class F0ElRhoBase;
class F0ElSpeedSound;
class F0ElRho;
class F0ElRhoClamp;
class F0ElRhoLin;
......
dart/src/wF0El.cpp
+ 40
− 40
View file @ 4ac7fa02
......@@ -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
}
/**
......
dart/src/wF0El.h
+ 17
− 14
View file @ 4ac7fa02
......@@ -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;
......
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment