渦格法代碼及解釋

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:/ 網(wǎng)格信息double x1,z1,x2,z2;/左右自由渦的坐標(biāo)double x3,z3,x4,z4;/3/4弦線處的坐標(biāo)double x,z;/控制點(diǎn)的坐標(biāo) ,3/4弦

2、線中點(diǎn)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)/迭代次數(shù)次數(shù)限制 , 精度要求, 此處可修改, 是迭帶開關(guān)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 這是一個(gè)用渦格法計(jì)算機(jī)翼升力的程序!endl;cout 請輸入翼型個(gè)參數(shù): 展長 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! 請重新輸入翼型數(shù)據(jù)endl;cout 請輸入來流馬赫數(shù)和攻角Ma;a=a*PI/180;coutMtaendl;cout 請輸入根弦上的節(jié)點(diǎn)數(shù), 前緣上的節(jié)點(diǎn)數(shù):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;左半邊機(jī)翼，右半邊機(jī)翼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;/展長每個(gè)分塊的長度/ 前緣節(jié)點(diǎn)的x 坐標(biāo)cout 前緣節(jié)點(diǎn)處的for(i=0;iNq+1;i+)x 坐標(biāo) 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);/將坐標(biāo)打出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);/將坐標(biāo)打出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;/ 存儲(chǔ)系數(shù)矩陣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*方程組系數(shù)矩陣*endl;for(i=0;i2*Ng*Nq;i+

14、)for(j=0;j2*Ng*Nq;j+)coutarrayij;coutendl;cout*線性方程組的右端項(xiàng)*endl;for(i=0;i2*Ng*Nq;i+)bi=-1*340*M*a;coutbiendl;cout*Gauss-seidel法解線性方程組迭代20 步的結(jié)果( 每個(gè)渦格的環(huán)量)*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;為了驗(yàn)證代碼的正確性，此處的算例采用的是空氣動(dòng)力學(xué)一書中關(guān)于渦格法的一道算例，書中給出了算例的過程和解。*運(yùn)行結(jié)果 *這是一個(gè)用渦格法計(jì)算機(jī)翼升力的程序!請輸入翼型個(gè)參數(shù): 展長 L,根弦 Bg, 前緣后掠角Xo,后緣后掠角Xc5145-455145-45請輸入來流馬赫數(shù)和攻角0.210.20.0174533請輸入根弦上的節(jié)點(diǎn)數(shù), 前緣上的節(jié)點(diǎn)數(shù):2525前緣節(jié)點(diǎn)處的x

16、坐標(biāo)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、數(shù)矩陣 *-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、方程組的右端項(xiàng)*-1.18682-1.18682-1.18682-1.18682-1.18682-1.18682-1.18682-1.18682*Gauss-seidel法解線性方程組迭代20 步的結(jié)果 ( 每個(gè)渦格的環(huán)量)*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