涡格法代码及解释.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<<"请输入来流马赫数和攻角"<>M>>a;
a=a*PI/180;
cout<cout<<"请输入根弦上的节点数,前缘上的节点数:
"<>Ng>>Nq;
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;jgirdleft[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;iairfoil.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
可见程序正确
欢迎您的下载,
资料仅供参考!
致力为企业和个人提供合同协议,策划案计划书,学习资料等等
打造全网一站式需求