forked from p35462178/gks2d-str
521 lines
14 KiB
C++
521 lines
14 KiB
C++
|
//this libarary includes some basic functions for gks
|
|||
|
|
#include"gks_basic.h"
|
|||
|
|
#define S 110.4
|
|||
|
|
|
|||
|
|
//artifical viscosity for euler problem
|
|||
|
|
double c1_euler = -1.0; //initialization
|
|||
|
|
double c2_euler = -1.0; //initialization
|
|||
|
|
bool is_Prandtl_fix = false; //initialization
|
|||
|
|
double Pr = 1.0; //initialization
|
|||
|
|
|
|||
|
|
|
|||
|
|
|
|||
|
|
// to prepare the basic element for moment calculation
|
|||
|
|
MMDF1d::MMDF1d() { u = 1.0; lambda = 1.0; };
|
|||
|
|
MMDF1d::MMDF1d(double u_in, double lambda_in)
|
|||
|
|
{
|
|||
|
|
u = u_in;
|
|||
|
|
lambda = lambda_in;
|
|||
|
|
calcualte_MMDF1d();
|
|||
|
|
}
|
|||
|
|
void MMDF1d::calcualte_MMDF1d()
|
|||
|
|
{
|
|||
|
|
uwhole[0] = 1;
|
|||
|
|
uwhole[1] = u;
|
|||
|
|
uplus[0] = 0.5 * Alpha(lambda, -u);
|
|||
|
|
uminus[0] = 0.5 * Alpha(lambda, u);
|
|||
|
|
uplus[1] = u * uplus[0] + 0.5 * Beta(lambda, u);
|
|||
|
|
uminus[1] = u * uminus[0] - 0.5 * Beta(lambda, u);
|
|||
|
|
for (int i = 2; i <= 9; i++)
|
|||
|
|
{
|
|||
|
|
uwhole[i] = u * uwhole[i - 1] + 0.5 * (i - 1) / lambda * uwhole[i - 2];
|
|||
|
|
}
|
|||
|
|
for (int i = 2; i <= 9; i++)
|
|||
|
|
{
|
|||
|
|
uplus[i] = u * uplus[i - 1] + 0.5 * (i - 1) / lambda * uplus[i - 2];
|
|||
|
|
uminus[i] = u * uminus[i - 1] + 0.5 * (i - 1) / lambda * uminus[i - 2];
|
|||
|
|
}
|
|||
|
|
xi2 = 0.5 * K / lambda;
|
|||
|
|
xi4 = 0.25 * (K * K + 2 * K) / (lambda * lambda);
|
|||
|
|
//xi6?? how to calculate
|
|||
|
|
xi6 = 0.5 * (K + 4) / lambda * xi4;
|
|||
|
|
|
|||
|
|
for (int i = 0; i < 10; i++)
|
|||
|
|
{
|
|||
|
|
for (int k = 0; k < 4; k++)
|
|||
|
|
{
|
|||
|
|
if ((i + 2 * k) <= 9)
|
|||
|
|
{
|
|||
|
|
if (k == 0)
|
|||
|
|
{
|
|||
|
|
upxi[i][k] = uplus[i];
|
|||
|
|
unxi[i][k] = uminus[i];
|
|||
|
|
}
|
|||
|
|
if (k == 1)
|
|||
|
|
{
|
|||
|
|
upxi[i][k] = uplus[i] * xi2;
|
|||
|
|
unxi[i][k] = uminus[i] * xi2;
|
|||
|
|
}
|
|||
|
|
if (k == 2)
|
|||
|
|
{
|
|||
|
|
upxi[i][k] = uplus[i] * xi4;
|
|||
|
|
unxi[i][k] = uminus[i] * xi4;
|
|||
|
|
}
|
|||
|
|
if (k == 3)
|
|||
|
|
{
|
|||
|
|
upxi[i][k] = uplus[i] * xi6;
|
|||
|
|
unxi[i][k] = uminus[i] * xi6;
|
|||
|
|
}
|
|||
|
|
}
|
|||
|
|
}
|
|||
|
|
}
|
|||
|
|
for (int i = 0; i < 10; i++)
|
|||
|
|
{
|
|||
|
|
for (int k = 0; k < 4; k++)
|
|||
|
|
{
|
|||
|
|
if ((i + 2 * k) <= 9)
|
|||
|
|
{
|
|||
|
|
if (k == 0)
|
|||
|
|
{
|
|||
|
|
uxi[i][k] = uwhole[i];
|
|||
|
|
}
|
|||
|
|
if (k == 1)
|
|||
|
|
{
|
|||
|
|
uxi[i][k] = uwhole[i] * xi2;
|
|||
|
|
}
|
|||
|
|
if (k == 2)
|
|||
|
|
{
|
|||
|
|
uxi[i][k] = uwhole[i] * xi4;
|
|||
|
|
}
|
|||
|
|
if (k == 3)
|
|||
|
|
{
|
|||
|
|
uxi[i][k] = uwhole[i] * xi6;
|
|||
|
|
}
|
|||
|
|
}
|
|||
|
|
}
|
|||
|
|
}
|
|||
|
|
}
|
|||
|
|
|
|||
|
|
//a general G function
|
|||
|
|
void G(int no_u, int no_xi, double* psi, double a[3], MMDF1d m)
|
|||
|
|
{
|
|||
|
|
psi[0] = a[0] * m.uxi[no_u][no_xi] + a[1] * m.uxi[no_u + 1][no_xi] + a[2] * 0.5 * (m.uxi[no_u + 2][no_xi] + m.uxi[no_u][no_xi + 1]);
|
|||
|
|
psi[1] = a[0] * m.uxi[no_u + 1][no_xi] + a[1] * m.uxi[no_u + 2][no_xi] + a[2] * 0.5 * (m.uxi[no_u + 3][no_xi] + m.uxi[no_u + 1][no_xi + 1]);
|
|||
|
|
psi[2] = 0.5 * (a[0] * (m.uxi[no_u + 2][no_xi] + m.uxi[no_u][no_xi + 1]) +
|
|||
|
|
a[1] * (m.uxi[no_u + 3][no_xi] + m.uxi[no_u + 1][no_xi + 1]) +
|
|||
|
|
a[2] * 0.5 * (m.uxi[no_u + 4][no_xi] + m.uxi[no_u][no_xi + 2] + 2 * m.uxi[no_u + 2][no_xi + 1]));
|
|||
|
|
|
|||
|
|
}
|
|||
|
|
void GL(int no_u, int no_xi, double* psi, double a[3], MMDF1d m)
|
|||
|
|
{
|
|||
|
|
psi[0] = a[0] * m.upxi[no_u][no_xi] + a[1] * m.upxi[no_u + 1][no_xi] + a[2] * 0.5 * (m.upxi[no_u + 2][no_xi] + m.upxi[no_u][no_xi + 1]);
|
|||
|
|
psi[1] = a[0] * m.upxi[no_u + 1][no_xi] + a[1] * m.upxi[no_u + 2][no_xi] + a[2] * 0.5 * (m.upxi[no_u + 3][no_xi] + m.upxi[no_u + 1][no_xi + 1]);
|
|||
|
|
psi[2] = 0.5 * (a[0] * (m.upxi[no_u + 2][no_xi] + m.upxi[no_u][no_xi + 1]) +
|
|||
|
|
a[1] * (m.upxi[no_u + 3][no_xi] + m.upxi[no_u + 1][no_xi + 1]) +
|
|||
|
|
a[2] * 0.5 * (m.upxi[no_u + 4][no_xi] + m.upxi[no_u][no_xi + 2] + 2 * m.upxi[no_u + 2][no_xi + 1]));
|
|||
|
|
|
|||
|
|
}
|
|||
|
|
void GR(int no_u, int no_xi, double* psi, double a[3], MMDF1d m)
|
|||
|
|
{
|
|||
|
|
psi[0] = a[0] * m.unxi[no_u][no_xi] + a[1] * m.unxi[no_u + 1][no_xi] + a[2] * 0.5 * (m.unxi[no_u + 2][no_xi] + m.unxi[no_u][no_xi + 1]);
|
|||
|
|
psi[1] = a[0] * m.unxi[no_u + 1][no_xi] + a[1] * m.unxi[no_u + 2][no_xi] + a[2] * 0.5 * (m.unxi[no_u + 3][no_xi] + m.unxi[no_u + 1][no_xi + 1]);
|
|||
|
|
psi[2] = 0.5 * (a[0] * (m.unxi[no_u + 2][no_xi] + m.unxi[no_u][no_xi + 1]) +
|
|||
|
|
a[1] * (m.unxi[no_u + 3][no_xi] + m.unxi[no_u + 1][no_xi + 1]) +
|
|||
|
|
a[2] * 0.5 * (m.unxi[no_u + 4][no_xi] + m.unxi[no_u][no_xi + 2] + 2 * m.unxi[no_u + 2][no_xi + 1]));
|
|||
|
|
|
|||
|
|
}
|
|||
|
|
|
|||
|
|
//moments of the maxwellian distribution function
|
|||
|
|
MMDF::MMDF() { u = 1.0; v = 1.0; lambda = 1.0; };
|
|||
|
|
MMDF::MMDF(double u_in, double v_in, double lambda_in)
|
|||
|
|
{
|
|||
|
|
u = u_in;
|
|||
|
|
v = v_in;
|
|||
|
|
lambda = lambda_in;
|
|||
|
|
calcualte_MMDF();
|
|||
|
|
}
|
|||
|
|
void MMDF::calcualte_MMDF()
|
|||
|
|
{
|
|||
|
|
//the moment of u in whole domain
|
|||
|
|
uwhole[0] = 1;
|
|||
|
|
uwhole[1] = u;
|
|||
|
|
//the moment of u in half domain (>0 && <0)
|
|||
|
|
uplus[0] = 0.5 * Alpha(lambda, -u);
|
|||
|
|
uminus[0] = 1.0 - uplus[0];
|
|||
|
|
uplus[1] = u * uplus[0] + 0.5 * Beta(lambda, u);
|
|||
|
|
uminus[1] = u - uplus[1];
|
|||
|
|
double overlambda = 1.0 / lambda;
|
|||
|
|
for (int i = 2; i <= 6; i++)
|
|||
|
|
{
|
|||
|
|
uwhole[i] = u * uwhole[i - 1] + 0.5 * (i - 1) * overlambda * uwhole[i - 2];
|
|||
|
|
}
|
|||
|
|
for (int i = 2; i <= 6; i++)
|
|||
|
|
{
|
|||
|
|
uplus[i] = u * uplus[i - 1] + 0.5 * (i - 1) * overlambda * uplus[i - 2];
|
|||
|
|
uminus[i] = uwhole[i] - uplus[i];
|
|||
|
|
}
|
|||
|
|
//the moment of v in whole domain
|
|||
|
|
//the moment of v in half domain is useless
|
|||
|
|
vwhole[0] = 1;
|
|||
|
|
vwhole[1] = v;
|
|||
|
|
for (int i = 2; i <= 6; i++)
|
|||
|
|
{
|
|||
|
|
vwhole[i] = v * vwhole[i - 1] + 0.5 * (i - 1) * overlambda * vwhole[i - 2];
|
|||
|
|
}
|
|||
|
|
//the moment of xi (kesi)
|
|||
|
|
xi2 = 0.5 * K * overlambda;
|
|||
|
|
xi4 = 0.25 * (K * K + 2 * K) * overlambda * overlambda;
|
|||
|
|
//xi6?? how to calculate
|
|||
|
|
//xi6 = 0.5*(K + 4) / lambda*xi4; // no use
|
|||
|
|
|
|||
|
|
//get other variables by the above mement variables
|
|||
|
|
//i, j, k, means i power of u, j power of v, and 2k power of xi
|
|||
|
|
for (int i = 0; i < 7; i++)
|
|||
|
|
{
|
|||
|
|
for (int j = 0; j < 7; j++)
|
|||
|
|
{
|
|||
|
|
for (int k = 0; k < 3; k++)
|
|||
|
|
{
|
|||
|
|
if ((i + j + 2 * k) <= 6)
|
|||
|
|
{
|
|||
|
|
if (k == 0)
|
|||
|
|
{
|
|||
|
|
// u, half domain > 0; v, whole domain; xi, equals 0
|
|||
|
|
upvxi[i][j][k] = uplus[i] * vwhole[j];
|
|||
|
|
// u, half domain < 0; v, whole domain; xi, equals 0
|
|||
|
|
unvxi[i][j][k] = uminus[i] * vwhole[j];
|
|||
|
|
// u, whole domain ; v, whole domain; xi, equals 0
|
|||
|
|
uvxi[i][j][k] = upvxi[i][j][k] + unvxi[i][j][k];
|
|||
|
|
// = (uplus[i] + uminus[i]) * vwhole[j] = uwhole[i] * vwhole[j]
|
|||
|
|
}
|
|||
|
|
if (k == 1)
|
|||
|
|
{
|
|||
|
|
// u, half domain > 0; v, whole domain; xi, equals 2
|
|||
|
|
upvxi[i][j][k] = uplus[i] * vwhole[j] * xi2;
|
|||
|
|
// u, half domain < 0; v, whole domain; xi, equals 2
|
|||
|
|
unvxi[i][j][k] = uminus[i] * vwhole[j] * xi2;
|
|||
|
|
// u, whole domain ; v, whole domain; xi, equals 2
|
|||
|
|
uvxi[i][j][k] = upvxi[i][j][k] + unvxi[i][j][k];
|
|||
|
|
}
|
|||
|
|
if (k == 2)
|
|||
|
|
{
|
|||
|
|
// u, half domain > 0; v, whole domain; xi, equals 4
|
|||
|
|
upvxi[i][j][k] = uplus[i] * vwhole[j] * xi4;
|
|||
|
|
// u, half domain < 0; v, whole domain; xi, equals 4
|
|||
|
|
unvxi[i][j][k] = uminus[i] * vwhole[j] * xi4;
|
|||
|
|
// u, whole domain ; v, whole domain; xi, equals 4
|
|||
|
|
uvxi[i][j][k] = upvxi[i][j][k] + unvxi[i][j][k];
|
|||
|
|
}
|
|||
|
|
|
|||
|
|
}
|
|||
|
|
}
|
|||
|
|
}
|
|||
|
|
}
|
|||
|
|
}
|
|||
|
|
|
|||
|
|
MMDF1st::MMDF1st() { u = 1.0; v = 1.0; lambda = 1.0; };
|
|||
|
|
MMDF1st::MMDF1st(double u_in, double v_in, double lambda_in)
|
|||
|
|
{
|
|||
|
|
u = u_in;
|
|||
|
|
v = v_in;
|
|||
|
|
lambda = lambda_in;
|
|||
|
|
calcualte_MMDF1st();
|
|||
|
|
}
|
|||
|
|
void MMDF1st::calcualte_MMDF1st()
|
|||
|
|
{
|
|||
|
|
uwhole[0] = 1.0;
|
|||
|
|
uwhole[1] = u;
|
|||
|
|
|
|||
|
|
//return erfc(sqrt(lambda)*u);
|
|||
|
|
|
|||
|
|
//return exp(-lambda*u*u) / sqrt(pi*lambda);
|
|||
|
|
uminus[0] = 0.5 * Alpha(lambda, u);
|
|||
|
|
uplus[0] = 1.0 - uminus[0];
|
|||
|
|
double beta = Beta(lambda, u);
|
|||
|
|
//uplus[1] = u*uplus[0] + 0.5*Beta(lambda, u);
|
|||
|
|
//uminus[1] = u*uminus[0] - 0.5*Beta(lambda, u);
|
|||
|
|
uplus[1] = u * uplus[0] + 0.5 * beta;
|
|||
|
|
uminus[1] = u * uminus[0] - 0.5 * beta;
|
|||
|
|
double overlambda = 1.0 / lambda;
|
|||
|
|
for (int i = 2; i <= 3; i++)
|
|||
|
|
{
|
|||
|
|
uwhole[i] = u * uwhole[i - 1] + 0.5 * (i - 1) * overlambda * uwhole[i - 2];
|
|||
|
|
}
|
|||
|
|
for (int i = 2; i <= 3; i++)
|
|||
|
|
{
|
|||
|
|
uplus[i] = u * uplus[i - 1] + 0.5 * (i - 1) * overlambda * uplus[i - 2];
|
|||
|
|
uminus[i] = uwhole[i] - uplus[i];
|
|||
|
|
}
|
|||
|
|
vwhole[0] = 1;
|
|||
|
|
vwhole[1] = v;
|
|||
|
|
for (int i = 2; i <= 2; i++)
|
|||
|
|
{
|
|||
|
|
vwhole[i] = v * vwhole[i - 1] + 0.5 * (i - 1) * overlambda * vwhole[i - 2];
|
|||
|
|
}
|
|||
|
|
xi2 = 0.5 * K * overlambda;
|
|||
|
|
}
|
|||
|
|
|
|||
|
|
|
|||
|
|
|
|||
|
|
void Collision(double* w0, double left, double right, MMDF& m2, MMDF& m3)
|
|||
|
|
{
|
|||
|
|
// get the equilibrium variables by collision
|
|||
|
|
w0[0] = left * m2.uplus[0] + right * m3.uminus[0];
|
|||
|
|
w0[1] = left * m2.uplus[1] + right * m3.uminus[1];
|
|||
|
|
w0[2] = left * m2.vwhole[1] * m2.uplus[0] + right * m3.vwhole[1] * m3.uminus[0];
|
|||
|
|
w0[3] = 0.5 * left * (m2.uplus[2] + m2.uplus[0] * m2.vwhole[2] + m2.uplus[0] * m2.xi2) +
|
|||
|
|
0.5 * right * (m3.uminus[2] + m3.uminus[0] * m3.vwhole[2] + m3.uminus[0] * m3.xi2);
|
|||
|
|
}
|
|||
|
|
void Collision(double* w0, double left, double right, MMDF1st& m2, MMDF1st& m3)
|
|||
|
|
{
|
|||
|
|
w0[0] = left * m2.uplus[0] + right * m3.uminus[0];
|
|||
|
|
w0[1] = left * m2.uplus[1] + right * m3.uminus[1];
|
|||
|
|
w0[2] = left * m2.vwhole[1] * m2.uplus[0] + right * m3.vwhole[1] * m3.uminus[0];
|
|||
|
|
w0[3] = 0.5 * left * (m2.uplus[2] + m2.uplus[0] * m2.vwhole[2] + m2.uplus[0] * m2.xi2) +
|
|||
|
|
0.5 * right * (m3.uminus[2] + m3.uminus[0] * m3.vwhole[2] + m3.uminus[0] * m3.xi2);
|
|||
|
|
}
|
|||
|
|
|
|||
|
|
|
|||
|
|
|
|||
|
|
void A(double* a, double der[4], double prim[4])
|
|||
|
|
{
|
|||
|
|
//Note: solve the a coefficient (slope information), equivalent to solving the matrix to get a
|
|||
|
|
//input: der[4], the gradient of conservative variables, density, momentumX, momentumY, energy, partialW/patialX
|
|||
|
|
//input: prim[4], the primary variables (density, u, v, lambda) at one point or one interface
|
|||
|
|
//output: a, solve the results
|
|||
|
|
|
|||
|
|
double R4, R3, R2;
|
|||
|
|
double overden = 1.0 / prim[0];
|
|||
|
|
//double overlambda = 1.0 / prim[3];
|
|||
|
|
R4 = der[3] * overden - 0.5 * (prim[1] * prim[1] + prim[2] * prim[2] + 0.5 * (K + 2) / prim[3]) * der[0] * overden;
|
|||
|
|
R3 = (der[2] - prim[2] * der[0]) * overden;
|
|||
|
|
R2 = (der[1] - prim[1] * der[0]) * overden;
|
|||
|
|
a[3] = (4.0 / (K + 2)) * prim[3] * prim[3] * (2 * R4 - 2 * prim[1] * R2 - 2 * prim[2] * R3);
|
|||
|
|
a[2] = 2 * prim[3] * R3 - prim[2] * a[3];
|
|||
|
|
a[1] = 2 * prim[3] * R2 - prim[1] * a[3];
|
|||
|
|
a[0] = der[0] * overden - prim[1] * a[1] - prim[2] * a[2] - 0.5 * a[3] * (prim[1] * prim[1] + prim[2] * prim[2] + 0.5 * (K + 2) / prim[3]);
|
|||
|
|
|
|||
|
|
//R4 = der[3] / prim[0] - 0.5*(prim[1] * prim[1] + prim[2] * prim[2] + 0.5*(K + 2) / prim[3])*der[0] / prim[0];
|
|||
|
|
//R3 = (der[2] - prim[2]*der[0]) / prim[0];
|
|||
|
|
//R2 = (der[1] - prim[1]*der[0]) / prim[0];
|
|||
|
|
//a[3] = 4.0 * prim[3]*prim[3] / (K + 2)*(2.0 * R4 - 2.0 * prim[1]*R2 - 2.0 * prim[2]*R3);
|
|||
|
|
//a[2] = 2.0 * prim[3]*R3 - prim[2]*a[3];
|
|||
|
|
//a[1] = 2.0 * prim[3]*R2 - prim[1]*a[3];
|
|||
|
|
//a[0] = der[0] / prim[0] - prim[1]*a[1] - prim[2]*a[2] - 0.5*a[3] * (prim[1]*prim[1] + prim[2]*prim[2] + 0.5*(K + 2) / prim[3]);
|
|||
|
|
|
|||
|
|
}
|
|||
|
|
|
|||
|
|
|
|||
|
|
|
|||
|
|
double Get_Tau_NS(double density0, double lambda0)
|
|||
|
|
{
|
|||
|
|
if (tau_type == Euler)
|
|||
|
|
{
|
|||
|
|
return 0.0;
|
|||
|
|
}
|
|||
|
|
else if (tau_type == NS)
|
|||
|
|
{
|
|||
|
|
if (Mu > 0.0)
|
|||
|
|
{
|
|||
|
|
//cout << "here" << endl;
|
|||
|
|
return 2.0 * Mu * lambda0 / density0;
|
|||
|
|
}
|
|||
|
|
else if (Nu > 0.0)
|
|||
|
|
{
|
|||
|
|
return 2 * Nu * lambda0;
|
|||
|
|
}
|
|||
|
|
else
|
|||
|
|
{
|
|||
|
|
return 0.0;
|
|||
|
|
}
|
|||
|
|
}
|
|||
|
|
else if (tau_type == Sutherland)
|
|||
|
|
{
|
|||
|
|
return TauNS_Sutherland(density0, lambda0);
|
|||
|
|
}
|
|||
|
|
else
|
|||
|
|
{
|
|||
|
|
return TauNS_power_law(density0, lambda0);
|
|||
|
|
}
|
|||
|
|
|
|||
|
|
}
|
|||
|
|
double Get_Tau(double density_left, double density_right, double density0, double lambda_left, double lambda_right, double lambda0, double dt)
|
|||
|
|
{
|
|||
|
|
if (tau_type == Euler)
|
|||
|
|
{
|
|||
|
|
if (c1_euler <= 0 && c2_euler <= 0)
|
|||
|
|
{
|
|||
|
|
return 0.0;
|
|||
|
|
}
|
|||
|
|
else
|
|||
|
|
{
|
|||
|
|
double C = c2_euler * abs(density_left / lambda_left - density_right / lambda_right) / abs(density_left / lambda_left + density_right / lambda_right);
|
|||
|
|
//C+= c2_euler * abs(density_left - density_right) / abs(density_left + density_right);
|
|||
|
|
//if (C < 10)
|
|||
|
|
//{
|
|||
|
|
return c1_euler * dt + dt * C;
|
|||
|
|
//}
|
|||
|
|
//else
|
|||
|
|
//{
|
|||
|
|
// return c1_euler*dt + 10*dt;
|
|||
|
|
// //return c1_euler*dt + dt*C;
|
|||
|
|
//}
|
|||
|
|
}
|
|||
|
|
}
|
|||
|
|
else if (tau_type == NS)
|
|||
|
|
{
|
|||
|
|
double tau_n = c2_euler * abs(density_left / lambda_left - density_right / lambda_right) / abs(density_left / lambda_left + density_right / lambda_right) * dt;
|
|||
|
|
if (tau_n != tau_n)
|
|||
|
|
{
|
|||
|
|
tau_n = 0.0;
|
|||
|
|
}
|
|||
|
|
if ((Mu > 0.0 && Nu > 0.0) || (Mu < 0.0 && Nu < 0.0))
|
|||
|
|
{
|
|||
|
|
return 0.0;
|
|||
|
|
}
|
|||
|
|
else
|
|||
|
|
{
|
|||
|
|
if (Mu > 0.0)
|
|||
|
|
{
|
|||
|
|
return tau_n + 2.0 * Mu * lambda0 / density0;
|
|||
|
|
}
|
|||
|
|
else if (Nu > 0.0)
|
|||
|
|
{
|
|||
|
|
return tau_n + 2 * Nu * lambda0;
|
|||
|
|
}
|
|||
|
|
else
|
|||
|
|
{
|
|||
|
|
return 0.0;
|
|||
|
|
}
|
|||
|
|
}
|
|||
|
|
}
|
|||
|
|
else if (tau_type == Sutherland)
|
|||
|
|
{
|
|||
|
|
double tau_n = c2_euler * abs(density_left / lambda_left - density_right / lambda_right) / abs(density_left / lambda_left + density_right / lambda_right) * dt;
|
|||
|
|
if (tau_n != tau_n)
|
|||
|
|
{
|
|||
|
|
tau_n = 0.0;
|
|||
|
|
}
|
|||
|
|
if ((Mu > 0.0 && Nu > 0.0) || (Mu < 0.0 && Nu < 0.0))
|
|||
|
|
{
|
|||
|
|
|
|||
|
|
return 0.0;
|
|||
|
|
}
|
|||
|
|
else
|
|||
|
|
{
|
|||
|
|
if (Mu > 0.0)
|
|||
|
|
{
|
|||
|
|
return tau_n + TauNS_Sutherland(density0, lambda0);
|
|||
|
|
|
|||
|
|
}
|
|||
|
|
else if (Nu > 0.0)
|
|||
|
|
{
|
|||
|
|
return tau_n + 2 * Nu * lambda0;
|
|||
|
|
}
|
|||
|
|
else
|
|||
|
|
{
|
|||
|
|
return 0.0;
|
|||
|
|
}
|
|||
|
|
}
|
|||
|
|
}
|
|||
|
|
else if (tau_type == Power_law)
|
|||
|
|
{
|
|||
|
|
|
|||
|
|
double tau_n = c2_euler * abs(density_left / lambda_left - density_right / lambda_right) / abs(density_left / lambda_left + density_right / lambda_right) * dt;
|
|||
|
|
if (tau_n != tau_n)
|
|||
|
|
{
|
|||
|
|
tau_n = 0.0;
|
|||
|
|
}
|
|||
|
|
if ((Mu > 0.0 && Nu > 0.0) || (Mu < 0.0 && Nu < 0.0))
|
|||
|
|
{
|
|||
|
|
return 0.0;
|
|||
|
|
}
|
|||
|
|
else
|
|||
|
|
{
|
|||
|
|
if (Mu > 0.0)
|
|||
|
|
{
|
|||
|
|
return tau_n + TauNS_power_law(density0, lambda0);
|
|||
|
|
}
|
|||
|
|
else if (Nu > 0.0)
|
|||
|
|
{
|
|||
|
|
return tau_n + 2 * Nu * lambda0;
|
|||
|
|
}
|
|||
|
|
else
|
|||
|
|
{
|
|||
|
|
return 0.0;
|
|||
|
|
}
|
|||
|
|
}
|
|||
|
|
|
|||
|
|
}
|
|||
|
|
else
|
|||
|
|
{
|
|||
|
|
return 0.0;
|
|||
|
|
}
|
|||
|
|
}
|
|||
|
|
double TauNS_Sutherland(double density0, double lambda0)
|
|||
|
|
{
|
|||
|
|
double T0 = 1.0 / 2.0 / (lambda0 * R_gas);
|
|||
|
|
double mu = Mu * pow(T0 / T_inf, 1.5) * (T_inf + S) / (T0 + S);
|
|||
|
|
return 2.0 * mu * lambda0 / density0;
|
|||
|
|
}
|
|||
|
|
double TauNS_power_law(double density0, double lambda0)
|
|||
|
|
{
|
|||
|
|
double T0 = 1.0 / 2.0 / (lambda0 * R_gas);
|
|||
|
|
double mu = Mu * pow(T0 / T_inf, 0.76);
|
|||
|
|
return 2.0 * mu * lambda0 / density0;
|
|||
|
|
}
|
|||
|
|
|
|||
|
|
|
|||
|
|
|
|||
|
|
//计算λ
|
|||
|
|
double Lambda(double density, double u, double densityE)
|
|||
|
|
{
|
|||
|
|
return (K + 1.0) * 0.25 * (density / (densityE - 0.5 * density * (u * u)));
|
|||
|
|
}
|
|||
|
|
double Lambda(double density, double u, double v, double densityE)
|
|||
|
|
{
|
|||
|
|
return (K + 2.0) * 0.25 * (density / (densityE - 0.5 * density * (u * u + v * v)));
|
|||
|
|
}
|
|||
|
|
|
|||
|
|
void Convar_to_ULambda_1d(double* primvar, double convar[3])
|
|||
|
|
{
|
|||
|
|
primvar[0] = convar[0];
|
|||
|
|
primvar[1] = U(convar[0], convar[1]);
|
|||
|
|
primvar[2] = Lambda(convar[0], primvar[1], convar[2]);
|
|||
|
|
}
|
|||
|
|
|
|||
|
|
//计算α，α=erfc（±√λ*U）
|
|||
|
|
double Alpha(double lambda, double u)
|
|||
|
|
{
|
|||
|
|
return erfc(sqrt(lambda) * u);
|
|||
|
|
//return erfcmy(sqrt(lambda)*u);
|
|||
|
|
}
|
|||
|
|
//计算β，β=e^(-λ*U^2)/(√πλ)
|
|||
|
|
double Beta(double lambda, double u)
|
|||
|
|
{
|
|||
|
|
double pi = 3.14159265358979323846;
|
|||
|
|
return exp(-lambda * u * u) / sqrt(pi * lambda);
|
|||
|
|
}
|
|||
|
|
|
|||
|
|
|
|||
|
|
//solution of matrix equation b=Ma 1D
|
|||
|
|
void Microslope(double* a, double der[3], double prim[3])
|
|||
|
|
{
|
|||
|
|
double R4, R2;
|
|||
|
|
R4 = der[2] / prim[0] - 0.5 * (prim[1] * prim[1] + 0.5 * (K + 1) / prim[2]) * der[0] / prim[0];
|
|||
|
|
R2 = (der[1] - prim[1] * der[0]) / prim[0];
|
|||
|
|
a[2] = 4 * prim[2] * prim[2] / (K + 1) * (2 * R4 - 2 * prim[1] * R2);
|
|||
|
|
a[1] = 2 * prim[2] * R2 - prim[1] * a[2];
|
|||
|
|
a[0] = der[0] / prim[0] - prim[1] * a[1] - 0.5 * a[2] * (prim[1] * prim[1] + 0.5 * (K + 1) / prim[2]);
|
|||
|
|
}
|
|||
|
|
|
|||
|
|
void Convar_to_ULambda_2d(double* primvar, double convar[4]) //density, U, V, lambda
|
|||
|
|
{
|
|||
|
|
primvar[0] = convar[0];
|
|||
|
|
primvar[1] = U(convar[0], convar[1]);
|
|||
|
|
primvar[2] = V(convar[0], convar[2]);
|
|||
|
|
//here primvar[3] refers lambda
|
|||
|
|
primvar[3] = Lambda(convar[0], primvar[1], primvar[2], convar[3]);
|
|||
|
|
}
|
|||
|
|
|