涡格法代码及解释_物理_自然科学_专业资料

1、#include "iostream.h"#includc "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 x 1 ,z 1 ,x2,z2;左右自由涡的坐标double x3,z3,x4,z

2、4;3/4 弦线处的坐标double x,z;控制点的坐标,3/4弦线中点gird() x 1 =0.0f,x2=0.0f,z 1 =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 x 1,double zl,double x2,double z2 )double a,b,c,d,c;a=l/(x2-x)*(zl-z)-(x 1 -x)*(z2-z);b=(x2-x l)*(x i -x)+(z2-z l)*(zl -z)/sqrt(pow(x 1 x

3、),2)+pow(z 1 -z),2);c=(x2-x i)*(x2-x)+(z2-z 1 )*(z2z)/sqrt(pow(x2x),2)+pow(z2z),2);d=(l -(x 1 -x)/sqrt(pow(x 1 -x),2)+pow(z 1 -z),2)/(z 1 -z);c=(l -(x2-x)/sqrt(pow(x2-x),2)+pow(z2-z),2)/(z2-z);return (a*(b-c)+d-e)/4/pi;void gaussscidcl(int n,double *m,doubledouble *x,double *b)高斯塞得尔迭带法int t二0,i,j;迭代

4、次数while(t<20)/次数限制，精度要求，此处可修改，是迭带开关for(i=0;i<n;i+)mi = 0; for(j=0;j<n;j+)if(i 匸j)mi+=aij*x|j;xi = (bi - mi)/aii;/迭代 )cout«+t;for(i=0;i<n;i+)if(i%5=0) cout«cndl;) cout«n "«xij;cout«endl;void main()airfoil airfoil;int ng,nq,i,j,k,l,m,n,x,y;double y=0.0,m,a,ep=

5、 1 e-10,p=1.22505,cy=0.0; /p 为海平面空气密度 comvv “这是一个用涡格法计算机翼升力的程序!"«endl;cout«nih输入翼型个参数:展长l,根弦bg,前缘后掠角xo,丿匸缘后掠角xch«endl;whilc(l) cin»airfoil.l»airfoil.bg»airfoil.xo»airfoil.xc;if(airfoil.bg-airfoil.l*(tan(airfoil.xo*pi/180)+tan(airfoil.xc*pi/180)/2>0)cout

6、1;airfoil.l«h u«airfoil.bg«h h«airfoil.xo«u“vvairfoil.xcvv “"«endl;break;elsecout«*翼型的稍弦为0!请重新输入翼型数据“vvendl;cout«h请输入来流马赫数和攻角n«endl;cin»m»a;a=a*pi/180;cout«m«,t,«a«endl;coutvv”请输入根弦上的节点数，前缘上的节点数:” vvendl; cin»ng

7、7;nq;cout«ng«m n«nq«h h«endl;nq;ng;变成分多少块double *baseq=new doublefnq+l;double *bascb=ncw doublenq+l;double *result=new double2*nq*ng;double *b=new double2*nq*n創;double *ml=new double2*nq*ng;gird *girdleft,*girdright;/左半边机翼，右半边机翼girdlcft=ncw gird*ng;for(i=0;i<ng;i+)girdlef

8、ti=new girdfnq;girdright=new gird*ng;for(i=0;i<ng;i+)girdrightfi=new girdnq;double width二airfoil.l/nq/2;展长每个分块的长度前缘节点的x坐标cout«n前缘节点处的x坐标"«endl;for(i=0;i<nq+1 ;i+)baseqi=0+i*width*tan(airfoil.xo*pi/180); cout«baseqij«" "«endl;每一条平行于根弦的弦的长度cout«h每一条平行

9、于根弦的弦的长度n«endl;for(i=0;i<nq+l；i+) basebi=airfoil.bg-i*(tan(airfoil.xo*pi/180)+tan(airfoil.xc*pi/180)*width; cout«basebfi«" "«endl;for(i=0;i<ng;i+)for(j=0；j<nq;j+)girdleftij.xl=baseqj+basebj/4/ng+i*basebj/ng; girdrightij.xl=girdlefti|j.x 1;girdlefti j .x3=girdle

10、fti j .x 1 +baseb j/2/ng; girdrightij.x3=girdleftij.x3;girdlefti j.zl =0+j *width; girdrightli (j j .z 1 =-1 *girdlefti |j .z 1;girdlefti j .z3=girdlefti |j .z 1; girdright ij .z3=-1 *girdleftifj .z3; girdlcftij.z2=girdlcfti|j.zl+width; girdrightij.z2=-l*girdleftij.z2;girdlefti|j.z4=girdleftijljj.z2

11、;girdright i j j .z4=-1 *girdlefti j .z4;girdleftij.x2=baseqj+l+baseb j+l/4/ng+i*baseb j+11/ng; girdrightfi j .x2=girdlefti j .x2;girdleftij.x4=girdleftij.x2+basebj+l/2/ng;girdrighti j .x4=girdlefti j *x4; girdleftij.x=(girdleftij.x3+girdleftijj.x4)/2; girdrightfifjl.x=girdleftfij.x;girdleftij.z=(gi

12、rdleftij.z3+girdleftij.z4)/2; girdrighti j .z=-1 *girdlefti j .z;cout«h(x 1 ,z l):n«,("«girdleftij.x 1 «,«girdleftirj.z 1 «n),«n”;将坐标打cout«,(x2,z2):,«,(,，«girdleftij.x2«,；,«girdleftij.z2«,)n«endl; coutvv”(x3,z3):”vv“(”vvgirdl

13、eftijx3vv”,“vvgirdleftijz3vv”)“vv”“；cout«,(x4,z4):,«,，c，«girdleftlijljj.x4«,；,«girdleftlij|jj.z4«n)n«,*”;cout«,(x,z):,«uc，«girdleftij.x«n,n«girdleftij.z«,),«endl;cout«n(x 1 ,z 1): ,«n(h«girdrighti j .x 1 «u, m&

14、#171;girdrighti j .z 1“; 将坐标打岀打岀cout«,(x2,z2):,«,(，«girdrightij.x2«,m«girdrightijlj.z2«m),«endl; cout«,(x3,z3):,«',(h«girdrightlij|jj.x3«，v，«girdrightlijj.z3«,),«h“；cout«,(x4,z4):n«,("«girdrightij.x4«,

15、n«girdrightij.z4«,),«h“；cout«,(x,z):,«,(，«girdrighti|j.x«,«girdrighti|j.z«,),«endl;存储系数矩阵double *array;array=new double*l2*ng*nqj;for(i=0;i<2*ng*nq;i+)arrayi=ncw doubic2*ng*nq;k=i%nq;l=i/nq; for(j=0;j <n q*ng ;j+)m=j%nq;n=j/nq;x=2*i;y=2*j;array

16、 x y =vec(girdleftl k .x,girdleftl k .z,girdleftn m .x 1 ,girdleftn m .z 1,girdleft n m .x 2,girdleftn m .z2);array x y+1 =vec(girdleftl k .x,girdleftl k .z,girdrightn m .x 1 ,girdrightn m .zl ,girdright n m .x2,girdrightn m .z2);arrayx+ljy=vec(girdrightlllkj.x,girdrightlljkj.z,girdleftnjlm.xl,girdl

17、eftnjm.zl,girdleftn m .x2,girdleftn m .z2);array x+1 y+1 =vec(girdrightl k .x,girdrightl k ,z,girdrightn m .x 1 ,girdrightn m .z 1 ,girdr ightnm.x2,girdrightnm.z2);cout«n * * 方程组系数矢冃阵 *for(i=0;i<2*ng*nq;i+)for(j=0;j<2*ng*nq;j+4-)cout«arrayij«" ”；cout«endl;cou1*线性方程组的右端

18、项*冶*"«endl;for(i=0;i<2*ng*nq;i+)bil=-l*340*m*a;cout«bi«cndl;cout«h*gauss-seidel法解线性方程纟h.迭代20步的结果（每个涡格的环量)* * * *軻«cndl;for(i=0;i<2*ng*nq;i+)resultlij=o.o;gaussscidcl(2*nq*ng,m 1,array,result,b);for(i=0;i<ng*nq;i+)airfoil. y=airfoil.y+2*p*m*340*width*result2*i j

19、; airfoil.s=(baseb0+basebnq)*airfoil.l/2; airfoil.cy=2*airfoil.y/p/pow(m*340,2)/airfoil.s; cout«,y=n«airfoil.y«'t,«ncy=',«airfoil.cy«cndl;为了验证代码的正确性，此处的算例采用的是空气动力学一书中关于涡格法的一道算例，书中给出了算例的过程和解。丄寸.丄丄丄丄丄丄寸.丄丄亠丄寸.叶.叶.叶.吟.叶.吟.叶.叫.吟.彳.叶.吟.彳.叶.吟.叶.运行结果立立彳吟彳吟彳吟吟彳吟彳彳吟这是一个用

20、涡格法计算机翼升力的程序！请输入翼型个参数:展长l,根弦bg,前缘后掠介xo,后缘后掠角xc5 145-455145-45请输入来流马赫数和攻角0.20.20.0174533请输入根弦上的节点数，前缘上的节点数:2 525询缘节点处的x人怡标00.6251.251.8752.5每一条平行于根弦的弦的长度1(x 1 ,zl):(0.25,0)(x2,z2):(0.875,0.625)(x3,z3):(0.75,0)(x4,z4):(l .375,0.625)(x,z):(l .0625,0.3125)rj r| rjw rj r|rj* rjrj* rj rjw rj rjw rjw rj rj

21、w r| rjwrj* rj rj* rj* rj* rj rjw rj rjw rj(x 1 ,zl ):(0.25,0)(x2,z2):(0.875,-0.625)(x3,z3):(0.75,0)(x4,z4):(l .375,-0.625)(x,z):( 1.0625,-0.3125)y. 上 7厶 y. % 上土 y. 7厶上土 y. 上土 y. 丄 卜土 上 7厶上土(“ 土 土 上土 -if 上 *1* y. 土 *1* 土 上土 *1* 土土<twi(x 1 ,z 1 ):(0.875,0.625)(x2,z2):( 1.5,1.25)(x3,z3):(l .375,0.6

22、25)(x4,z4):(2,1.25)(x,z):(l .6875,0.9375)* * right * * *(x 1 ,z 1 ):(0.875,-0625)(x2,z2):(1.5,-l .25)(x3,z3):( 1.375,-0.625)(x4,z4):(2,1.25)(x,z):( 1.6875,-0.9375)tw tw »tw »tw 丫 *tw »tw *. rt» *t» »tw rtw rt» *tw »tw rtw rt» <tw rt» rt« rt&#

23、187; tw i 盲 #* rtw »tw rtw *tw »tw <t» rtw <!*t» rtw rtw »tw 打.: rtw(x 1 ,z 1):(1.5,1.25)(x2,z2):(2 j 25,1-875)(x3,z3):(2,1.25)(x4,z4):(2.625,1.875)(x,z):(2.3125,1.5625)丄“ 士 丄“ *1* 士士士 *1*- 丄“ 士 丄“ 士 丄“ 士 土 士 土 丄e 丄“卄"卄“代“代“卄“"“门ght“"不"""&

24、quot;"""水“(x 1 ,z 1):(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)tw丫rtw rt» rtw <tw rt« tw i 盲 #* rtwrtwrtwrtw rtw 打. rtw(x 1 ,zl ):(225,1.875)(x2,z2):(2.75,2.5)(x3,z3):(2.625,1.875)(x4,z4):(3.25,2.5)(x,z):(29375,2 1875)丄“

25、j士 丄“ *1* 士士士 *1* *1 1 - 丄“ 士 丄“ 士 丄“ 士 土 士 土 丄e *! 丄“卄"卄“代“代“卄“"“门ght“"不"""""""水“(x 1 ,z 1 ):(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

26、.009369350.00600848-0.004230970.2946751.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.01868781.13826().00780396-0.1797380.00398332-0.01711960.0617391 -0.02463680.32177-0.0115021-1.13826-0

27、.006009450.1797380.003467210.0246368 -0.06173910.0115021-0.321770-006009451.138260.00346721-0-1797380.0259969 -0.0136999 0.0617391-0.007693410.32177-0.00460806-1.13826-0.002921990.0136999 -0.0259969 0.00769341-0.06173910.002921991.138260.00460806-0.32177*线件方程纟甘的右端项*-1.18682-1.18682-1.18682-1.18682-1

28、.18682-1.18682-1.18682-1.18682*gaussseidel法解线性方程纽迭代20步的结果(每个涡格的环量)林*l042674.31261.40375-1.489461.53951-1.57981.61008-1.632431.69757-2.024681.93371-2.108042.00797-2.128352.02866-2.133562.03405-2.134912.03545-2.135262.03581-2.135352.03591-2.135372.03593-2.135382.03594-2.135382.03594-2.135382.03594-2.135382.03594-2.135382.03594-2.13538-1.808621.799811.92227-1.80888-1.961492.00467-1.9713