基于simple算法的流场模拟计算.docx
《基于simple算法的流场模拟计算.docx》由会员分享，可在线阅读，更多相关《基于simple算法的流场模拟计算.docx（12页珍藏版）》请在冰点文库上搜索。
基于simple算法的流场模拟计算
基于sｉmｐle算法的流场模拟计算
————————————————————————————————作者：
———————————————————————————————— 日期：
ﻩ
1、问题描述
图1为20℃的水在长度为１50ｍm,宽度为10mm的管道中流动,流入管道速度为0.2ｍ/s，流出管道背压为1atm，基于simple算法对整个流场进行计算，计算管长10０mm处流速并与fluenｔ计算结果对比。
图1流动系统
流动计算:
该流动问题为二维定常无内热源不可压缩层流流动。
2、控制方程离散
在二维直角坐标系中，对流—扩散方程的通用形式为：
图２网格编号
针对本问题，其连续性方程和动量方程为:
交错网格下动量方程的离散：
上式积分得到：
对流项采用一阶迎风格式，扩散项采用中心差分,
否
ﻩ
是
ﻩ
TDＭA算法
在上式中,假定
和
是边界上的值，为己知。
上式中任一方程都可写成:
除第一及最后一个方程外，其余方程可写为:
这些方程可通过消元和回代两个过程求解。
现引入记号：
则，
即：
在边界点，ｊ=1与j=n+１,
为了求解方程组,首先要对方程组按的形式编排，并明确其中的系数
和
。
从ｊ=２起，计算出
和
，直到j=n。
由于在边界位置（n+１）
的数值是已知的,因此,可连续计算出
。
3、计算结果对比
本文基于sｉmpｌe算法用编写MATLAB程序对整个流场进行计算，另外借助fluｅnt计算流体力学软件数值模拟，通过对比分析计算结果,得出整个流场的流速矢量图和速度云图。
上图为两种计算方法的计算结果,发现两者速度分布趋势一致,编程计算结果数值偏大，原因在于计算采用一阶迎风格式，节点数值趋近于其上游节点值。
下图为速度云图,入口段效应的影响导致流动未达到充分发展。
MＡTLAＢ程序：
%％%%%%%%%%%%％%％%%％％%％%％%
dｘ=5e－4;
ｄｙ＝5ｅ-４;
den=998；
dyｎａ=100１.6e-6；
ａ=0．０1/dy+３；%x节点％
b=0．１5/dx+１;%y节点%
%％%%％%%%%%%%%%％%%%%%％%％%％%%％％%%%%％％%%%%％％%%％%％％%%%%％%%%%％%％%％%%％%
%边界条件的设置、初始值的设置%
%%％%%%%%%％%%％％％%％%%％%%％%%%%%%％%%%%%%％％%%%％%％％%%%%%%%%％%%％%%％%％%%%
T=０.２*oｎｅs（a,ｂ）;
U1_old=zeｒoｓ（a,b）；
P1_old=zeros（ａ,ｂ）；％初值%
U1_old（３:
ａ-2,:
）=T（3:
ａ-2,:
）；%加速收敛，边界条件%
uw_sum（1,b）＝0;
ue_sum（1,b）＝0；
de1_sum（1,b）=０；
m＝1;
A（1:
a-４,2:
b－1）=０．2;%初值为1%
e22＝１;
whilee22>0．01％整个流场迭代％
if m＞4
break
ｅnd
m=m+1;
ｆorj=2:
1：
b-1；
ｉfj＞２
break
ｅnd
n=1;%迭代次数%
e=1；%初次迭代%
Ｂ＝［０；0;0.0８２345５;０.149211;0．196６５2;0.224194;０.2３67４；
0.241087;０.242138;0.24２21４；0.２４2125;0.24２0８3；
0.242125;0．2422１4;0．242１38；0．２4１087;0.２3674；
0．22４195;0.1９６６５２;0.149212;0.082３459;0；0]；
U１_olｄ（３:
a-2，1）=B（3:
a－2,1）;
U1_old（3:
a-2,２）=B（3:
a－２,1）;
whilee>0．００1％迭代,算两列,其他为已知%
％％%%%％%%％%%%%%%%%％%％％%%%％％%％%%%%%%%%%%%%％％%%%%%%%％%%％%%％%%％%%%%％%%%%%％%%％%%
%对（i,j）点求uxin%
%％％%％%%%％%％％％%%%％%%％%%%％%%%%%%%％％％％%%%%%％%％%%%%%%％%%％%%%%％%%%％%%％％%％%%％％%%%
fori=３:
1:
a-2
Fe1（i,j）=ｄen＊（U1_old（i,j）＋U1＿olｄ（i,ｊ+1））*dy／2;％对流强度％
Ｆw1（i,j）=ｄeｎ＊（Ｕ1_old（i,j－1）+U１_old（i，j））＊dｙ/2;
ap＿ｕ1（ｉ，j）=ｍax（Fe1（i，j）,0）+max（－Fw１（i,j），０）+4*dynａ;%离散方程对应系数%
ae_u1（i，ｊ）=dyna+mａｘ（-Ｆe1（i,j）,0）;
aw＿u1（i，j）＝dyna+ｍａx（Ｆｗ１（i,ｊ）,0）;
ａn_ｕ1（i，j）=ｄyna;
as_u1（i,j）=dyｎa;
A1_u（i,i-1）＝-ａn＿u1（i,j）;
A１_u（ｉ,i）=ａp＿ｕ1（i,ｊ）;
A1_ｕ（i，ｉ+1）＝－ａs_u1（i，ｊ）；
b1_u（i,j）=aw＿u1（i,j）*U１_ｏｌｄ（i,j-1）+ａe_ｕ１（i,j）*U1_ｏｌd（i,j+１）...
+（P１＿ｏld（i,j）-P１＿oｌd（ｉ,j＋1））*dy/2;%（i，j+1）点的源项%
end
A1_ｕ（1,1）=1;％附加点%
A1_u（2,2）=1;％边界点％
A1＿ｕ（a－1，ａ-1）=1;
A1＿ｕ（a，a）＝１;
b1_ｕ（1，j）＝0;％附加点%
b1_u（2,ｊ）=0;%边界点%
ｂ１_u（a-1，j）=0；
b1＿u（a,j）＝0;
Ｕ1_new（:
j）＝inv（A１_u）*b１＿u（:
ｊ）；%得到（i，j＋1）点的速度计算值%
％对（i,j）点求p%
uw＿sum（1,j）＝0；
uｅ_sｕm（１,j）=0；
ｄe1_sｕｍ（１,j）=dｙ／ap_u1（2，j）;
for ｉ＝3：
1:
a-2
uw_ｓum（1,ｊ）＝uｗ_ｓum（1,j）+U1＿olｄ（i,j-１）;
uｅ＿sｕｍ（1,ｊ）=ue＿suｍ（１,j）＋U1_nｅw（i,j）;
de1_sum（1,j）=ｄe1＿sum（１,j）+dｙ／ａp_u1（i,j）;
ｅnd
P1＿fiｘ（１,j）=（uw_sum（1,j）-ｕe＿ｓum（1,j））／ｄe1_suｍ（1,j）;
％%%%%%%％%%%％%％%%%％%%%%％%％％%%%％％%％%%%%％%%%%%％%%%%%％％%%％％%%%%%%%%％%%％%%%%%%%%
%对（i,ｊ）点进行速度修正%
%%%%%％%%%%%%%%%%％%％%%%％％%%％%%%%%％％％%%％%%%%%%%%%%%%%%%%%%%%％％%％％%%%%%%%％%%％％
ｆori=3:
1:
a-2
U1＿ｆix（i,j）=（dy/ap_u1（i,j））＊P1_fix（1,ｊ）;％速度修正%
eｎｄ
U1_oｌd（３：
a-２,j）=Ｕ1_nｅw（3:
a-2,j）+U1_ｆｉx（３:
a-２,ｊ）;
Ｐ1_ｏｌd（2:
a-1,j）=P1＿oｌd（2:
a－１,j）+0．4＊P1_fix（1,j）；
ee＝uｗ_suｍ（１,j）-ue_ｓum（1,j）;
e=max（mａx（abｓ（ee）））;％判断的是否合理？
％
ｎ=n＋１;
Ｎ（１,ｊ）＝n;
ｉfn>1０00
breａk
enｄ
A（（a-4）*ｍ+1:
（a-4）＊（m+1），2:
b-1）=Ｕ１＿ｏlｄ（３:
a-2,2：
b-1）;%赋值,即ｍ从１开始%
eｎd
end
E（１:
ａ-4,2:
j）=A（（a-4）*ｍ+1:
（ａ－4）*（ｍ+1）,2:
j）-Ａ（（a－4）*（m-1）+１:
（a-4）＊m,2:
j）;
e22＝maｘ（maｘ（abs（E）））;
eｎd%整个流场迭代结束%
y=0:
1:
１０；%结果的输出%
u＝Ｕ1_old（2:
2：
a-1）;
u
ｐlｏt（y,u）
升级会员 