涡格法代码及解释

#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 弦线中点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(t20)/ 次数限制 ,精度要求 ,此处可修改 ,是迭带开关 for(i=0;in;i+) Mi = 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; 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;Xcendl;airfoil.Xccoutairfoil.Lairfoil.Bgairfoil.Xoairfoil.Xc;if(airfoil.Bg-airfoil.L*(tan(airfoil.Xo*PI/180)+tan(airfoil.Xc*PI/180)/20) coutairfoil.Lairfoil.Bgairfoil.Xoendl;break;elsecout 翼型的稍弦为 0!请重新输入翼型数据 endl;cout 请输入来流马赫数和攻角 Ma;a=a*PI/180;coutMtaendl;cout 请输入根弦上的节点数,前缘上的节点数 ：NgNq;coutNg Nq endl;Nq-;Ng-;/ 变成分多少块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=new GIRDNq;double width=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+) 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=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+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):(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,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 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,girdleftnm.z1,girdleftnm.x 2,girdleftnm.z2);arrayxy+1=vec(girdleftlk.x,girdleftlk.z,girdrightnm.x1,girdrightnm.z1,girdright nm.x2,girdrightnm.z2);arrayx+1y=vec(girdrightlk.x,girdrightlk.z,girdleftnm.x1,girdleftnm.z1,girdleftn m.x2,girdleftnm.z2);arrayx+1y+1=vec(girdrightlk.x,girdrightlk.z,girdrightnm.x1,girdrightnm.z1,girdr ightnm.x2,girdrightnm.z2);cout*方程组系数矩阵*endl;for(i=0;i2*Ng*Nq;i+) 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.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,后缘后掠角Xc145-455 1 45 -45 请输入来流马赫数和攻角0.210.2 0.0174533请输入根弦上的节点数 ,前缘上的节点数252 5前缘节点处的 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)* right* *(x1,z1):(0.25,0)(x2,z2):(0.875,-0.625)(x3,z3):(0.75,0)(x4,z4):(1.375,-0.625)* left* *(x,z):(1.0625,0.3125)(x,z):(1.0625,-0.3125)(x1,z1):(0.875,0.625) (x2,z2):(1.5,1.25)(x3,z3):(1.375,0.625) (x4,z4):(2,1.25) * right* * (x1,z1):(0.875,-0.625)(x2,z2):(1.5,-1.25)(x3,z3):(1.375,-0.625)(x4,z4):(2,-1.25)* left* *(x,z):(1.6875,0.9375)(x,z):(1.6875,-0.9375)(x1,z1):(1.5,1.25) (x2,z2):(2.125,1.875) (x3,z3):(2,1.25) (x4,z4):(2.625,1.875) * right* * (x1,z1):(1.5,-1.25) (x2,z2):(2.125,-1.875) (x3,z3):(2,-1.25) (x4,z4):(2.625,-1.875) * left* *(x,z):(2.3125,1.5625)(x,z):(2.3125,-1.5625)(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.294675 0.179738-0.03263340.0171196-0.009369350.00600848-0.004230970.294675 1.13826 0.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.0186878 1.138260.00780396-0.1797380.00398332-0.01711960.0617391 -0.02463680.32177 -0.0115021-1.13826-0.006009450.179738-0.003467210.0246368 -0.06173910.0115021 -0.321770.00600945 1.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.00292199 1.13826*线性方程组的右端项 *-1.18682 -1.18682 -1.18682 -1.18682 -1.18682 -1.18682 -1.18682-1.18682*Gauss-seidel法解线性方程组迭代20 步的结果 （每个涡格的环量 ）*11.04267-1.5798-1.31261.610081.40375-1.63243-1.489461.5395121.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-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.13526 1.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.035932.14261-2.142612.13538-2.135381.86144-1.86144112.03594-2.035942.14261-2.142612.13538-2.135381.86144-1.86144122.03594-2.035942.14261-2.142612.13538-2.135381.86144-1.86144132.03594-2.035942.14261-2.142612.13538-2.135381.86144-1.86144142.03594-2.035942.14261-2.142612.13538-2.135381.86144-1.86144152.03594-2.035942.14261-2.142612.13538-2.135381.86144-1.86144162.03594-2.035942.14261-2.142612.13538-2.135381.86144-1.86144172.03594-2.035942.14261-2.142612.13538-2.135381.86144-1.86144182.03594-2.035942.14261-2.142612.13538-2.135381.86144-1.86144192.03594-2.035942.14261-2.142612.13538-2.135381.86144-1.86144202.03594-2.035942.14261-2.142612.13538-2.135381.86144-1.86144Y=851.296 Cy=0.0601131在同样的条件下，空气动力学书中给出的结果为：Y=848.554 Cy=0.059919可见程序正确