基于simple算法的流場模擬計算

1、問題描述圖1為20℃的水在長度為150mm，寬度為10mm的管道中流動，流入管道速度為0.2m/s，流出管道背壓為1atm，基于simple算法對整個流場進行計算，計算管長100mm處流速并與fluent計算結(jié)果對比。圖1流動系統(tǒng)流動計算：該流動問題為二維定常無內(nèi)熱源不可壓縮層流流動?？刂品匠屉x散在二維直角坐標系中，對流—擴散方程的通用形式為：圖2網(wǎng)格編號針對本問題，其連續(xù)性方程和動量方程為：交錯網(wǎng)格下動量方程的離散：上式積分得到：對流項采用一階迎風格式，擴散項采用中心差分，MATLAB程序：%%%%%%%%%%%%%%%%%%%%%%%%dx=5e-4;dy=5e-4;den=998;dyna=1001.6e-6;a=0.01/dy+3;%x節(jié)點%b=0.15/dx+1;%y節(jié)點%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%邊界條件的設(shè)置、初始值的設(shè)置%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%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;%迭代次數(shù)%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;%離散方程對應(yīng)系數(shù)%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