涡格法代码及解释.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
可见程序正确
欢迎您的下载，
资料仅供参考!
致力为企业和个人提供合同协议，策划案计划书，学习资料等等
打造全网一站式需求
升级会员 