📄 main.cpp
字号:
{
nx /= nr;
ny /= nr;
nz /= nr;
}
else
{
nx /= -nr;
ny /= -nr;
nz /= -nr;
}
normal_x[num] = nx;
normal_y[num] = ny;
normal_z[num] = nz;
// 计算面元面积
a1 = sqrt((x2-x1)*(x2-x1)+(y2-y1)*(y2-y1)+(z2-z1)*(z2-z1));
b1 = sqrt((x3-x2)*(x3-x2)+(y3-y2)*(y3-y2)+(z3-z2)*(z3-z2));
a2 = sqrt((x4-x1)*(x4-x1)+(y4-y1)*(y4-y1)+(z4-z1)*(z4-z1));
b2 = sqrt((x3-x4)*(x3-x4)+(y3-y4)*(y3-y4)+(z3-z4)*(z3-z4));
c = sqrt((x3-x1)*(x3-x1)+(y3-y1)*(y3-y1)+(z3-z1)*(z3-z1));
S1 = 0.5 * (a1 + b1 + c);
S2 = 0.5 * (a2 + b2 + c);
A_face[num] = sqrt(S1*(S1-a1)*(S1-b1)*(S1-c)) + sqrt(S2*(S2-a2)*(S2-b2)*(S2-c));
//计算基于当地来流速度的各面元切向
u1 = u + q * r0z - r * r0y;
v1 = v + r * r0x - b * r0z;
w1 = w + b * r0y - q * r0x;
tx = normal_y[num]*(normal_x[num]*v1-normal_y[num]*u1)-normal_z[num]*(normal_z[num]*u1-normal_x[num]*w1);
ty = normal_z[num]*(normal_y[num]*w1-normal_z[num]*v1)-normal_x[num]*(normal_x[num]*v1-normal_y[num]*u1);
tz = normal_x[num]*(normal_z[num]*u1-normal_x[num]*w1)-normal_y[num]*(normal_y[num]*w1-normal_z[num]*v1);
tr = sqrt(tx*tx + ty*ty + tz*tz);
if (fabs(tr)<1e-6)
{
printf("计算有误");
return -1;
}
if (tx < 0)
{
tx /= -tr;
ty /= -tr;
tz /= -tr;
}
else
{
tx /= tr;
ty /= tr;
tz /= tr;
}
parallel_x[num] = tx;
parallel_y[num] = ty;
parallel_z[num] = tz;
//计算当地来流和面元夹角正弦、余弦值
//---------------------------------------------------------------------------
test1 = u1 * normal_x[num] + v1 * normal_y[num] + w1 * normal_z[num];
test2 = u1 * parallel_x[num] + v1 * parallel_y[num] + w1 * parallel_z[num];
V1mo = sqrt(u1*u1 + v1*v1 + w1*w1); //当地来流速度
Ma_face[num] = V1mo/sqrt(gama * R * *T); //当地马赫数
S_face[num] = V1mo/sqrt(2 * R * *T); //当地速度比
cos_face[num] = test2/V1mo;
if (test1>0)
{
sin_face[num] = test1/V1mo;
}
else if (test1<0)
{
sin_face[num] = test1/V1mo;
}
else if (test1 == 0)
{
if(fabs(u1+normal_x[num])>fabs(u1))
{
sin_face[num] = -1;
}
else
{
sin_face[num] = 1;
}
}
//-----------------------------------------------------------------------------------------------------------------------
dF_x[num] = 0;
dF_y[num] = 0;
dF_z[num] = 0;
dM_x[num] = 0;
dM_y[num] = 0;
dM_z[num] = 0;
if (Kn<=0.001)
{
if (Ma_face[num]>=7)//牛顿公式
{
if(test1>0)
{
cp_face[num] = 0;
}
else if(test1<0)
{
cp_face[num] = 2*(sin_face[num])*(sin_face[num]);
}
else if (test1==0)
{
if(fabs(u1+normal_x[num])>fabs(u1))
{
cp_face[num] = 0;
}
else
{
cp_face[num] = 2;
}
}
dF_x[num] = (0.5* *rho *(V1mo * V1mo) * cp_face[num]+*p)*A_face[num] * normal_x[num];
dF_y[num] = (0.5* *rho *(V1mo * V1mo) * cp_face[num]+*p) *A_face[num] * normal_y[num];
dF_z[num] = (0.5* *rho *(V1mo * V1mo) * cp_face[num]+*p) *A_face[num] * normal_z[num];
F[0] = F[0] + dF_x[num];
F[1] = F[1] + dF_y[num];
F[2] = F[2] + dF_z[num];
dM_x[num] = (r0z * dF_y[num] - r0y * dF_z[num]);
dM_y[num] = (r0z * dF_x[num] - r0x * dF_z[num]);
dM_z[num] = (r0x * dF_y[num] - r0y * dF_x[num]);
M[0] =M[0] + dM_x[num];
M[1] =M[1] + dM_y[num];
M[2] =M[2] + dM_z[num];
Cd = (F[0]*u + F[1]*v + F[2]*w)/(0.5 * *rho * vmo * vmo * ref_A * L);
CMx = M[0] / (0.5 * *rho * vmo * vmo * ref_A * L);
CMy = M[1] / (0.5 * *rho * vmo * vmo * ref_A * L);
CMz = M[2] / (0.5 * *rho * vmo * vmo * ref_A * L);
}
else //修正的牛顿公式
{
if(test1>0)
{
cp_face[num] = 0;
}
if(test1<0)
{
cp_face[num] = 2 / gama / pow(Ma_face[num],2) * (pow(pow(gama+1,2)*pow(Ma_face[num],2)/(4*gama*pow(Ma_face[num],2)-2*gama+2), gama/(gama-1))*(1-gama+2*gama*pow(Ma_face[num],2))/(gama+1)-1) * pow(sin_face[num],2);
}
if (test1==0)
{
if(fabs(u1+normal_x[num])>fabs(u1))
{
cp_face[num] = 0;
}
else
{
cp_face[num] = 2 / gama / pow(Ma_face[num],2) * (pow(pow(gama+1,2)*pow(Ma_face[num],2)/(4*gama*pow(Ma_face[num],2)-2*gama+2), gama/(gama-1))*(1-gama+2*gama*pow(Ma_face[num],2))/(gama+1)-1);
}
}
F[0] = F[0]-1 * (0.5* *rho *(V1mo * V1mo) * cp_face[num]+*p)*A_face[num] * normal_x[num];
F[1] = F[1]-1 * (0.5* *rho *(V1mo * V1mo) * cp_face[num]+*p)*A_face[num] * normal_y[num];
F[2] = F[2]-1 * (0.5* *rho *(V1mo * V1mo) * cp_face[num]+*p)*A_face[num] * normal_z[num];
M[0] =M[0] + (r0z * F[1]-r0y * F[2]);
M[1] =M[1] + (r0z * F[0]-r0x * F[2]);
M[2] =M[2] + (r0x * F[1]-r0y * F[0]);
Cd = F[0] / (0.5 * *rho * vmo * vmo * ref_A);
Cl = F[1] / (0.5 * *rho * vmo * vmo * ref_A);
Cz = F[2] / (0.5 * *rho * vmo * vmo * ref_A);
CMx = M[0] / (0.5 * *rho * vmo * vmo * ref_A * L);
CMy = M[1] / (0.5 * *rho * vmo * vmo * ref_A * L);
CMz = M[2] / (0.5 * *rho * vmo * vmo * ref_A * L);
}
}
//--------------------------------------------------------------------------------------------------------------------------------
//自由分子流
else if (Kn>100)
{
p_face[num] = *rho * (V1mo * V1mo) / 2 / pow(S_face[num],2) * (((2-xigema1)/sqrt(PI)*S_face[num]*sin_face[num]+xigema1/2*sqrt(Tw / *T))*exp(-pow(S_face[num]*sin_face[num],2)) + ((2-xigema1)*(0.5+pow(S_face[num]*sin_face[num],2))+xigema1/2*sqrt(PI*Tw / *T)*S_face[num]*sin_face[num])*(1+erf(S_face[num]*sin_face[num])));
tao_face[num] = -xigema * *rho * (V1mo * V1mo) * cos_face[num] / 2 / sqrt(PI) / S_face[num] * (exp(-pow(S_face[num]*sin_face[num],2)) + sqrt(PI)*S_face[num]*sin_face[num]*(1+erf(S_face[num]*sin_face[num])));
dF_x[num] = (-p_face[num] * normal_x[num] - tao_face[num] * parallel_x[num])*A_face[num];
dF_y[num] = (-p_face[num] * normal_y[num] - tao_face[num] * parallel_y[num])*A_face[num];
dF_z[num] = (-p_face[num] * normal_z[num] - tao_face[num] * parallel_z[num])*A_face[num];
F[0] = F[0] + dF_x[num];
F[1] = F[1] + dF_y[num];
F[2] = F[2] + dF_z[num];
dM_x[num] = (r0z * dF_y[num] - r0y * dF_z[num]);
dM_y[num] = (r0z * dF_x[num] - r0x * dF_z[num]);
dM_z[num] = (r0x * dF_y[num] - r0y * dF_x[num]);
M[0] =M[0] + dM_x[num];
M[1] =M[1] + dM_y[num];
M[2] =M[2] + dM_z[num];
Cd = F[0] / (0.5 * *rho * vmo * vmo * ref_A);
Cl = F[1] / (0.5 * *rho * vmo * vmo * ref_A);
Cz = F[2] / (0.5 * *rho * vmo * vmo * ref_A);
CMx = M[0] / (0.5 * *rho * vmo * vmo * ref_A * L);
CMy = M[1] / (0.5 * *rho * vmo * vmo * ref_A * L);
CMz = M[2] / (0.5 * *rho * vmo * vmo * ref_A * L);
}
//桥函数部分
//------------------------------------------------------------------------------------------------------------
else
{
//临界努森数
double Kn_cont = 0.001;
double Kn_fm = 100;
double h_cont;
double h_fm;
//桥函数系数
double Pb;
Pb = sin(PI * ( 0.375 + 0.125 * log10(Kn) ));
//分别用牛顿法和自由分子流理论计算Kn_cont = 0.001对应的C_cont和Kn_fm = 10对应的C_fm
double lamda_cont = Kn_cont * L;
double lamda_fm = Kn_fm * L;
//------------------------------------------------------------------------------------------------------------------------
//------------------------------------------------------------------------------------------------------------------------
env_calc_h(lamda_cont,&h_cont); //调用高度反算函数
env_calc_atmos(h_cont, rho, p, T, lamda, Me); //计算大气参数
// 来流马赫数
R = 8.31 * 1000 / *Me;
Ma = vmo/sqrt(gama * R * *T);
Ma_face[num] = V1mo/sqrt(gama * R * *T); //当地马赫数
S_face[num] = V1mo/sqrt(2 * R * *T); //当地速度比
if (Ma_face[num]>=7)//牛顿公式
{
if(test1>0)
{
Cp_facecont[num] = 0;
}
else if(test1<0)
{
Cp_facecont[num] = 2*(sin_face[num])*(sin_face[num]);
}
else if (test1==0)
{
if(fabs(u1+normal_x[num])>fabs(u1))
{
Cp_facecont[num] = 0;
}
else
{
Cp_facecont[num] = 2;
}
}
Ctao_facecont[num] = 0;
dF_x[num] = (0.5* *rho *(V1mo * V1mo) * Cp_facecont[num]+*p)*A_face[num] * normal_x[num];
dF_y[num] = (0.5* *rho *(V1mo * V1mo) * Cp_facecont[num]+*p)*A_face[num] * normal_y[num];
dF_z[num] = (0.5* *rho *(V1mo * V1mo) * Cp_facecont[num]+*p)*A_face[num] * normal_z[num];
F[0] = F[0] + dF_x[num];
F[1] = F[1] + dF_y[num];
F[2] = F[2] + dF_z[num];
dM_x[num] = (r0z * dF_y[num] - r0y * dF_z[num]);
dM_y[num] = (r0z * dF_x[num] - r0x * dF_z[num]);
dM_z[num] = (r0x * dF_y[num] - r0y * dF_x[num]);
M[0] =M[0] + dM_x[num];
M[1] =M[1] + dM_y[num];
M[2] =M[2] + dM_z[num];
Cd_cont = F[0] / (0.5 * *rho * vmo * vmo * ref_A);
Cl_cont = F[1] / (0.5 * *rho * vmo * vmo * ref_A);
Cz_cont = F[2] / (0.5 * *rho * vmo * vmo * ref_A);
CMx_cont = M[0] / (0.5 * *rho * vmo * vmo * ref_A * L);
CMy_cont = M[1] / (0.5 * *rho * vmo * vmo * ref_A * L);
CMz_cont = M[2] / (0.5 * *rho * vmo * vmo * ref_A * L);
}
else //修正的牛顿公式
{
if(test1>0)
{
Cp_facecont[num] = 0;
}
if(test1<0)
{
Cp_facecont[num] = 2 / gama / pow(Ma_face[num],2) * (pow(pow(gama+1,2)*pow(Ma_face[num],2)/(4*gama*pow(Ma_face[num],2)-2*gama+2), gama/(gama-1))*(1-gama+2*gama*pow(Ma_face[num],2))/(gama+1)-1) * pow(sin_face[num],2);
}
if (test1==0)
{
if(fabs(u1+normal_x[num])>fabs(u1))
{
Cp_facecont[num] = 0;
}
else
{
Cp_facecont[num] = 2 / gama / pow(Ma_face[num],2) * (pow(pow(gama+1,2)*pow(Ma_face[num],2)/(4*gama*pow(Ma_face[num],2)-2*gama+2), gama/(gama-1))*(1-gama+2*gama*pow(Ma_face[num],2))/(gama+1)-1);
}
}
Ctao_facecont[num] = 0;
dF_x[num] = (0.5* *rho *(V1mo * V1mo) * cp_face[num]+*p)*A_face[num] * normal_x[num];
dF_y[num] = (0.5* *rho *(V1mo * V1mo) * cp_face[num]+*p)*A_face[num] * normal_y[num];
dF_z[num] = (0.5* *rho *(V1mo * V1mo) * cp_face[num]+*p)*A_face[num] * normal_z[num];
F[0] = F[0] + dF_x[num];
F[1] = F[1] + dF_y[num];
F[2] = F[2] + dF_z[num];
dM_x[num] = (r0z * dF_y[num] - r0y * dF_z[num]);
dM_y[num] = (r0z * dF_x[num] - r0x * dF_z[num]);
dM_z[num] = (r0x * dF_y[num] - r0y * dF_x[num]);
M[0] =M[0] + dM_x[num];
⌨️ 快捷键说明
复制代码
Ctrl + C
搜索代码
Ctrl + F
全屏模式
F11
切换主题
Ctrl + Shift + D
显示快捷键
?
增大字号
Ctrl + =
减小字号
Ctrl + -