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

涡格法代码及解释

#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

{

M[i] = 0;

for(j=0;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

{

if(i%5==0){cout<

cout<<" "<

}

cout<

}

}

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<<"这是一个用涡格法计算机翼升力的程序!"<

cout<<"请输入翼型个参数:展长L, 根弦Bg,前缘后掠角Xo,后缘后掠角Xc"<

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<

break;

}

else

{

cout<<"翼型的稍弦为0!请重新输入翼型数据"<

}

}

cout<<"请输入来流马赫数和攻角"<

cin>>M>>a;

a=a*PI/180;

cout<

cout<<"请输入根弦上的节点数,前缘上的节点数:"<

cin>>Ng>>Nq;

cout<

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

{

girdleft[i]=new GIRD[Nq];

}

girdright=new GIRD*[Ng];

for(i=0;i

{

girdright[i]=new GIRD[Nq];

}

double width=airfoil.L/Nq/2;//展长每个分块的长度

//前缘节点的x坐标