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

1、问题描述图1为20℃的水在长度为150mm，宽度为10mm的管道中流动，流入管道速度为0.2m/s，流出管道背压为1atm，基于simple算法对整个流场进行计算，计算管长100mm处流速并与fluent计算结果对比。图1流动系统流动计算：363650998.21/;1001.6110;998.210.21010Re23001001.6110dmmkgmPasud该流动问题为二维定常无内热源不可压缩层流流动。2、控制方程离散在二维直角坐标系中，对流—扩散方程的通用形式为：uStxyxxyy图2网格编号针对本问题，其连续性方程和动量方程为：ypyvyxvxyvvxvuxpyuyxuxyuvxuuyvxu)()()()(0)()(交错网格下动量方程的离散：()()neneneneswswswswneneuswswpuuudxdyvudxdydxdyudxdyxyxxxuudxdySdxdyyy上式积分得到：1()()((2((eewwnnssweeEPwPWnNPsPSCPPFuFuyFuFuxppyDuuDuuDuuDuuSSuxy）-））-）+对流项采用一阶迎风格式，扩散项采用中心差分，(),(),(),(),we，，WE1()2,11,22()()1y()2PPWWEESSNN否是赋值：0pp假设一个速度初场（其它变量的初场是否需要视情况而定）假设一个压力场，即给定压力猜测值开始根据当前的已知量，计算动量离散方程等方程中的系数和常数项步骤1：依次求解动量离散方程********()()eenbnbePEnnnbnbnPNauauAppbavavAppb步骤2：根据速度,uv求解压力修正方程PPWWEESSNNapapapapapb''''''****[()()][()()]wesnbuuyvvx'步骤3：对压力和速度进行修正*****()()eeeeePEnnnnnPNpppuuuudppvvvvdpp'''''''步骤4：求解其他变量的离散方程（视需要进行）ppnbnbaab收敛否？结束赋值：000ppuuTDMA算法11212223232333434344454CDCDCDC1111nnnnnnnnnDCC在上式中，假定1和1n是边界上的值，为己知。上式中任一方程都可写成：11jjjjjjjDC除第一及最后一个方程外，其余方程可写为：222231222333342333444453444CDDDCDDDCDDD11nnnnnnnnnCDDD这些方程可通过消元和回代两个过程求解。223132233422333322++CCDDDDDD现引入记号：'222221222,+CACDDD则，'332334332332+CCDADA''332333332332+,CCACDADA即：'1jjjjAC'1'11+,jjjjjjjjjjjjCCACDADA在边界点，j=1与j=n+1，'1110,AC'1110,nnnAC为了求解方程组，首先要对方程组按的形式编排，并明确其中的系数jjjD和jC。从j=2起，计算出'jC和jA，直到j=n。由于在边界位置（n+1）的数值是已知的，因此，可连续计算出j。3、计算结果对比本文基于simple算法用编写MATLAB程序对整个流场进行计算，另外借助fluent计算流体力学软件数值模拟，通过对比分析计算结果，得出整个流场的流速矢量图和速度云图。0.00000.00250.00500.00750.01000.01250.000.050.100.150.200.250.30uy数值计算matlab计算上图为两种计算方法的计算结果，发现两者速度分布趋势一致，编程计算结果数值偏大，原因在于计算采用一阶迎风格式，节点数值趋近于其上游节点值。下图为速度云图，入口段效应的影响导致流动未达到充分发展。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;whilee220.01%整个流场迭代%ifm4breakendm=m+1;forj=2:1:b-1;ifj2breakendn=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%迭代，算两列，其他为已知%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%对（i,j）点求uxin%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%fori=3:1:a-2Fe1(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)点的源项%endA1_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-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);%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%对（i，j）点进行速度修正%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%fori=3:1:a-2U1_fix(i,j)=(dy/ap_u1(i,j))*P1_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;ifn1000breakendA((a-4)*m+1:(a-4)*(m+1),2:b-1)=U1_old(3:a-2,2:b-1);%赋值，即m从1开始%endendE(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);uplot(y,u)