渦格法代碼及解釋

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/網(wǎng)格信息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;/迭代次數(shù)while(t<20)/次數(shù)限制,

4、精度要求,此處可修改,是迭帶開關(guān) 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<<"請輸入翼型個參數(shù):展長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!請重新輸入翼型數(shù)據(jù)"<<endl;cout<<"請輸入來流馬赫數(shù)和攻角"<<endl;cin>

7、;>M>>a;a=a*PI/180;cout<<M<<'t'<<a<<endl;cout<<"請輸入根弦上的節(jié)點數(shù),前緣上的節(jié)點數(shù):"<<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、塊的長度/前緣節(jié)點的x坐標cout<<"前緣節(jié)點處的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;/存儲系數(shù)矩陣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<<"*方程組系數(shù)矩陣*"<<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步的結(jié)果(每個渦格的環(huán)量)*"<<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;為了驗證代碼的正確性，此處的算例采用的是空氣動力學(xué)一書中關(guān)于渦格法的一道算例，書中給出了算例的過程和解。*運行結(jié)果*這是一個用渦格法計算機翼升力的程序!請輸入翼型個參數(shù):展長L, 根弦Bg,前緣后掠角Xo,后緣后掠角Xc5145-455 1 45 -45請輸入來流馬赫數(shù)和攻角0.210.2 0.0174533請輸入根弦上的節(jié)點數(shù),前緣上的節(jié)點數(shù):252 5前緣節(jié)點處

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)*方程組系數(shù)矩陣*-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步的結(jié)果(每個渦格的環(huán)量)*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