涡格法代码及解释

1、#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弦线中点GI

2、RD()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

3、(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+)Mi = 0; f

4、or(j=0;jn;j+)if(i!=j)Mi+=aij*xj;xi = (bi - Mi)/aii; /迭代 cout+t;for(i=0;in;i+)if(i%5=0)coutendl;cout xi;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,后缘后掠角Xcairfoil.Lairfoil

5、.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;elsecout翼型的稍弦为0!请重新输入翼型数据endl;cout请输入来流马赫数和攻角Ma;a=a*PI/180;coutMtaendl;cout请输入根弦上的节点数,前缘上的节点数:NgNq;coutNg Nq endl;Nq-;Ng-;/变成分多少块double *baseq=n

6、ew doubleNq+1;double *baseB=new doubleNq+1;double *result=new double2*Nq*Ng;double *b=new double2*Nq*Ng;double *M1=new double2*Nq*Ng;GIRD *girdleft,*girdright;/左半边机翼，右半边机翼girdleft=new GIRD*Ng;for(i=0;iNg;i+)girdlefti=new GIRDNq; girdright=new GIRD*Ng;for(i=0;iNg;i+)girdrighti=new GIRDNq;double width

7、=airfoil.L/Nq/2;/展长每个分块的长度/前缘节点的x坐标cout前缘节点处的x坐标endl;for(i=0;iNq+1;i+)baseqi=0+i*width*tan(airfoil.Xo*PI/180);coutbaseqi endl;/每一条平行于根弦的弦的长度cout每一条平行于根弦的弦的长度endl;for(i=0;iNq+1;i+)baseBi=airfoil.Bg-i*(tan(airfoil.Xo*PI/180)+tan(airfoil.Xc*PI/180)*width;coutbaseBi endl;for(i=0;iNg;i+)for(j=0;jNq;j+)gi

8、rdleftij.x1=baseqj+baseBj/4/Ng+i*baseBj/Ng;girdrightij.x1=girdleftij.x1;girdleftij.x3=girdleftij.x1+baseBj/2/Ng;girdrightij.x3=girdleftij.x3;girdleftij.z1=0+j*width;girdrightij.z1=-1*girdleftij.z1;girdleftij.z3=girdleftij.z1;girdrightij.z3=-1*girdleftij.z3;girdleftij.z2=girdleftij.z1+width;girdright

9、ij.z2=-1*girdleftij.z2;girdleftij.z4=girdleftij.z2;girdrightij.z4=-1*girdleftij.z4;girdleftij.x2=baseqj+1+baseBj+1/4/Ng+i*baseBj+1/Ng;girdrightij.x2=girdleftij.x2;girdleftij.x4=girdleftij.x2+baseBj+1/2/Ng;girdrightij.x4=girdleftij.x4;girdleftij.x=(girdleftij.x3+girdleftij.x4)/2;girdrightij.x=girdlef

10、tij.x;girdleftij.z=(girdleftij.z3+girdleftij.z4)/2;girdrightij.z=-1*girdleftij.z;cout*left*endl;cout(x1,z1):(girdleftij.x1,girdleftij.z1) ;/将坐标打出cout(x2,z2):(girdleftij.x2,girdleftij.z2)endl;cout(x3,z3):(girdleftij.x3,girdleftij.z3) ;cout(x4,z4):(girdleftij.x4,girdleftij.z4) ;cout(x,z):(girdleftij.x

11、,girdleftij.z)endl;cout*right*endl;cout(x1,z1):(girdrightij.x1,girdrightij.z1) ;/将坐标打出cout(x2,z2):(girdrightij.x2,girdrightij.z2)endl;cout(x3,z3):(girdrightij.x3,girdrightij.z3) ;cout(x4,z4):(girdrightij.x4,girdrightij.z4) ;cout(x,z):(girdrightij.x,girdrightij.z)endl;/存储系数矩阵double *array;array=new d

12、ouble*2*Ng*Nq;for(i=0;i2*Ng*Nq;i+)arrayi=new double2*Ng*Nq;for(i=0;iNq*Ng;i+)k=i%Nq;l=i/Nq;for(j=0;jNq*Ng;j+)m=j%Nq;n=j/Nq;x=2*i;y=2*j;arrayxy=vec(girdleftlk.x,girdleftlk.z,girdleftnm.x1,girdleftnm.z1,girdleftnm.x2,girdleftnm.z2);arrayxy+1=vec(girdleftlk.x,girdleftlk.z,girdrightnm.x1,girdrightnm.z1,

13、girdrightnm.x2,girdrightnm.z2);arrayx+1y=vec(girdrightlk.x,girdrightlk.z,girdleftnm.x1,girdleftnm.z1,girdleftnm.x2,girdleftnm.z2);arrayx+1y+1=vec(girdrightlk.x,girdrightlk.z,girdrightnm.x1,girdrightnm.z1,girdrightnm.x2,girdrightnm.z2);cout*方程组系数矩阵*endl;for(i=0;i2*Ng*Nq;i+)for(j=0;j2*Ng*Nq;j+)coutarr

14、ayij ;coutendl;cout*线性方程组的右端项*endl;for(i=0;i2*Ng*Nq;i+)bi=-1*340*M*a;coutbiendl;cout*Gauss-seidel法解线性方程组迭代20步的结果(每个涡格的环量)*endl;for(i=0;i2*Ng*Nq;i+)resulti=0.0;Gaussseidel(2*Nq*Ng,M1,array,result,b);for(i=0;iNg*Nq;i+)airfoil.Y=airfoil.Y+2*p*M*340*width*result2*i;airfoil.S=(baseB0+baseBNq)*airfoil.L/2

15、;airfoil.Cy=2*airfoil.Y/p/pow(M*340,2)/airfoil.S;coutY=airfoil.YtCy=airfoil.Cyendl;为了验证代码的正确性，此处的算例采用的是空气动力学一书中关于涡格法的一道算例，书中给出了算例的过程和解。*运行结果*这是一个用涡格法计算机翼升力的程序!请输入翼型个参数:展长L, 根弦Bg,前缘后掠角Xo,后缘后掠角Xc5145-455 1 45 -45请输入来流马赫数和攻角0.210.2 0.0174533请输入根弦上的节点数,前缘上的节点数:252 5前缘节点处的x坐标00.6251.251.8752.5每一条平行于根弦的弦的

16、长度11111*left*(x1,z1):(0.25,0) (x2,z2):(0.875,0.625)(x3,z3):(0.75,0) (x4,z4):(1.375,0.625) (x,z):(1.0625,0.3125)*right*(x1,z1):(0.25,0) (x2,z2):(0.875,-0.625)(x3,z3):(0.75,0) (x4,z4):(1.375,-0.625) (x,z):(1.0625,-0.3125)*left*(x1,z1):(0.875,0.625) (x2,z2):(1.5,1.25)(x3,z3):(1.375,0.625) (x4,z4):(2,1.

17、25) (x,z):(1.6875,0.9375)*right*(x1,z1):(0.875,-0.625) (x2,z2):(1.5,-1.25)(x3,z3):(1.375,-0.625) (x4,z4):(2,-1.25) (x,z):(1.6875,-0.9375)*left*(x1,z1):(1.5,1.25) (x2,z2):(2.125,1.875)(x3,z3):(2,1.25) (x4,z4):(2.625,1.875) (x,z):(2.3125,1.5625)*right*(x1,z1):(1.5,-1.25) (x2,z2):(2.125,-1.875)(x3,z3):

18、(2,-1.25) (x4,z4):(2.625,-1.875) (x,z):(2.3125,-1.5625)*left*(x1,z1):(2.125,1.875) (x2,z2):(2.75,2.5)(x3,z3):(2.625,1.875) (x4,z4):(3.25,2.5) (x,z):(2.9375,2.1875)*right*(x1,z1):(2.125,-1.875) (x2,z2):(2.75,-2.5)(x3,z3):(2.625,-1.875) (x4,z4):(3.25,-2.5) (x,z):(2.9375,-2.1875)*方程组系数矩阵*-1.13826 -0.29

19、4675 0.179738 -0.0326334 0.0171196 -0.00936935 0.00600848 -0.004230970.294675 1.13826 0.0326334 -0.179738 0.00936935 -0.0171196 0.00423097 -0.006008480.32177 -0.0575242 -1.13826 -0.0186878 0.179738 -0.00780396 0.0171196 -0.003983320.0575242 -0.32177 0.0186878 1.13826 0.00780396 -0.179738 0.00398332

20、-0.01711960.0617391 -0.0246368 0.32177 -0.0115021 -1.13826 -0.00600945 0.179738 -0.003467210.0246368 -0.0617391 0.0115021 -0.32177 0.00600945 1.13826 0.00346721 -0.1797380.0259969 -0.0136999 0.0617391 -0.00769341 0.32177 -0.00460806 -1.13826 -0.002921990.0136999 -0.0259969 0.00769341 -0.0617391 0.00

21、460806 -0.32177 0.00292199 1.13826*线性方程组的右端项*-1.18682-1.18682-1.18682-1.18682-1.18682-1.18682-1.18682-1.18682*Gauss-seidel法解线性方程组迭代20步的结果(每个涡格的环量)*1 1.04267 -1.3126 1.40375 -1.48946 1.53951 -1.5798 1.61008 -1.632432 1.69757 -1.80862 1.92227 -1.96149 2.00467 -2.02468 1.79981 -1.808883 1.93371 -1.9713

22、1 2.08496 -2.09885 2.10122 -2.10804 1.845 -1.848064 2.00797 -2.01931 2.12727 -2.13144 2.12628 -2.12835 1.85706 -1.857985 2.02866 -2.03164 2.13859 -2.13972 2.13299 -2.13356 1.86029 -1.860546 2.03405 -2.03482 2.14156 -2.14186 2.13476 -2.13491 1.86114 -1.861217 2.03545 -2.03565 2.14234 -2.14242 2.13522 -2.13526 1.86136 -1.861388 2.03581 -2.03586 2.14254 -2.14256 2.13534 -2.13535 1.86142 -1.861429 2.03591 -2.03592 2.14259 -2.1426 2.13537 -2.13537 1.86143 -1.8614310 2.03593 -2.03593 2.14261 -2.14261 2.13538 -2.13538 1.86144 -1.8614411 2.03594 -2.03594 2.14261 -2.14261 2.13538 -2.13