📄 rjfuncs.cpp
字号:
// RJFuncs.cpp: implementation of the RJFuncs class.
//
//////////////////////////////////////////////////////////////////////
#include "stdio.h"
#include "math.h"
#include "RJFuncs.h"
//////////////////////////////////////////////////////////////////////
// Construction/Destruction
//////////////////////////////////////////////////////////////////////
RJFuncs::RJFuncs()
{
}
RJFuncs::~RJFuncs()
{
}
double RJFuncs::PI(double LM,double K)
//PI函数，表示流体静压与总压的关系，
//是速度系数LM和比热比K的函数
{
double PI_O;
PI_O=pow(1.0-(K-1.0)*LM*LM/(K+1.0),K/(K-1.0));
return (PI_O);
}
double RJFuncs::TOR(double LM,double K)
//TOR函数，表达静温与总温的关系，
//是速度系数LM与比热比K的函数
{
double TOR_O;
TOR_O=1.0-(K-1.0)*LM*LM/(K+1.0);
return (TOR_O);
}
double RJFuncs::APS(double LM,double K)
//APS函数，表达与的关系，
//是速度系数LM与比热比K的函数
{
double APS_O;
APS_O=pow(1.0-(K-1.0)*LM*LM/(K+1.0),1.0/(K-1.0));
return (APS_O);
}
double RJFuncs::QL( double LM,double K)
{
double QL_O;
QL_O=LM*pow((K+1.0)/2.0,1.0/(K-1.0))*APS(LM,K);
return (QL_O);
}
double RJFuncs::YL(double LM,double K)
{
double YL_O;
YL_O=LM*pow((K+1.0)/2.0,1.0/(K-1))/(1.0-(K-1.0)*LM*LM/(K+1.0));
return (YL_O);
}
double RJFuncs::FL(double LM,double K)
{
double FL_O;
FL_O=(1+LM*LM)*pow(1.0-(K-1.0)*LM*LM/(K+1.0),1.0/(K-1.0));
return (FL_O);
}
double RJFuncs::RL(double LM,double K)
{
double RL_O;
RL_O=(1.0-(K-1.0)*LM*LM/(K+1.0))/(1+LM*LM);
return (RL_O);
}
double RJFuncs::ZL(double LM)
{
double ZL_O;
ZL_O=LM+1.0/LM;
return (ZL_O);
}
double RJFuncs::Ma_To_LM(double Ma,double K)
{
double LM,Tmp_V;
Tmp_V=Ma*Ma*(K+1.0)/(2.0+(K-1)*Ma*Ma);
LM=sqrt(Tmp_V);
return (LM);
}
double RJFuncs::QL_To_LM(double QL_O,double K,int N)
{
double LM,A,B,C;
if ((1-QL_O)<0.001)
{
LM=1.0;
return(LM);
}
else
{
if (N==0)
{A=1.0;B=0.0;}
else {A=sqrt((K+1.0)/(K-1.0))-0.001;B=1.0;}
do
{
C=B-(QL(B,K)-QL_O)*(B-A)/(QL(B,K)-QL(A,K));
A=B;
B=C;
}
while(fabs(QL(C,K)-QL_O)>0.0001);
LM=C;
return(LM);
}
}
double RJFuncs::BZ(double R,double K)
{
double BZ_O;
BZ_O=sqrt(2.0*R*(K+1.0)/K);
return (BZ_O);
}
double RJFuncs::QM(double R,double K)
{
double QM_O;
QM_O=sqrt(pow(2.0/(K+1.0),(K+1.0)/(K-1.0))*K/R);
return (QM_O);
}
double RJFuncs::PI_To_LM(double PI_O,double K)
{
double LM;
LM=sqrt((1-pow(PI_O,(K-1.0)/K))*(K+1.0)/(K-1.0));
return (LM);
}
void RJFuncs::Air(double HE,double &TH,double &PH,double &RU,double &AV,double &GH)
{
double HG,TS,PS,RUS,GHS;//T,P,Ru,g of heignt 0;
double W; //coefficient
HG=HE/(1.0+HE/6356766.0);
TS=288.15;
PS=101325.0;
RUS=1.225;
GHS=9.80665;
if(HG<11019.1)
{
W=1.0-HG/44330.8;
TH=TS*W;
PH=PS*pow(W,5.2559);
RU=RUS*pow(W,4.25559);
}
else
{
if(HG<20063.1)
{
W=exp((14964.7-HG)/6341.6);
TH=216.65;
PH=PS*0.11953*W;
RU=RUS*0.15898*W;
}
else
{
W=1.0+(HG-24902.1)/221552.0;
TH=221.552*W;
PH=PS*0.025158*pow(W,-34.1629);
RU=RUS*0.032722*pow(W,-35.1629);
}
}
AV=20.0468*sqrt(TH);
GH=GHS/pow(1.0+HG/6356766,2.0);
}
double RJFuncs::CAZ1(double x[],double y[],int n,double t)
{ int i,j,k,m;
double z,s;
z=0.0;
if (n<1) return(z);
if (n==1) { z=y[0]; return(z);}
if (n==2)
{ z=(y[0]*(t-x[1])-y[1]*(t-x[0]))/(x[0]-x[1]);
return(z);
}
if (t<=x[1]) { k=0; m=2;}
else if (t>=x[n-2]) { k=n-3; m=n-1;}
else
{ k=1; m=n;
while (m-k!=1)
{ i=(k+m)/2;
if (t<x[i-1]) m=i;
else k=i;
}
k=k-1; m=m-1;
if (fabs(t-x[k])<fabs(t-x[m])) k=k-1;
else m=m+1;
}
z=0.0;
for (i=k;i<=m;i++)
{ s=1.0;
for (j=k;j<=m;j++)
if (j!=i) s=s*(t-x[j])/(x[i]-x[j]);
z=z+s*y[i];
}
return(z);
}
double RJFuncs::CAZ2(double x[],double y[],double z[],int n,int m,double u,double v)
{ int nn,mm,ip,iq,i,j,k,l;
double b[3],h,w;
nn=3;
if (n<=3) { ip=0; nn=n;}
else if (u<=x[1]) ip=0;
else if (u>=x[n-2]) ip=n-3;
else
{ i=1; j=n;
while (((i-j)!=1)&&((i-j)!=-1))
{ l=(i+j)/2;
if (u<x[l-1]) j=l;
else i=l;
}
if (fabs(u-x[i-1])<fabs(u-x[j-1])) ip=i-2;
else ip=i-1;
}
mm=3;
if (m<=3) { iq=0; mm=m;}
else if (v<=y[1]) iq=0;
else if (v>=y[m-2]) iq=m-3;
else
{ i=1; j=m;
while (((i-j)!=1)&&((i-j)!=-1))
{ l=(i+j)/2;
if (v<y[l-1]) j=l;
else i=l;
}
if (fabs(v-y[i-1])<fabs(v-y[j-1])) iq=i-2;
else iq=i-1;
}
for (i=0;i<=nn-1;i++)
{ b[i]=0.0;
for (j=0;j<=mm-1;j++)
{ k=m*(ip+i)+(iq+j);
h=z[k];
for (k=0;k<=mm-1;k++)
if (k!=j)
h=h*(v-y[iq+k])/(y[iq+j]-y[iq+k]);
b[i]=b[i]+h;
}
}
w=0.0;
for (i=0;i<=nn-1;i++)
{ h=b[i];
for (j=0;j<=nn-1;j++)
if (j!=i)
h=h*(u-x[ip+j])/(x[ip+i]-x[ip+j]);
w=w+h;
}
return(w);
}
double RJFuncs::AQL2( double QL_O,double K)
{ //求QL(LM)之大于1.0之LM
double a=1.0;
double b=sqrt((K+1.0)/(K-1.0))-0.001;
double f1,f2,f3;
double LM;
if(fabs(QL_O-1.0)<0.001){LM=1.0;}
else
{
do
{
f1=QL_O-QL(a,K);
f2=QL_O-QL(b,K);
f3=QL_O-QL((a+b)/2.0,K);
if (f1*f3<0) b=(a+b)/2.0;
else a=(a+b)/2.0;
if (fabs(f3)<1e-10){LM=(a+b)/2.0;break;}
}while(1);
}
return (LM);
}
double RJFuncs::CAZ22(double x[],double y[],double z[],int n,int m,double u,double v)
{
int n1,n2,m1,m2,i;
double z11,z12,z21,z22;
if(x[1]>x[0])
{
if(u<=x[0])
{
n1=0;n2=1;
}
else
{
if(u>=x[n-1])
{
n1=n-2;
n2=n-1;
}
else
{
for(i=1;i<n;i++)
{
if(u<x[i]) {n1=i-1;n2=i;break;}
}
}
}
}
else
{
if(u>=x[0])
{
n1=0;n2=1;
}
else
{
if(u<=x[n-1])
{
n1=n-2;
n2=n-1;
}
else
{
for(i=1;i<n;i++)
{
if(u>x[i]) {n1=i-1;n2=i;break;}
}
}
}
}
if(y[1]>y[0])
{
if(v<=y[0])
{
m1=0;m2=1;
}
else
{
if(v>=y[m-1])
{
m1=m-2;
m2=m-1;
}
else
{
for(i=1;i<m;i++)
{
if(v<y[i]) {m1=i-1;m2=i;break;}
}
}
}
}
else
{
if(v>=y[0])
{
m1=0;m2=1;
}
else
{
if(v<=y[m-1])
{
m1=m-2;
m2=m-1;
}
else
{
for(i=1;i<m;i++)
{
if(v>y[i]) {m1=i-1;m2=i;break;}
}
}
}
}
z11=z[n1*m+m1];
z12=z[n1*m+m2];
z21=z[n2*m+m1];
z22=z[n2*m+m2];
////////////////
// Znm
// n1 n2
// m1 z11,z21 z1
// m2 z12,z22 z2
//////////////
double z1,z2,z3;
z1=z11+(z21-z11)*(u-x[n1])/(x[n2]-x[n1]);
z2=z12+(z22-z12)*(u-x[n1])/(x[n2]-x[n1]);
z3=z1+(z2-z1)*(v-y[m1])/(y[m2]-y[m1]);
return(z3);
}
double RJFuncs::CAZ11(double x[],double y[],int n,double t)
{
int n1,n2,i;
double z;
if(x[1]>x[0])
{
if(t<=x[0])
{
n1=0;n2=1;
}
else
{
if(t>=x[n-1])
{
n1=n-2;
n2=n-1;
}
else
{
for(i=1;i<n;i++)
{
if(t<x[i]) {n1=i-1;n2=i;break;}
}
}
}
}
else
{
if(t>=x[0])
{
n1=0;n2=1;
}
else
{
if(t<=x[n-1])
{
n1=n-2;
n2=n-1;
}
else
{
for(i=1;i<n;i++)
{
if(t>x[i]) {n1=i-1;n2=i;break;}
}
}
}
}
z=y[n1]+(y[n2]-y[n1])*(t-x[n1])/(x[n2]-x[n1]);
return (z);
}
double RJFuncs::LM_To_Ma(double LM,double K)
{
double M;
M=sqrt(2.0/(K+1.0))*LM*pow(1.0-(K-1.0)*LM*LM/(K+1.0),-0.5);
return (M);
}
double RJFuncs::ShockAngleOfMDA(double Ma,double k)
//某马赫数最大气流偏转角对应激波角
{
double sin_em;
sin_em=1.0/(k*Ma*Ma)*((k+1.0)*Ma*Ma/4.0-1.0+sqrt((k+1.0)*(1+(k-1.0)*Ma*Ma/2.0+(k+1)*pow(Ma,4)/16.0)));
double ShockAngle=asin(sqrt(sin_em))*RAD;
return (ShockAngle);
}
double RJFuncs::MaxDeflectAngle(double Ma,double k)
//对应某马赫数的最大气流偏转角
{
double ShockAngle=ShockAngleOfMDA(Ma,k)/RAD;
double DeflectAngle=RAD*atan(1.0/(((k+1)*Ma*Ma/(2*Ma*Ma*pow(sin(ShockAngle),2)-2.0)-1.0)*tan(ShockAngle)));
return (DeflectAngle);
}
double RJFuncs::ShockAngleFunc(double Ma,double DeflectAngle,double ShockAngle,double k)
//激波角与马赫数，气流偏转角的函数关系
{
double erro;
erro=tan(DeflectAngle/RAD)-(pow(Ma,2.0)*pow(sin(ShockAngle/RAD),2.0)-1.0)/((pow(Ma,2.0)*((k+1.0)/2.0-pow(sin(ShockAngle/RAD),2.0))+1.0)*tan(ShockAngle/RAD));
return (erro);
}
double RJFuncs::WeakShockAngle(double Ma,double DeflectAngle,double k)
{
double s_angle1,s_angle2,s_angle;
double erro1,erro2,erro;
//二分法求激波角
s_angle1=1e-6;
s_angle2=ShockAngleOfMDA(Ma,k);
Begin:
erro1=ShockAngleFunc(Ma,DeflectAngle,s_angle1,k);
erro2=ShockAngleFunc(Ma,DeflectAngle,s_angle2,k);
s_angle=(s_angle1+s_angle2)/2.0;
erro=ShockAngleFunc(Ma,DeflectAngle,s_angle,k);
if(fabs(erro)<1e-10) return(s_angle);
else
{
if(erro*erro1<0) s_angle2=s_angle;
else s_angle1=s_angle;
goto Begin;
}
}
double RJFuncs::StrongShockAngle(double Ma,double DeflectAngle,double k)
{
double s_angle1,s_angle2,s_angle;
double erro1,erro2,erro;
//二分法求激波角
s_angle1=ShockAngleOfMDA(Ma,k);
s_angle2=90;
for(;;)
{
erro1=ShockAngleFunc(Ma,DeflectAngle,s_angle1,k);
erro2=ShockAngleFunc(Ma,DeflectAngle,s_angle2,k);
s_angle=(s_angle1+s_angle2)/2.0;
erro=ShockAngleFunc(Ma,DeflectAngle,s_angle,k);
if(fabs(erro)<1e-10) break;
else
{
if(erro*erro1<0) s_angle2=s_angle;
else s_angle1=s_angle;
}
}
return(s_angle);
}
double RJFuncs::MachAfterShock(double Ma,double ShockAngle,double k)
//激波后马赫数
{
double msin;
msin=pow(Ma,2)*pow(sin(ShockAngle/RAD),2);
double Ma2;
Ma2=sqrt((Ma*Ma+2.0/(k-1))/(2*k*msin/(k-1.0)-1.0)+(Ma*Ma-msin)/((k-1)*msin*0.5+1.0));
return (Ma2);
}
double RJFuncs::PresRatioAfterShock(double Ma,double ShockAngle,double k)
//激波后压强关系
{
double P2_P1;
P2_P1=2.0*k*Ma*Ma*sin(ShockAngle/RAD)*sin(ShockAngle/RAD)/(k+1.0)-(k-1.0)/(k+1.0);
return(P2_P1);
}
double RJFuncs::DensRatioAfterShock(double Ma,double ShockAngle,double k)
//激波后密度关系
{
double msin;
msin=pow(Ma,2)*pow(sin(ShockAngle/RAD),2);
return((k+1.0)*msin/(2+(k-1.0)*msin));
}
double RJFuncs::PresRecoverCoef(double Ma,double ShockAngle,double k)
//激波后总压恢复系数
{
double RatioPres=PresRatioAfterShock(Ma,ShockAngle,k);
double RatioDens=DensRatioAfterShock(Ma,ShockAngle,k);
return(pow(RatioDens,k/(k-1.0))*pow(1.0/RatioPres,1.0/(k-1.0)));
}
⌨️ 快捷键说明
复制代码
Ctrl + C
搜索代码
Ctrl + F
全屏模式
F11
切换主题
Ctrl + Shift + D
显示快捷键
?
增大字号
Ctrl + =
减小字号
Ctrl + -