wogefa.cpp

涡格法是空气动力学中的常用方法之一

#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(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;

}