涡格法代码及解释培训讲学

涡格法代码及解释#include "iostream.h"#include "stdio.h"#include "math.h"#define PI 3.1415926class 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(t<20)//次数限制,精度要求,此处可修改,是迭带开关{for(i=0;i<n;i++){M[i] = 0;for(j=0;j<n;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;i<n;i++){if(i%5==0){cout<<endl;}cout<<" "<<x[i];}cout<<endl;}}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,后缘后掠角Xc"<<endl;while(1){cin>>airfoil.L>>airfoil.Bg>>airfoil.Xo>>airfoil.Xc;if(airfoil.Bg-airfoil.L*(tan(airfoil.Xo*PI/180)+tan(airfoil.Xc*PI/180))/2>0){cout<<airfoil.L<<" "<<airfoil.Bg<<" "<<airfoil.Xo<<" "<<airfoil.Xc<<" "<<endl;break;}else{cout<<"翼型的稍弦为0!请重新输入翼型数据"<<endl;}}cout<<"请输入来流马赫数和攻角"<<endl;cin>>M>>a;a=a*PI/180;cout<<M<<'\t'<<a<<endl;cout<<"请输入根弦上的节点数,前缘上的节点数:"<<endl;cin>>Ng>>Nq;cout<<Ng<<" "<<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];GIRD **girdleft,**girdright;//左半边机翼，右半边机翼girdleft=new GIRD*[Ng];for(i=0;i<Ng;i++){girdleft[i]=new GIRD[Nq];}girdright=new GIRD*[Ng];for(i=0;i<Ng;i++){girdright[i]=new GIRD[Nq];}double width=airfoil.L/Nq/2;//展长每个分块的长度//前缘节点的x坐标cout<<"前缘节点处的x坐标"<<endl;for(i=0;i<Nq+1;i++){baseq[i]=0+i*width*tan(airfoil.Xo*PI/180);cout<<baseq[i]<<" "<<endl;}//每一条平行于根弦的弦的长度cout<<"每一条平行于根弦的弦的长度"<<endl;for(i=0;i<Nq+1;i++){baseB[i]=airfoil.Bg-i*(tan(airfoil.Xo*PI/180)+tan(airfoil.Xc*PI/180))*width;cout<<baseB[i]<<" "<<endl;}for(i=0;i<Ng;i++){for(j=0;j<Nq;j++){girdleft[i][j].x1=baseq[j]+baseB[j]/4/Ng+i*baseB[j]/Ng;girdright[i][j].x1=girdleft[i][j].x1;girdleft[i][j].x3=girdleft[i][j].x1+baseB[j]/2/Ng;girdright[i][j].x3=girdleft[i][j].x3;girdleft[i][j].z1=0+j*width;girdright[i][j].z1=-1*girdleft[i][j].z1;girdleft[i][j].z3=girdleft[i][j].z1;girdright[i][j].z3=-1*girdleft[i][j].z3;girdleft[i][j].z2=girdleft[i][j].z1+width;girdright[i][j].z2=-1*girdleft[i][j].z2;girdleft[i][j].z4=girdleft[i][j].z2;girdright[i][j].z4=-1*girdleft[i][j].z4;girdleft[i][j].x2=baseq[j+1]+baseB[j+1]/4/Ng+i*baseB[j+1]/Ng;girdright[i][j].x2=girdleft[i][j].x2;girdleft[i][j].x4=girdleft[i][j].x2+baseB[j+1]/2/Ng;girdright[i][j].x4=girdleft[i][j].x4;girdleft[i][j].x=(girdleft[i][j].x3+girdleft[i][j].x4)/2;girdright[i][j].x=girdleft[i][j].x;girdleft[i][j].z=(girdleft[i][j].z3+girdleft[i][j].z4)/2;girdright[i][j].z=-1*girdleft[i][j].z;cout<<"************************left**********************"<<endl;cout<<"(x1,z1):"<<"("<<girdleft[i][j].x1<<","<<girdleft[i][j].z1<<")"<<"";//将坐标打出cout<<"(x2,z2):"<<"("<<girdleft[i][j].x2<<","<<girdleft[i][j].z2<<")"<<endl;cout<<"(x3,z3):"<<"("<<girdleft[i][j].x3<<","<<girdleft[i][j].z3<<")"<<"";cout<<"(x4,z4):"<<"("<<girdleft[i][j].x4<<","<<girdleft[i][j].z4<<")"<<"";cout<<"(x,z):"<<"("<<girdleft[i][j].x<<","<<girdleft[i][j].z<<")"<<endl;cout<<"************************right**********************"<<endl;cout<<"(x1,z1):"<<"("<<girdright[i][j].x1<<","<<girdright[i][j].z1<<")"<<" ";//将坐标打出cout<<"(x2,z2):"<<"("<<girdright[i][j].x2<<","<<girdright[i][j].z2<<")"<<endl;cout<<"(x3,z3):"<<"("<<girdright[i][j].x3<<","<<girdright[i][j].z3<<")"<<" ";cout<<"(x4,z4):"<<"("<<girdright[i][j].x4<<","<<girdright[i][j].z4<<")"<<" ";cout<<"(x,z):"<<"("<<girdright[i][j].x<<","<<girdright[i][j].z<<")"<<endl;}}//存储系数矩阵double **array;array=new double*[2*Ng*Nq];for(i=0;i<2*Ng*Nq;i++){array[i]=new double[2*Ng*Nq];}for(i=0;i<Nq*Ng;i++){k=i%Nq;l=i/Nq;for(j=0;j<Nq*Ng;j++){m=j%Nq;n=j/Nq;x=2*i;y=2*j;array[x][y]=vec(girdleft[l][k].x,girdleft[l][k].z,girdleft[n][m].x1,girdleft[n][m].z1, girdleft[n][m].x2,girdleft[n][m].z2);array[x][y+1]=vec(girdleft[l][k].x,girdleft[l][k].z,girdright[n][m].x1,girdright[n][m ].z1,girdright[n][m].x2,girdright[n][m].z2);array[x+1][y]=vec(girdright[l][k].x,girdright[l][k].z,girdleft[n][m].x1,girdleft[n][m ].z1,girdleft[n][m].x2,girdleft[n][m].z2);array[x+1][y+1]=vec(girdright[l][k].x,girdright[l][k].z,girdright[n][m].x1,girdright [n][m].z1,girdright[n][m].x2,girdright[n][m].z2);}}cout<<"****************方程组系数矩阵***************"<<endl;for(i=0;i<2*Ng*Nq;i++){for(j=0;j<2*Ng*Nq;j++){cout<<array[i][j]<<" ";}cout<<endl;}cout<<"***************线性方程组的右端项*************"<<endl;for(i=0;i<2*Ng*Nq;i++){b[i]=-1*340*M*a;cout<<b[i]<<endl;}cout<<"***************Gauss-seidel法解线性方程组迭代20步的结果(每个涡格的环量)*************"<<endl;for(i=0;i<2*Ng*Nq;i++){result[i]=0.0;}Gaussseidel(2*Nq*Ng,M1,array,result,b);for(i=0;i<Ng*Nq;i++)airfoil.Y=airfoil.Y+2*p*M*340*width*result[2*i];airfoil.S=(baseB[0]+baseB[Nq])*airfoil.L/2;airfoil.Cy=2*airfoil.Y/p/pow(M*340,2)/airfoil.S;cout<<"Y="<<airfoil.Y<<'\t'<<"Cy="<<airfoil.Cy<<endl;}为了验证代码的正确性，此处的算例采用的是《空气动力学》一书中关于涡格法的一道算例，书中给出了算例的过程和解。