基于simple算法的流场模拟计算.docx
《基于simple算法的流场模拟计算.docx》由会员分享，可在线阅读，更多相关《基于simple算法的流场模拟计算.docx（8页珍藏版）》请在冰豆网上搜索。
基于simple算法的流场模拟计算
基于simple算法的流场模拟计算
1、问题描述
图1为20?
的水在长度为150mm，宽度为10mm的管道中流动，流入管道
速度为0.2m/s，流出管道背压为1atm，基于simple算法对整个流场进行计算，
计算管长100mm处流速并与fluent计算结果对比。
图1流动系统
流动计算:
dmm,50
3,,998.21/;kgm
6,，,1001.6110;Pas,
3ud998.210.21010，，，,Re2300,,,,61001.6110，,
该流动问题为二维定常无内热源不可压缩层流流动。
2、控制方程离散
在二维直角坐标系中，对流—扩散方程的通用形式为:
,,，，,,,,,,,,,,,,uS，，，，，，,,，,，,,,,,,,,,,,,txyxxyy,,,,,,,,,,,
图2网格编号
针对本问题，其连续性方程和动量方程为:
,,,u,v（）（），,0,,x,y,
,,,,,uu,uv,,u,,u,p（）（）,,,,,,,，,，,,,,,,,x,y,x,x,y,y,x,,,,,
,,,,,vu,vv,,v,,v,p（）（）,,,,,，,，,,,,,,,,x,y,x,x,y,y,y,,,,,,
交错网格下动量方程的离散:
nenenene,,,,,pu,,（）（）,,uudxdyvudxdydxdyudxdy，,,，,,,,,,,,,,,,,,,,,,,xyxxx,,swswswsw
nene,,,,u，，udxdySdxdyu,,,,,,,,yy,,swsw
上式积分得到:
1FuFuyFuFuxppyDuuDuu,,，,,,,,，,,）-）（）（）（（，，eewwnnssweeEPwPW2
，,,，,,DuuDuuSSuxy）-）+（（，，nNPsPSCPP
对流项采用一阶迎风格式，扩散项采用中心差分，
,,,yywe（）,（）,FuyDFuyD,,,,,,,,,wwweeexx,,WPPE
,xy,,ns（）,（）,FvxDFuyD,,,,,,,,,nnnsssyx,,PNSP
1,auauauauauppyb,，，，，,,，（）PPWWEESSNNWE,2,aDFaD,，,,WwwEe,
11,aDFaDF,，,,,,SssNnn22,
aaaaaFFFF,，，，，,，,（）（）,PWESNewns,1,bPP,,,y（）WE,2,
开始
变量的初场是否需要视情况而定）假设一个速度初场（其它,
假设一个压力场，即给定压力猜测值
根据当前的已知量，计算动量离散方程等方程中的系数和常数项
0赋值:
pp,
步骤1:
依次求解动量离散方程
****auauAppb,，,，（）,eenbnbePE****avavAppb,，,，（）,nnnbnbnPN
赋值:
,0uv,步骤2:
根据速度求解压力修正方程pp,
0''''''uu,apapapapapb,，，，，PPWWEESSNN0'****,,,buuyvvx,,,，,,[（）（）][（）（）],,,,wesn
步骤3:
对压力和速度进行修正
*'ppp,，**'''（）uuuudpp,，,，,eeeeePE
**'''（）vvvvdpp,，,，,nnnnnPN
aab,,,，,步骤4:
求解其他变量的离散方程（视需要进行）,ppnbnb
否
收敛否,
是
结束
TDMA算法
,C11
，,,DC,,,,,2122232,，,,DC,,,,,3233343
，,,DC4344454,,,,,
，,,,,,,,DCnnnnnnn,，11
C,nn，，11
在上式中，假定和是边界上的值，为己知。
上式中任一方程都可写成:
,1n，1
，,,,,,,,DCjjjjjjj,，11除第一及最后一个方程外，其余方程可写为:
,C222,，，,,,231DDD222
C,,333,，，,,,342DDD333
C444,,,，，453,,,DDD444
,Cnnn,，，,,,nnn，,11DDDnnn
这些方程可通过消元和回代两个过程求解。
,,,,C,,22++C,,,,313,,,,DD,22,,3,,,,,，,,34,,,,22,,DD,,,,3333,,,,DD,,22,,现引入记号:
,C'222AC,,,+,221DDD222
',,,,,,CC+3323则，,，,,,,,,34DADA,,,,332332,,,,
',,CC+'3323,,AC,33,,DADA,,332332
',,,，AC即:
，jjjj1
',,CC+jjjj,1',,AC,jj,,DADA,,jjjjjj,,11
在边界点，j=1与j=n+1，
'AC,,0,,111
'AC,,0,,nnn，，，111
为了求解方程组，首先要对方程组按的形式编排，并明确其中的系数和。
从D,,Cjjjj
'j=2起，计算出和，直到j=n。
由于在边界位置（n+1）的数值是已知的，因此，C,Ajj
可连续计算出。
j
3、计算结果对比
本文基于simple算法用编写MATLAB程序对整个流场进行计算，另外借助fluent计算流体力学软件数值模拟，通过对比分析计算结果，得出整个流场的流速矢量图和速度云图。
0.30
数值计算0.25matlab计算
0.20
0.15
u
0.10
0.05
0.000.00000.00250.00500.00750.01000.0125
y
上图为两种计算方法的计算结果，发现两者速度分布趋势一致，编程计算结果数值偏大，原因在于计算采用一阶迎风格式，节点数值趋近于其上游节点值。
下图为速度云图，入口段效应的影响导致流动未达到充分发展。
MATLAB程序:
%%%%%%%%%%%%%%%%%%%%%%%%dx=5e-4;
dy=5e-4;
den=998;
dyna=1001.6e-6;
a=0.01/dy+3;%x节点%
b=0.15/dx+1;%y节点%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%边界条件的设置、初始值的设置%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
T=0.2*ones（a,b）;
U1_old=zeros（a,b）;
P1_old=zeros（a,b）;%初值%U1_old（3:
a-2,:
）=T（3:
a-2,:
）;%加速收敛，边界条件%
uw_sum（1,b）=0;
ue_sum（1,b）=0;
de1_sum（1,b）=0;
m=1;
A（1:
a-4,2:
b-1）=0.2;%初值为1%e22=1;
whilee22>0.01%整个流场迭代%
ifm>4
break
end
m=m+1;
forj=2:
1:
b-1;
ifj>2
break
end
n=1;%迭代次数%
e=1;%初次迭代%
B=[0;0;0.0823455;0.149211;0.196652;0.224194;0.23674;0.241087;0.242138;0.242214;0.242125;0.242083;0.242125;0.242214;0.242138;0.241087;0.23674;0.224195;0.196652;0.149212;0.0823459;0;0];
U1_old（3:
a-2,1）=B（3:
a-2,1）;
U1_old（3:
a-2,2）=B（3:
a-2,1）;
whilee>0.001%迭代，算两列，其他为已知%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%
%对（i,j）点求uxin%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%
fori=3:
1:
a-2
Fe1（i,j）=den*（U1_old（i,j）+U1_old（i,j+1））*dy/2;%对流强度%
Fw1（i,j）=den*（U1_old（i,j-1）+U1_old（i,j））*dy/2;
ap_u1（i,j）=max（Fe1（i,j）,0）+max（-Fw1（i,j）,0）+4*dyna;%离散方程对应系数%
ae_u1（i,j）=dyna+max（-Fe1（i,j）,0）;
aw_u1（i,j）=dyna+max（Fw1（i,j）,0）;
an_u1（i,j）=dyna;
as_u1（i,j）=dyna;
A1_u（i,i-1）=-an_u1（i,j）;
A1_u（i,i）=ap_u1（i,j）;
A1_u（i,i+1）=-as_u1（i,j）;
b1_u（i,j）=aw_u1（i,j）*U1_old（i,j-1）+ae_u1（i,j）*U1_old（i,j+1）...
+（P1_old（i,j）-P1_old（i,j+1））*dy/2;%（i,j+1）点的源项%
end
A1_u（1,1）=1;%附加点%
A1_u（2,2）=1;%边界点%
A1_u（a-1,a-1）=1;
A1_u（a,a）=1;
b1_u（1,j）=0;%附加点%
b1_u（2,j）=0;%边界点%
b1_u（a-1,j）=0;
b1_u（a,j）=0;
U1_new（:
j）=inv（A1_u）*b1_u（:
j）;%得到（i,j+1）点的速度计算值%%对（i，j）点求p%
uw_sum（1,j）=0;
ue_sum（1,j）=0;
de1_sum（1,j）=dy/ap_u1（2,j）;
fori=3:
1:
a-2
uw_sum（1,j）=uw_sum（1,j）+U1_old（i,j-1）;
ue_sum（1,j）=ue_sum（1,j）+U1_new（i,j）;
de1_sum（1,j）=de1_sum（1,j）+dy/ap_u1（i,j）;
end
P1_fix（1,j）=（uw_sum（1,j）-ue_sum（1,j））/de1_sum（1,j）;%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%
%对（i，j）点进行速度修正%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%
fori=3:
1:
a-2
U1_fix（i,j）=（dy/ap_u1（i,j））*P1_fix（1,j）;%速度修正%
end
U1_old（3:
a-2,j）=U1_new（3:
a-2,j）+U1_fix（3:
a-2,j）;
P1_old（2:
a-1,j）=P1_old（2:
a-1,j）+0.4*P1_fix（1,j）;
ee=uw_sum（1,j）-ue_sum（1,j）;
e=max（max（abs（ee）））;%判断的是否合理,%
n=n+1;
N（1,j）=n;
ifn>1000
break
end
A（（a-4）*m+1:
（a-4）*（m+1）,2:
b-1）=U1_old（3:
a-2,2:
b-1）;%赋值，即m从1开始%
end
end
E（1:
a-4,2:
j）=A（（a-4）*m+1:
（a-4）*（m+1）,2:
j）-A（（a-4）*（m-1）+1:
（a-4）*m
2:
j）;
e22=max（max（abs（E）））;
end%整个流场迭代结束%
y=0:
1:
10;%结果的输出%
u=U1_old（2:
2:
a-1）;
u
plot（y,u）
升级会员 