涡格法代码及解释

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 GIRDpublic:/ 网格信息double x1,z1,x2,z2;/左右自由涡的坐标double x3,z3,x4,z4;/3/4弦线处的坐标double x,z;/控制点的坐标 ,3/4弦

2、线中点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)

3、/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(t20)/迭代次数次数限制 , 精度要求, 此处可修改, 是迭带开关for(i=0;in;i+)M

4、i = 0;for(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;coutxi;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; /pcout 这是一个用涡格法计算机翼升力的程序!endl;cout 请输入翼型个参数: 展长 L,根弦 Bg, 前缘后掠角while(1)为海平面空气密度Xo,后缘后掠角

5、Xcairfoil.Lairfoil.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;coutNgNqendl;Nq-;Ng-

6、;/变成分多少块double *baseq=new 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=n

7、ew GIRDNq;double width=airfoil.L/Nq/2;/展长每个分块的长度/ 前缘节点的x 坐标cout 前缘节点处的for(i=0;iNq+1;i+)x 坐标 endl;baseqi=0+i*width*tan(airfoil.Xo*PI/180);coutbaseqiendl;/ 每一条平行于根弦的弦的长度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

8、(i=0;iNg;i+)for(j=0;jNq;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.z1=-1*girdleftij.z1;girdleftij.z3=girdleftij.z1;girdrightij.z3=-1*girdleftij.z3;girdleftij.z2=

9、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=girdleftij.x2+baseBj+1/2/Ng;girdrightij.x4=girdleftij.x4;girdleftij.x=(girdleftij.x3+girdlefti

10、j.x4)/2;girdrightij.x=girdleftij.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.z

11、4);cout(x,z):(girdleftij.x,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;/ 存储系数矩阵d

12、ouble *array;array=new double*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,girdleftn m.z1,girdleftnm.x2,girdleftnm.z2);arrayxy+1=vec(girdleftlk.x,girdleftlk.z,girdr

13、ightnm.x1,girdrig htnm.z1,girdrightnm.x2,girdrightnm.z2);arrayx+1y=vec(girdrightlk.x,girdrightlk.z,girdleftnm.x1,girdle ftnm.z1,girdleftnm.x2,girdleftnm.z2);arrayx+1y+1=vec(girdrightlk.x,girdrightlk.z,girdrightnm.x1,gir drightnm.z1,girdrightnm.x2,girdrightnm.z2);cout*方程组系数矩阵*endl;for(i=0;i2*Ng*Nq;i+

14、)for(j=0;j2*Ng*Nq;j+)coutarrayij;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.

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

16、坐标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)(x3,z3):(

17、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.25) (x2,z2):(

18、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.9375,-2.1875)*方程组系

19、数矩阵 *-1.13826-0.2946750.179738-0.03263340.0171196-0.009369350.00600848-0.004230970.2946751.138260.0326334-0.1797380.00936935-0.01711960.00423097-0.006008480.32177-0.0575242-1.13826-0.01868780.179738-0.007803960.0171196-0.003983320.0575242-0.321770.01868781.138260.00780396-0.1797380.00398332-0.017119

20、60.0617391-0.02463680.32177-0.0115021-1.13826-0.006009450.179738-0.003467210.0246368-0.06173910.0115021-0.321770.006009451.138260.00346721-0.1797380.0259969-0.01369990.0617391-0.007693410.32177-0.00460806-1.13826-0.002921990.0136999-0.02599690.00769341-0.06173910.00460806-0.321770.002921991.13826*线性

21、方程组的右端项*-1.18682-1.18682-1.18682-1.18682-1.18682-1.18682-1.18682-1.18682*Gauss-seidel法解线性方程组迭代20 步的结果 ( 每个涡格的环量)*11.04267-1.31261.40375-1.489461.53951-1.57981.61008-1.6324321.69757-1.808621.92227-1.961492.00467-2.024681.79981-1.8088831.93371-1.971312.08496-2.098852.10122-2.108041.845-1.8480642.00797

22、-2.019312.12727-2.131442.12628-2.128351.85706-1.8579852.02866-2.031642.13859-2.139722.13299-2.133561.86029-1.8605462.03405-2.034822.14156-2.141862.13476-2.134911.86114-1.8612172.03545-2.035652.14234-2.142422.13522-2.135261.86136-1.8613882.03581-2.035862.14254-2.142562.13534-2.135351.86142-1.8614292.03591-2.035922.14259-2.14262.13537-2.135371.86143-1.86143102.03593-2