-
Paul Dechamps authored
Ensure all numbers are doubles in Closures, Fluxes and Driver & debug info update in BoundaryLayer & ensure nIter in 3D rae test
Paul Dechamps authoredEnsure all numbers are doubles in Closures, Fluxes and Driver & debug info update in BoundaryLayer & ensure nIter in 3D rae test
DFluxes.cpp 8.12 KiB
#include "DFluxes.h"
#include "DBoundaryLayer.h"
#include <Eigen/Dense>
using namespace blast;
using namespace Eigen;
/**
* @brief Compute residual vector of the integral boundary layer equations.
*
* @param iPoint Marker id.
* @param bl Boundary layer region.
* @return VectorXd
*/
VectorXd Fluxes::computeResiduals(size_t iPoint, BoundaryLayer *bl)
{
size_t nVar = bl->getnVar();
std::vector<double> u(nVar, 0.);
for (size_t k = 0; k < nVar; ++k)
u[k] = bl->u[iPoint * nVar + k];
return blLaws(iPoint, bl, u);
}
/**
* @brief Compute Jacobian of the integral boundary layer equations.
*
* @param iPoint Marker id.
* @param bl Boundary layer region.
* @param sysRes Residual of the IBL at point iPoint.
* @param eps Pertubation from current solution used to compute the Jacobian.
* @return MatrixXd
*/
MatrixXd Fluxes::computeJacobian(size_t iPoint, BoundaryLayer *bl, VectorXd const &sysRes, double eps)
{
size_t nVar = bl->getnVar();
MatrixXd JacMatrix = MatrixXd::Zero(5, 5);
std::vector<double> uUp(nVar, 0.);
for (size_t k = 0; k < nVar; ++k)
uUp[k] = bl->u[iPoint * nVar + k];
double varSave = 0.;
for (size_t iVar = 0; iVar < nVar; ++iVar)
{
varSave = uUp[iVar];
uUp[iVar] += eps;
JacMatrix.col(iVar) = (blLaws(iPoint, bl, uUp) - sysRes) / eps;
uUp[iVar] = varSave;
}
return JacMatrix;
}
/**
* @brief Integral boundary layer equation.
*
* @param iPoint Marker id.
* @param bl Boundary layer region.
* @param u Solution vector at point iPoint.
* @return VectorXd
*/
VectorXd Fluxes::blLaws(size_t iPoint, BoundaryLayer *bl, std::vector<double> u)
{
size_t nVar = bl->getnVar();
MatrixXd timeMatrix = MatrixXd::Zero(5, 5);
VectorXd spaceVector = VectorXd::Zero(5);
// Dissipation ratio.
double dissipRatio;
if (bl->name == "wake")
dissipRatio = 0.9;
else
dissipRatio = 1.;
bl->Ret[iPoint] = std::max(u[3] * u[0] * Re, 1.0); // Reynolds number based on theta.
bl->deltaStar[iPoint] = u[0] * u[1]; // Displacement thickness.
// Boundary layer closure relations.
if (bl->regime[iPoint] == 0)
{
// Laminar closure relations.
std::vector<double> lamParam(8, 0.);
bl->closSolver->laminarClosures(lamParam, u[0], u[1], bl->Ret[iPoint], bl->blEdge->getMe(iPoint), bl->name);
bl->Hstar[iPoint] = lamParam[0];
bl->Hstar2[iPoint] = lamParam[1];
bl->Hk[iPoint] = lamParam[2];
bl->cf[iPoint] = lamParam[3];
bl->cd[iPoint] = lamParam[4];
bl->delta[iPoint] = lamParam[5];
bl->ctEq[iPoint] = lamParam[6];
bl->us[iPoint] = lamParam[7];
}
else if (bl->regime[iPoint] == 1)
{
// Turbulent closure relations.
std::vector<double> turbParam(8, 0.);
bl->closSolver->turbulentClosures(turbParam, u[0], u[1], u[4], bl->Ret[iPoint], bl->blEdge->getMe(iPoint), bl->name);
bl->Hstar[iPoint] = turbParam[0];
bl->Hstar2[iPoint] = turbParam[1];
bl->Hk[iPoint] = turbParam[2];
bl->cf[iPoint] = turbParam[3];
bl->cd[iPoint] = turbParam[4];
bl->delta[iPoint] = turbParam[5];
bl->ctEq[iPoint] = turbParam[6];
bl->us[iPoint] = turbParam[7];
}
else
{
std::cout << "Wrong regime\n"
<< std::endl;
throw;
}
// Gradients computation.
double dx = (bl->mesh->getLoc(iPoint) - bl->mesh->getLoc(iPoint - 1));
double dxExt = (bl->mesh->getExt(iPoint) - bl->mesh->getExt(iPoint - 1));
double dt_dx = (u[0] - bl->u[(iPoint - 1) * nVar + 0]) / dx;
double dH_dx = (u[1] - bl->u[(iPoint - 1) * nVar + 1]) / dx;
double due_dx = (u[3] - bl->u[(iPoint - 1) * nVar + 3]) / dx;
double dct_dx = (u[4] - bl->u[(iPoint - 1) * nVar + 4]) / dx;
double dHstar_dx = (bl->Hstar[iPoint] - bl->Hstar[iPoint - 1]) / dx;
double dueExt_dx = (bl->blEdge->getVt(iPoint) - bl->blEdge->getVt(iPoint - 1)) / dxExt;
double ddeltaStarExt_dx = (bl->blEdge->getDeltaStar(iPoint) - bl->blEdge->getDeltaStar(iPoint - 1)) / dxExt;
double c = 4 / (M_PI * dx);
double cExt = 4 / (M_PI * dxExt);
// Amplification routine.
double dN_dx = 0.0;
double ax = 0.0;
if (bl->xtrF == -1.0 && bl->regime[iPoint] == 0)
{
// No forced transition and laminar regime.
ax = amplificationRoutine(bl->Hk[iPoint], bl->Ret[iPoint], u[0]);
dN_dx = (u[2] - bl->u[(iPoint - 1) * nVar + 2]) / dx;
}
double Mea = bl->blEdge->getMe(iPoint);
double Hstara = bl->Hstar[iPoint];
double Hstar2a = bl->Hstar2[iPoint];
double Hka = bl->Hk[iPoint];
double cfa = bl->cf[iPoint];
double cda = bl->cd[iPoint];
double deltaa = bl->delta[iPoint];
double ctEqa = bl->ctEq[iPoint];
double usa = bl->us[iPoint];
// Space part.
spaceVector(0) = dt_dx + (2. + u[1] - Mea * Mea) * (u[0] / u[3]) * due_dx - cfa / 2.;
spaceVector(1) = u[0] * dHstar_dx + (2 * Hstar2a + Hstara * (1 - u[1])) * u[0] / u[3] * due_dx - 2. * cda + Hstara * cfa / 2.;
spaceVector(2) = dN_dx - ax;
spaceVector(3) = due_dx - c * (u[1] * dt_dx + u[0] * dH_dx) - dueExt_dx + cExt * ddeltaStarExt_dx;
if (bl->regime[iPoint] == 1)
spaceVector(4) = (2. * deltaa / u[4]) * dct_dx - 5.6 * (ctEqa - u[4] * dissipRatio) - 2. * deltaa * (4. / (3. * u[1] * u[0]) * (cfa / 2. - (((Hka - 1.) / (6.7 * Hka * dissipRatio)) * ((Hka - 1.) / (6.7 * Hka * dissipRatio)))) - 1. / u[3] * due_dx);
// Time part.
timeMatrix(0, 0) = u[1] / u[3];
timeMatrix(0, 1) = u[0] / u[3];
timeMatrix(0, 3) = u[0] * u[1] / (u[3] * u[3]);
timeMatrix(1, 0) = (1. + u[1] * (1. - Hstara)) / u[3];
timeMatrix(1, 1) = (1. - Hstara) * u[0] / u[3];
timeMatrix(1, 3) = (2. - Hstara * u[1]) * u[0] / (u[3] * u[3]);
timeMatrix(2, 2) = 1.;
timeMatrix(3, 0) = -c * u[1];
timeMatrix(3, 1) = -c * u[0];
timeMatrix(3, 3) = 1.;
if (bl->regime[iPoint] == 1)
{
timeMatrix(4, 3) = 2. * deltaa / (usa * u[3] * u[3]);
timeMatrix(4, 4) = 2. * deltaa / (u[3] * u[4] * usa);
}
else
timeMatrix(4, 4) = 1.0;
return timeMatrix.partialPivLu().solve(spaceVector);
}
/**
* @brief Amplification routines of the e^N method for transition capturing.
*
* @param Hk Kinematic shape parameter.
* @param Ret Reynolds number based on the momentum thickness.
* @param theta Momentum thickness.
* @return double
*/
double Fluxes::amplificationRoutine(double Hk, double Ret, double theta) const
{
double dgr = 0.08;
double Hk1 = 3.5;
double Hk2 = 4.;
double Hmi = (1. / (Hk - 1.));
double logRet_crit = 2.492 * std::pow(1 / (Hk - 1.), 0.43) + 0.7 * (std::tanh(14 * Hmi - 9.24) + 1.0); // value at which the amplification starts to grow
double ax = 0.;
if (Ret <= 0.)
Ret = 1.;
if (std::log10(Ret) < logRet_crit - dgr)
ax = 0.;
else
{
double rNorm = (std::log10(Ret) - (logRet_crit - dgr)) / (2. * dgr);
double Rfac = 0.;
if (rNorm <= 0)
Rfac = 0.;
if (rNorm >= 1)
Rfac = 1.;
else
Rfac = 3. * rNorm * rNorm - 2 * rNorm * rNorm * rNorm;
// Normal envelope amplification rate
double m = -0.05 + 2.7 * Hmi - 5.5 * Hmi * Hmi + 3. * Hmi * Hmi * Hmi + 0.1 * std::exp(-20. * Hmi);
double arg = 3.87 * Hmi - 2.52;
double dndRet = 0.028 * (Hk - 1.) - 0.0345 * std::exp(-(arg * arg));
ax = (m * dndRet / theta) * Rfac;
// Set correction for separated profiles
if (Hk > Hk1)
{
double Hnorm = (Hk - Hk1) / (Hk2 - Hk1);
double Hfac = 0.;
if (Hnorm >= 1.)
Hfac = 1.;
else
Hfac = 3. * Hnorm * Hnorm - 2. * Hnorm * Hnorm * Hnorm;
double ax1 = ax;
double gro = 0.3 + 0.35 * std::exp(-0.15 * (Hk - 5.));
double tnr = std::tanh(1.2 * (std::log10(Ret) - gro));
double ax2 = (0.086 * tnr - 0.25 / std::pow((Hk - 1.), 1.5)) / theta;
if (ax2 < 0.)
ax2 = 0.;
ax = Hfac * ax2 + (1. - Hfac) * ax1;
if (ax < 0.)
ax = 0.;
}
}
return ax;
}