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

1、问题描述图1为20℃的水在长度为150mm，宽度为10mm的管道中流动，流入管道速度为0.2m/s，流出管道背压为1atm，基于simple算法对整个流场进行计算，计算管长100mm处流速并与fluent计算结果对比。图1流动系统流动计算：该流动问题为二维定常无内热源不可压缩层流流动。控制方程离散在二维直角坐标系中，对流—扩散方程的通用形式为：图2网格编号针对本问题，其连续性方程和动量方程为：交错网格下动量方程的离散：上式积分得到：对流项采用一阶迎风格式，扩散项采用中心差分，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>4breakendm=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);whilee>0.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;ifn>1000breakendA((a-4)*m+1:(a-4)*(m+1),2:b-1)=U1_old(3:a-2,2:b-1);%赋值，即m从1开始%e