涡格法代码及解释

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

2、/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);

3、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)/次数限制,

4、精度要求,此处可修改,是迭带开关 for(i=0;i<n;i+)Mi = 0; for(j=0;j<n;j+)if(i!=j)Mi+=aij*xj;xi = (bi - Mi)/aii; /迭代 cout<<+t;for(i=0;i<n;i+)if(i%5=0)cout<<endl;cout<<" "<<xi;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.2

5、2505,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/18

6、0)/2>0)cout<<airfoil.L<<" "<<airfoil.Bg<<" "<<airfoil.Xo<<" "<<airfoil.Xc<<" "<<endl;break;elsecout<<"翼型的稍弦为0!请重新输入翼型数据"<<endl;cout<<"请输入来流马赫数和攻角"<<endl;cin>

7、;>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 doubleNq+1;double *baseB=ne

8、w 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;i<Ng;i+)girdlefti=new GIRDNq; girdright=new GIRD*Ng;for(i=0;i<Ng;i+)girdrighti=new GIRDNq;double width=airfoil.L/Nq/2;/展长每个分

9、块的长度/前缘节点的x坐标cout<<"前缘节点处的x坐标"<<endl;for(i=0;i<Nq+1;i+)baseqi=0+i*width*tan(airfoil.Xo*PI/180);cout<<baseqi<<" "<<endl;/每一条平行于根弦的弦的长度cout<<"每一条平行于根弦的弦的长度"<<endl;for(i=0;i<Nq+1;i+)baseBi=airfoil.Bg-i*(tan(airfoil.Xo*PI/180)

10、+tan(airfoil.Xc*PI/180)*width;cout<<baseBi<<" "<<endl;for(i=0;i<Ng;i+)for(j=0;j<Nq;j+)girdleftij.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.z

11、1=-1*girdleftij.z1;girdleftij.z3=girdleftij.z1;girdrightij.z3=-1*girdleftij.z3;girdleftij.z2=girdleftij.z1+width;girdrightij.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=g

12、irdleftij.x2+baseBj+1/2/Ng;girdrightij.x4=girdleftij.x4;girdleftij.x=(girdleftij.x3+girdleftij.x4)/2;girdrightij.x=girdleftij.x;girdleftij.z=(girdleftij.z3+girdleftij.z4)/2;girdrightij.z=-1*girdleftij.z;cout<<"*left*"<<endl;cout<<"(x1,z1):"<<"("&

13、lt;<girdleftij.x1<<","<<girdleftij.z1<<")"<<" "/将坐标打出cout<<"(x2,z2):"<<"("<<girdleftij.x2<<","<<girdleftij.z2<<")"<<endl;cout<<"(x3,z3):"<<

14、"("<<girdleftij.x3<<","<<girdleftij.z3<<")"<<" "cout<<"(x4,z4):"<<"("<<girdleftij.x4<<","<<girdleftij.z4<<")"<<" "cout<<"(x,z):

15、"<<"("<<girdleftij.x<<","<<girdleftij.z<<")"<<endl;cout<<"*right*"<<endl;cout<<"(x1,z1):"<<"("<<girdrightij.x1<<","<<girdrightij.z1<<")&

16、quot;<<" "/将坐标打出cout<<"(x2,z2):"<<"("<<girdrightij.x2<<","<<girdrightij.z2<<")"<<endl;cout<<"(x3,z3):"<<"("<<girdrightij.x3<<","<<girdrightij.

17、z3<<")"<<" "cout<<"(x4,z4):"<<"("<<girdrightij.x4<<","<<girdrightij.z4<<")"<<" "cout<<"(x,z):"<<"("<<girdrightij.x<<","<

18、;<girdrightij.z<<")"<<endl;/存储系数矩阵double *array;array=new double*2*Ng*Nq;for(i=0;i<2*Ng*Nq;i+)arrayi=new double2*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;arrayxy=vec(girdleftlk.x,girdleftlk.z,girdleftnm.x1,girdleftnm.z1,gir

19、dleftnm.x2,girdleftnm.z2);arrayxy+1=vec(girdleftlk.x,girdleftlk.z,girdrightnm.x1,girdrightnm.z1,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,

20、girdrightnm.x2,girdrightnm.z2);cout<<"*方程组系数矩阵*"<<endl;for(i=0;i<2*Ng*Nq;i+)for(j=0;j<2*Ng*Nq;j+)cout<<arrayij<<" "cout<<endl;cout<<"*线性方程组的右端项*"<<endl;for(i=0;i<2*Ng*Nq;i+)bi=-1*340*M*a;cout<<bi<<endl;cout&

21、lt;<"*Gauss-seidel法解线性方程组迭代20步的结果(每个涡格的环量)*"<<endl;for(i=0;i<2*Ng*Nq;i+)resulti=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*result2*i;airfoil.S=(baseB0+baseBNq)*airfoil.L/2;airfoil.Cy=2*airfoil.Y/p/pow(M*340,2)/airfoil.S;

22、cout<<"Y="<<airfoil.Y<<'t'<<"Cy="<<airfoil.Cy<<endl;为了验证代码的正确性，此处的算例采用的是空气动力学一书中关于涡格法的一道算例，书中给出了算例的过程和解。*运行结果*这是一个用涡格法计算机翼升力的程序!请输入翼型个参数:展长L, 根弦Bg,前缘后掠角Xo,后缘后掠角Xc5145-455 1 45 -45请输入来流马赫数和攻角0.210.2 0.0174533请输入根弦上的节点数,前缘上的节点数:252 5前缘节点处

23、的x坐标00.6251.251.8752.5每一条平行于根弦的弦的长度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)(x

24、3,z3):(1.375,0.625) (x4,z4):(2,1.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.2

25、5) (x2,z2):(2.125,-1.875)(x3,z3):(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.937

26、5,-2.1875)*方程组系数矩阵*-1.13826 -0.294675 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

27、 0.00780396 -0.179738 0.00398332 -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.0

28、259969 0.00769341 -0.0617391 0.00460806 -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

29、1.79981 -1.808883 1.93371 -1.97131 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.0