基于simple算法的流场模拟计算

1、1、问题描述图1为20C的水在长度为150mm,宽度为10mm的管道中流动，流入管道速度为0.2m/s,流出管道背压为1atm，基于simple算法对整个流场进行计算，计算管长100mm处流速并与fluent计算结果对比。+±frayI流动计算：d=50mmp=998.21kg/m3；卩=1001.61x10-6Pa-s;pud998.21x0.2x10x10-3Re=<2300卩1001.61x10-6该流动问题为二维定常无内热源不可压缩层流流动。2、控制方程离散在二维直角坐标系中，对流一扩散方程的通用形式为:(pu9)+(pU9)=2axayaxv9ax针对本问题，其连续性

2、方程和动量方程为:d(pu)+d(pv)0dxdyd(puu)d(puv)d(du)v+=IH+dxdydxVdx丿"d(pvu)d(pvv)d(dv)dxdydxV"du、dy丿dv)<Hdpdxdyddx丿dy1dy丿dy交错网格下动量方程的离散：dpff匚(pu)udxdy+fff(pv)udxdy二-ff擎dxdy+邛ddxsw+嗚swdy'sw'du'Idy丿dxswdxdy+ffSdxdyusw'du、udxvdx丿swdxdy上式积分得到：(FuFu)Ay+(FuFu)Ax=(pp)Ay+D(uu)_D(uu)eewwnn

3、ss2weeEPwPW+D(uu)-D(uu)+(S+Su)AxAynNPsPSCPP对流项采用一阶迎风格式，扩散项采用中心差分，F=(pu)Ay,wwF=(pv)Ax,nnWPnF=(pu)Ay,eeF=(pu)Ay,ssHAye5xPESPPNau=au+au+au+au+(p一p)Ay+bPPWWEESSNN2WEa=D+F,a=DWwwEeva=D+F,a=DFSs2sNn2na=a+a+a+a+(FF)+(FF)PWESNewnsb=Ay(P一P)2WE步骤1：依次求解动量离散方程au*eeau*+A(p*-p*)+bnbnbePEav*+A(p*-p*)+bnbnbnPN4赋值：u

4、*v*步骤2：根据速度u求解压力修正方程p0Jpap'=ap'+ap'+ap'+ap'+b'PPWWEESSNNu0Jub'=(pu*)-(pu*)Ay+(pv*)-(pv*)Axwesn90Hav*nn步骤3：对压力和速度进行修正p二p*+p'u=u*+u'=u*+d(p'-p')eeeeePEv二v*+v'二v*+d(p'-p')nnn"八-nPN步骤4：求解其他变量°的离散方程（视需要进行）a°二'a，,+bppnbnb收敛否?是*结束TD

5、MA算法0=C11一卩0+D0一a0=C2122232一卩0+D0-a0=C3233343一卩0+D0-a0=C4344454B0+D0-a0=Cnn-1nnnn+1n0=Cn+1n+1在上式中，假定01和0+1是边界上的值，为己知。上式中任1n+1-B0+D0-a0=Cjj-1jjjj+1j除第一及最后一个方程外，其余方程可写为：aBC20+D32a30+D3a0+D54af0Dn+1n方程都可写成：02=030420+Di2虫0+D23巳0+D34+巳0Dn2-D2C3-3C4C+n-n-1Dn这些方程可通过消元和回代两个过程求解。(卩3JD03=04+B0+CX0+丄'1D丿22

6、丿a3D丿+C3现引入记号：aDC'2(则，03C»+2-iD20+4即：jD-卩A丿332a3,C'D-BA33320=A0+C'jjj+1j(Bc'+c)323jD-卩A丿332BC'+C323D-BA332a卩C+CA=z,C'=jljD卩AzD卩AjjjTjjjT在边界点，j=l与j=n+1,A二0,C'=0111A=0,C'=0n+1n+1n+1为了求解方程组，首先要对方程组按的形式编排，并明确其中的系数.0D和C。从jjjjj=2起，计算出C'z和Az，直到j=n。由于在边界位置（n+1）0的数值是已

7、知的，因此，可连续计算出0。j3、计算结果对比本文基于simple算法用编写MATLAB程序对整个流场进行计算，另外借助fluent计算流体力学软件数值模拟，通过对比分析计算结果，得出整个流场的流速矢量图和速度云图。上图为两种计算方法的计算结果，发现两者速度分布趋势一致，编程计算结果数值偏大，原因在于计算采用一阶迎风格式，节点数值趋近于其上游节点值。下图为速度云图，入口段效应的影响导致流动未达到充分发展。RO2370血D1£flLKWnri.a«z-129IS.zrzDICQDDBS00*1DOOTM4*-1|42i«11Wflu3Ri"1.047MAT

8、LAB程序：%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);Ul_old=zeros(a,b);Pl_old=zeros(a,b);%初值Ul_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;whilee220.01%整个流场迭代ifm>4breakend

9、m=m+1;forj=2:1:b-1;ifj>2breakendn=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);whilee0.001%

10、迭代，算两列，其他为已知%对(i,j)点求uxin%fori=3:1:a-2Fe1(i,j)=den*(Ul_old(i,j)+Ul_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-

11、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).+(Pl_old(i,j)-Pl_old(i,j+1)*dy/2;%(i,j+1)点的源项endAl_u(1,1)=1;%附加点Al_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;Ul_new(:,j)=in

12、v(Al_u)*b1_u(:,j);%得到(i,j+1)点的速度计算值%对(匚，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-2uw_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);endP1_fix(1,j)=(uw_sum(1,j)-ue_sum(1,j)/de1_sum(1,j);%对(匚j)点进行速度修正%fori=3:1:a-2Ul_fix(i,j)=(dy/ap_u1(i,j)*Pl_fix(1,j);%速度修正endU1_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>1000breakendA(a-4)*m