1、#includeiostream.h#includestdio.h#includemath.h#definePI3.1415926classAIRFOIL/用來存放翼型的信息public:doubleL,Bg,S;doubleXo,Xc;doubleY,Cy;AIRFOIL()Y=0.0f,S=0.0f,L=0.0f,Bg=0.0f,Xo=0.0f,Xc=0.0f;classGIRD/網(wǎng)格信息public:doublex1,z1,x2,z2;/左右自由渦的坐標(biāo)doublex3,z3,x4,z4;/3/4弦線處的坐標(biāo)doublex,z;/控制點的坐標(biāo),3/4弦線中點GIRD()x1=0.0f,x

2、2=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;doublevec(doublex,doublez,doublex1,doublez1,doublex2,doublez2)doublea,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

3、),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;voidGaussseidel(intn,double*M,double*a,double*x,double*切/高斯一塞得爾迭帶法intt=0,i,j；/迭代次數(shù)while(t20)/次數(shù)限制，精度要求，此處可修改，是迭帶開關(guān)for(i=0;in;i+)Mi=0;for(j=0;jn;j+)if(i!=j)(Mi+=aij*xj;)xi

4、=(bi-Mi)/aii;迭代)cout+t;for(i=0;in;i+)(if(i%5=0)coutendl;coutxi;coutendl;voidmain()AIRFOILairfoil;intNg,Nq,i,j,k,l,m,n,x,y;doubleY=0.0,M,a,ep=1e-10,p=1.22505,Cy=0.0;/p為海平面空氣密度cout”這是一個用渦格法計算機翼升力的程序!”endl;cout請輸入翼型個參數(shù):展長L,根弦Bg,前緣后掠角Xo,后緣后掠角Xcairfoil.Lairfoil.Bgairfoil.Xoairfoil.Xc;if(airfoil.Bg-airfoi

5、l.L*(tan(airfoil.Xo*PI/180)+tan(airfoil.Xc*PI/180)/20)coutairfoil.Lairfoil.Bgairfoil.Xoairfoil.Xcendl;break;else8a翼型的稍弦為0!請重新輸入翼型數(shù)據(jù)endl;cout請輸入來流馬赫數(shù)和攻角Ma;a=a*PI/180;coutMtaendl;cout請輸入根弦上的節(jié)點數(shù)，前緣上的節(jié)點數(shù):NgNq;coutNgNqendl;Nq-;Ng-;/變成分多少塊double*baseq=newdoubleNq+1;double*baseB=newdoubleNq+1;double*result

6、=newdouble2*Nq*Ng;double*b=newdouble2*Nq*Ng;double*M1=newdouble2*Nq*Ng;GIRD*girdleft,*girdright;/左半邊機翼，右半邊機翼girdleft=newGIRD*Ng;for(i=0;iNg;i+)girdlefti=newGIRDNq;girdright=newGIRD*Ng;for(i=0;iNg;i+)girdrighti=newGIRDNq;doublewidth=airfoil.L/Nq/2;/展長每個分塊的長度前緣節(jié)點的x坐標(biāo)cout前緣節(jié)點處的x坐標(biāo)endl;for(i=0;iNq+1;i+)

7、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;coutbaseBiendl;for(i=0;iNg;i+)for(j=0;jNq;j+)girdleftij.x1=baseqj+baseBj/4/Ng+i*baseBj/Ng;girdrightij.x1=girdleftij.x1;

8、girdleftij.x3=girdleftij.xl+baseBj/2/Ng;girdrightij.x3=girdleftij.x3;girdleftij.z1=0+j*width;girdrightij.zl=-l*girdleftij.zl;girdleftij.z3=girdleftij.zl;girdrightij.z3=-l*girdleftij.z3;girdleftij.z2=girdleftij.zl+width;girdrightij.z2=-l*girdleftij.z2;girdleftij.z4=girdleftij.z2;girdrightij.z4=-l*gir

9、dleftij.z4;girdleftij.x2=baseqj+l+baseBj+l/4/Ng+i*baseBj+l/Ng;girdrightij.x2=girdleftij.x2;girdleftij.x4=girdleftij.x2+baseBj+l/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=-l*gird

10、leftij.z;coutn*Wendi，cout(xl,zl):(girdleftij.xl,girdleftij.zl)，將坐標(biāo)打出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*”eiidl.cout(xl,zl):(girdrightij.xl,girdrightij.zl)

11、將坐標(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;存儲系數(shù)矩陣double*array;array=newdouble*2*Ng*Nq;for(i=0;iv2*Ng*Nq;i+)(arrayi=newdouble2*Ng*Nq;for(i=0;ivNq*Ng;i+)k=i%Nq;l=

12、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,girdrightnm.x2,girdrightnm.z2);arrayx+1y=vec(girdrightlk.x,girdrightlk.z,girdleftnm.x1

13、,girdleftnm.z1,girdleftnm.x2,girdleftnm.z2);arrayx+1y+1=vec(girdrightlk.x,girdrightlk.z,girdrightnm.x1,girdrightnm.z1,girdrightnm.x2,girdrightnm.z2);cout*方程組系數(shù)矩陣*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*Ga

14、uss-seidel法解線性方程組迭代20步的結(jié)果(每個渦格的環(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.S=(baseB0+baseBNq)*airfoil.L/2;airfoil.Cy=2*airfoil.Y/p/pow(M*340,2)/airfoil.S;coutY=airfoil.YtCy=airfoil.Cyendl;為了驗證

15、代碼的正確性，此處的算例采用的是空氣動力學(xué)一書中關(guān)于渦格法的一道算例，書中給出了算例的過程和解。*運行結(jié)果這是一個用渦格法計算機翼升力的程序!請輸入翼型個參數(shù):展長L,根弦Bg,前緣后掠角Xo,后緣后掠角Xc45-455145-45請輸入來流馬赫數(shù)和攻角0.210.20.0174533請輸入根弦上的節(jié)點數(shù),前緣上的節(jié)點數(shù):2525前緣節(jié)點處的x坐標(biāo)00.6251.251.8752.5每一條平行于根弦的弦的長度11111(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

16、.3125)*right*(x1,z1):(0.25,0)(x3,z3):(0.75,0)(x2,z2):(0.875,-0.625)(x4,z4):(1.375,-0.625)(x,z):(1.0625,-0.3125)(x1,z1):(0.875,0.625)(x3,z3):(1.375,0.625)(x2,z2):(1.5,1.25)(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)

17、:(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)(x,z):(2.3125,1.5625)*right*(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)(x1,z1):(2.125,1.875)(x3,z3):(2.625,1.875)(x2,z2):(2.75,2.5)(x4,z4):(3.25,2.5

18、)(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)*right*(x1,z1):(0.25,0)(x3,z3):(0.75,0)(x2,z2):(0.875,-0.625)(x4,z4):(1.375,-0.625)(x,z):(1.0625,-0.3125)(x1,z1):(0.875,0.625)(x3,z3):(1.375,0.625)(x2,z2):(1.5,1.25)(x4,z4):(2,1.25)(x,z):(1

19、.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)(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):(2,-1.25)(x4,z4):(2.625,-

20、1.875)(x,z):(2.3125,-1.5625)(x1,z1):(2.125,1.875)(x3,z3):(2.625,1.875)(x2,z2):(2.75,2.5)(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)*方程組系數(shù)矩陣*-1.13826-0.2946750.179738-0.03263340.0171196-0.009369

21、350.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.01711960.0617391-0.02463680.32177-0.0115021-1.13826-0.006009450.

22、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*線性方程組的右端項-1.18682-1.18682-1.18682-1.18682-1.18682-1.18682-1.

23、18682-1.18682*Gauss-seidel法解線性方程組迭代20步的結(jié)果(每個渦格的環(huán)量)*1.04267-1.31261.40375-1.489461.53951-1.57981.61008-1.632431.69757-2.02468-1.808621.799811.92227-1.80888-1.961492.004671.93371-2.10804-1.971311.8452.08496-1.84806-2.098852.101222.00797-2.12835-2.019311.857062.12727-1.85798-2.131442.126282.02866-2.13356-2.031641.860292.13859-1.86054-2.139722.132992.03405-2.13491-2.034821.861142.14156-1.86121-2.141862.134762.03545-2.13526-2.035651.861362.14234-1.86138-2.142422.135222.03581-2.13535-2.035861.861422.14254-1.86142-2.142562.135342.03591-2.13537-2.035921