涡格法代码及解释教程.doc

#include iostream.h #include stdio.h #include math.h #define PI 3.1415926 class AIRFOIL //用来存放翼型的信息 { public: double L,Bg,S; double Xo,Xc; double Y,Cy; AIRFOIL(){Y=0.0f,S=0.0f,L=0.0f,Bg=0.0f,Xo=0.0f,Xc=0.0f;} }; class GIRD //网格信息 { public: double x1,z1,x2,z2;//左右自由涡的坐标 double x3,z3,x4,z4;//3/4弦线处的坐标 double x,z;//控制点的坐标,3/4弦线中点 GIRD(){x1=0.0f,x2=0.0f,z1=0.0f,z2=0.0f,x3=0.0f,x4=0.0f,z3=0.0f,z4=0.0f,x=0.0f,z=0.0f;} }; double vec(double x,double z,double x1,double z1,double x2,double z2 ) { double a,b,c,d,e; a=1/((x2-x)*(z1-z)-(x1-x)*(z2-z)); b=((x2-x1)*(x1-x)+(z2-z1)*(z1-z))/sqrt(pow((x1-x),2)+pow((z1-z),2)); c=((x2-x1)*(x2-x)+(z2-z1)*(z2-z))/sqrt(pow((x2-x),2)+pow((z2-z),2)); d=(1-(x1-x)/sqrt(pow((x1-x),2)+pow((z1-z),2)))/(z1-z); e=(1-(x2-x)/sqrt(pow((x2-x),2)+pow((z2-z),2)))/(z2-z); return (a*(b-c)+d-e)/4/PI; } void Gaussseidel(int n,double *M,double **a,double *x,double *b)//高斯--塞得尔迭带法 { int t=0,i,j;//迭代次数 while(t20)//次数限制,精度要求,此处可修改,是迭带开关 { for(i=0;in;i++) { M[i] = 0; for(j=0;jn;j++) { if(i!=j) { M[i]+=a[i][j]*x[j]; } } x[i] = (b[i] - M[i])/a[i][i]; //迭代 } cout++t; for(i=0;in;i++) { if(i%5==0){coutendl;} cout x[i]; } coutendl; } } void main() { AIRFOIL airfoil; int Ng,Nq,i,j,k,l,m,n,x,y; double Y=0.0,M,a,ep=1e-10,p=1.22505,Cy=0.0; //p为海平面空气密度 cout这是一个用涡格法计算机翼升力的程序!endl; cout请输入翼型个参数:展长L, 根弦Bg,前缘后掠角Xo,后缘后掠角Xcendl; while(1) { cinairfoil.Lairfoil.Bgairfoil.Xoairfoil.Xc; if(airfoil.Bg-airfoil.L*(tan(airfoil.Xo*PI/180)+tan(airfoil.Xc*PI/180))/20) { coutairfoil.L airfoil.Bg airfoil.Xo airfoil.Xc endl; break; } else { cout翼型的稍弦为0!请重新输入翼型数据endl; } } cout请输入来流马赫数和攻角endl; cinMa; a=a*PI/180; coutM\taendl; cout请输入根弦上的节点数,前缘上的节点数:endl; cinNgNq; coutNg Nq endl; Nq--;Ng--;//变成分多少块 double *baseq=new double[Nq+1]; double *baseB=new double[Nq+1]; double *result=new double[2*Nq*Ng]; double *b=new double[2*Nq*Ng]; double *M1=new double[2*Nq*Ng