涡格法代码及解释.docx
- 文档编号:30760262
- 上传时间:2023-08-20
- 格式:DOCX
- 页数:15
- 大小:35.95KB
涡格法代码及解释.docx
《涡格法代码及解释.docx》由会员分享,可在线阅读,更多相关《涡格法代码及解释.docx(15页珍藏版)》请在冰豆网上搜索。
涡格法代码及解释
#include"iostream.h"
#include"stdio.h"
#include"math.h"
#definePI3.1415926
classAIRFOIL//用来存放翼型的信息
{
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//网格信息
{
public:
doublex1,z1,x2,z2;//左右自由涡的坐标
doublex3,z3,x4,z4;//3/4弦线处的坐标
doublex,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;}
};
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),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*b)//高斯--塞得尔迭带法
{
intt=0,i,j;//迭代次数while(t<20)//次数限制,精度要求,此处可修改,是迭带开关{
for(i=0;i { M[i]=0; for(j=0;j { if(i! =j) {M[i]+=a[i][j]*x[j]; } } x[i]=(b[i]-M[i])/a[i][i];//迭代 } cout<<++t; for(i=0;i {if(i%5==0){cout< } cout< } } 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<<"这是一个用涡格法计算机翼升力的程序! "< cout<<"请输入翼型个参数: 展长L,根弦Bg,前缘后掠角Xo,后缘后掠角Xc"< (1) { 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< break; } else {cout<<"翼型的稍弦为0! 请重新输入翼型数据"< } } cout<<"请输入来流马赫数和攻角"< a=a*PI/180; cout< cout<<"请输入根弦上的节点数,前缘上的节点数: "< cout< Nq--;Ng--;//变成分多少块 double*baseq=newdouble[Nq+1]; double*baseB=newdouble[Nq+1]; double*result=newdouble[2*Nq*Ng]; double*b=newdouble[2*Nq*Ng]; double*M1=newdouble[2*Nq*Ng]; GIRD**girdleft,**girdright;//左半边机翼,右半边机翼girdleft=newGIRD*[Ng]; for(i=0;i {girdleft[i]=newGIRD[Nq]; } girdright=newGIRD*[Ng]; for(i=0;i {girdright[i]=newGIRD[Nq]; } doublewidth=airfoil.L/Nq/2;//展长每个分块的长度 //前缘节点的x坐标 cout<<"前缘节点处的x坐标"< {baseq[i]=0+i*width*tan(airfoil.Xo*PI/180);cout< } //每一条平行于根弦的弦的长度 cout<<"每一条平行于根弦的弦的长度"< for(i=0;i {baseB[i]=airfoil.Bg-i*(tan(airfoil.Xo*PI/180)+tan(airfoil.Xc*PI/180))*width;cout< } for(i=0;i {for(j=0;j girdleft[i][j].x1=baseq[j]+baseB[j]/4/Ng+i*baseB[j]/Ng; girdright[i][j].x1=girdleft[i][j].x1; girdleft[i][j].x3=girdleft[i][j].x1+baseB[j]/2/Ng; girdright[i][j].x3=girdleft[i][j].x3; girdleft[i][j].z1=0+j*width; girdright[i][j].z1=-1*girdleft[i][j].z1; girdleft[i][j].z3=girdleft[i][j].z1; girdright[i][j].z3=-1*girdleft[i][j].z3; girdleft[i][j].z2=girdleft[i][j].z1+width; girdright[i][j].z2=-1*girdleft[i][j].z2; girdleft[i][j].z4=girdleft[i][j].z2; girdright[i][j].z4=-1*girdleft[i][j].z4;girdleft[i][j].x2=baseq[j+1]+baseB[j+1]/4/Ng+i*baseB[j+1]/Ng;girdright[i][j].x2=girdleft[i][j].x2; girdleft[i][j].x4=girdleft[i][j].x2+baseB[j+1]/2/Ng; girdright[i][j].x4=girdleft[i][j].x4;girdleft[i][j].x=(girdleft[i][j].x3+girdleft[i][j].x4)/2; girdright[i][j].x=girdleft[i][j].x;girdleft[i][j].z=(girdleft[i][j].z3+girdleft[i][j].z4)/2; girdright[i][j].z=-1*girdleft[i][j].z;cout<<"************************left**********************"< "<<"("< ";// 将坐标打出 cout<<"(x2,z2): "<<"("< cout<<"(x3,z3): "<<"("< J J cout<<"(x4,z4): "<<"("< cout<<"(x,z): "<<"("< cout<<"************************right**********************"< cout<<"(x1,z1): "<<"("< ";// 将坐标打出 cout<<"(x2,z2): "<<"("< cout<<"(x3,z3): "<<"("< J J cout<<"(x4,z4): "<<"("< cout<<"(x,z): "<<"("< } } //存储系数矩阵 double**array; array=newdouble*[2*Ng*Nq];for(i=0;i<2*Ng*Nq;i++) { array[i]=newdouble[2*Ng*Nq]; } for(i=0;i { k=i%Nq;l=i/Nq; for(j=0;j { m=j%Nq;n=j/Nq;x=2*i;y=2*j; array[x][y]=vec(girdleft[l][k].x,girdleft[l][k].z,girdleft[n][m].x1,girdleft[n][m].z1,girdleft[n][m].x2,girdleft[n][m].z2); array[x][y+1]=vec(girdleft[l][k].x,girdleft[l][k].z,girdright[n][m].x1,girdright[n][m].z1,girdright[n][m].x2,girdright[n][m].z2); array[x+1][y]=vec(girdright[l][k].x,girdright[l][k].z,girdleft[n][m].x1,girdleft[n][m].z1,girdleft[n][m].x2,girdleft[n][m].z2); array[x+1][y+1]=vec(girdright[l][k].x,girdright[l][k].z,girdright[n][m].x1,gir dright[n][m].z1,girdright[n][m].x2,girdright[n][m].z2); for(i=0;i<2*Ng*Nq;i++) {for(j=0;j<2*Ng*Nq;j++){ cout< } cout<<"***************Gauss-seidel 法解线性方程组迭代 20步的结果(每个涡格的环 "< for(i=0;i<2*Ng*Nq;i++) result[i]=0.0; } Gaussseidel(2*Nq*Ng,M1,array,result,b); for(i=0;i airfoil.Y=airfoil.Y+2*p*M*340*width*result[2*i];airfoil.S=(baseB[0]+baseB[Nq])*airfoil.L/2;airfoil.Cy=2*airfoil.Y/p/pow(M*340,2)/airfoil.S;cout<<"Y="< } 为了验证代码的正确性,此处的算例采用的是《空气动力学》一书中关于涡 格法的一道算例,书中给出了算例的过程和解 这是一个用涡格法计算机翼升力的程序! 请输入翼型个参数: 展长L,根弦Bg,前缘后掠角Xo,后缘后掠角Xc 45-45 5145-45请输入来流马赫数和攻角 0.2 1 0.20.0174533 请输入根弦上的节点数,前缘上的节点数 2 5 25 前缘节点处的x坐标 0 0.625 1.25 1.875 2.5每一条平行于根弦的弦的长度 1 ************************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): (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): (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) ************************AC: ********************** ****************方程组系数矩阵*************** -1.13826-0.294675 0.179738 -0.0326334 0.0171196 -0.00936935 -0.00423097 0.2946751.13826 0.0326334-0.179738 0.00936935 -0.0171196 -0.00600848 0.32177-0.0575242 -1.13826 -0.0186878 0.179738 -0.00780396 -0.00398332 0.0575242-0.32177 0.0186878 1.13826 0.00780396 -0.179738 -0.0171196 0.0617391-0.0246368 0.32177 -0.0115021 -1.13826 -0.00600945 -0.00346721 0.0246368-0.0617391 0.0115021 -0.32177 0.00600945 1.13826 -0.179738 0.0259969-0.0136999 0.0617391 -0.00769341 0.32177 -0.00460806 -0.00292199 0.0136999-0.0259969 0.00769341 -0.0617391 0.00460806 -0.32177 (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) 0.00600848 0.00423097 0.0171196 0.00398332 0.179738 0.00346721 -1.13826 0.00292199 1.13826 *************** 线性方程组的右端项************* -1.18682 -1.18682 -1.18682 -1.18682 -1.18682 -1.18682 -1.18682 -1.18682 ***************Gauss-seidel 法解线性方程组迭代20步的结果(每个涡格的环量)********* 1.04267-1.31261.40375-1.489461.53951 -1.57981.61008-1.63243 1.69757-1.80862 1.92227-1.96149 2.00467 -2.024681.79981 1.93371-1.97131 -1.80888 2.08496-2.09885 -2.108041.845-1.84806 2.00797-2.01931 -2.128351.85706 2.02866-2.03164 -2.133561.86029 2.03405-2.03482 -2.134911.86114 2.03545-2.03565 -2.135261.86136 2.03581-2.03586 -2.135351.86142 2.03591-2.03592 -2.135371.86143 10 2.03593-2.03593 -2.135381.86144 11 2.03594-2.03594 2.10122 2.12727-2.13144 -1.85798 2.13859-2.13972 -1.86054 2.14156-2.14186 -1.86121 2.14234-2.14242 -1.86138 2.14254-2.14256 -1.86142 2.12628 2.13299 2.13476 2.13522 2.13534 2.14259-2.14262.13537 -1.86143 2.14261-2.142612.13538 -1.86144 2.14261-2.142612.13538 -2.135381.86144-1.86144 12 2.03594-2.035942.14261-2.142612.13538 -2.135381.86144-1.86144 13 2.03594-2.035942.14261-2.142612.13538 -2.135381.86144-1.86144 14 2.03594-2.035942.14261-2.142612.13538 -2.135381.86144-1.86144 15 2.03594-2.035942.14261-2.142612.13538 -2.135381.86144-1.86144 16 2.03594-2.035942.14261-2.142612.13538 -2.135381.86144-1.86144 17 2.03594-2.035942.14261-2.142612.13538 -2.135381.86144-1.86144 18 2.03594-2.035942.14261-2.142612.13538 -2.135381.86144-1.86144 19 2.03594-2.035942.14261-2.142612.13538 -2.135381.86144-1.86144 20 2.03594-2.035942.14261-2.142612.13538 -2.135381.86144-1.86144 Y=851.296Cy=0.0601131 在同样的条件下,《空气动力学》书中给出的结果为: Y=848.554Cy=0.059919 可见程序正确 欢迎您的下载, 资料仅供参考! 致力为企业和个人提供合同协议,策划案计划书,学习资料等等 打造全网一站式需求
- 配套讲稿:
如PPT文件的首页显示word图标,表示该PPT已包含配套word讲稿。双击word图标可打开word文档。
- 特殊限制:
部分文档作品中含有的国旗、国徽等图片,仅作为作品整体效果示例展示,禁止商用。设计者仅对作品中独创性部分享有著作权。
- 关 键 词:
- 涡格法 代码 解释
![提示](https://static.bdocx.com/images/bang_tan.gif)