1、1、問題描述圖1為20的水在長度為150mm，寬度為10mm的管道中流動(dòng)，流入管道速度為0.2m/s，流出管道背壓為1atm，基于simple算法對(duì)整個(gè)流場進(jìn)行計(jì)算，計(jì)算管長100mm處流速并與fluent計(jì)算結(jié)果對(duì)比。圖1 流動(dòng)系統(tǒng)流動(dòng)計(jì)算：該流動(dòng)問題為二維定常無內(nèi)熱源不可壓縮層流流動(dòng)。2、 控制方程離散在二維直角坐標(biāo)系中，對(duì)流擴(kuò)散方程的通用形式為：圖2 網(wǎng)格編號(hào)針對(duì)本問題，其連續(xù)性方程和動(dòng)量方程為：交錯(cuò)網(wǎng)格下動(dòng)量方程的離散：上式積分得到：對(duì)流項(xiàng)采用一階迎風(fēng)格式，擴(kuò)散項(xiàng)采用中心差分，假設(shè)一個(gè)速度初場（其它變量的初場是否需要視情況而定）假設(shè)一個(gè)壓力場，即給定壓力猜測值開始根據(jù)當(dāng)前的已知量，計(jì)算

2、動(dòng)量離散方程等方程中的系數(shù)和常數(shù)項(xiàng)步驟1：依次求解動(dòng)量離散方程步驟2：根據(jù)速度求解壓力修正方程步驟3：對(duì)壓力和速度進(jìn)行修正步驟4：求解其他變量的離散方程（視需要進(jìn)行）收斂否？結(jié)束賦值：賦值：否是TDMA算法在上式中，假定和是邊界上的值，為己知。上式中任一方程都可寫成：除第一及最后一個(gè)方程外，其余方程可寫為： 這些方程可通過消元和回代兩個(gè)過程求解。現(xiàn)引入記號(hào)：則，即：在邊界點(diǎn)，j=1與j=n+1，為了求解方程組，首先要對(duì)方程組按的形式編排，并明確其中的系數(shù)和。從j=2起，計(jì)算出和，直到j(luò)=n。由于在邊界位置（n+1）的數(shù)值是已知的，因此，可連續(xù)計(jì)算出。3、 計(jì)算結(jié)果對(duì)比本文基于simple算法用

3、編寫MATLAB程序?qū)φ麄€(gè)流場進(jìn)行計(jì)算，另外借助fluent計(jì)算流體力學(xué)軟件數(shù)值模擬，通過對(duì)比分析計(jì)算結(jié)果，得出整個(gè)流場的流速矢量圖和速度云圖。上圖為兩種計(jì)算方法的計(jì)算結(jié)果，發(fā)現(xiàn)兩者速度分布趨勢一致，編程計(jì)算結(jié)果數(shù)值偏大，原因在于計(jì)算采用一階迎風(fēng)格式，節(jié)點(diǎn)數(shù)值趨近于其上游節(jié)點(diǎn)值。下圖為速度云圖，入口段效應(yīng)的影響導(dǎo)致流動(dòng)未達(dá)到充分發(fā)展。 MATLAB程序：%dx=5e-4;dy=5e-4;den=998;dyna=1001.6e-6;a=0.01/dy+3;%x節(jié)點(diǎn)%b=0.15/dx+1;%y節(jié)點(diǎn)%邊界條件的設(shè)置、初始值的設(shè)置%T=0.2*ones(a,b);U1_old=zeros(a,b)

4、;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;while e22>0.01%整個(gè)流場迭代% if m>4 break end m=m+1;for j=2:1:b-1; if j>2 break end n=1;%迭代次數(shù)% e=1;%初次迭代%B=0;0;0.0823455;0.149211;0.196652;0.224194;0.23674;0.

5、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);while e>0.001%迭代，算兩列，其他為已知%對(duì)（i,j）點(diǎn)求uxin% for i=3:1:a-2 Fe1(i,j)=den*(U1_old(i,j)+U1_old(i,j+1)*dy/2;%對(duì)流強(qiáng)度% Fw1(i,j)

6、=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;%離散方程對(duì)應(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

7、(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)點(diǎn)的源項(xiàng)% end A1_u(1,1)=1;%附加點(diǎn)% A1_u(2,2)=1;%邊界點(diǎn)% A1_u(a-1,a-1)=1; A1_u(a,a)=1; b1_u(1,j)=0;%附加點(diǎn)% b1_u(2,j)=0;%邊界點(diǎn)% b1_u(a-1,j)=0; b1_u(a,j)=0; U1_new(:,j)=inv(A1_u)*b1_u(:,j);%得到(i,j+1)點(diǎn)的速度計(jì)算值% %對(duì)（i，j）點(diǎn)求p% uw_sum(1,j)=0; ue_sum(

8、1,j)=0; de1_sum(1,j)=dy/ap_u1(2,j); for i=3:1:a-2 uw_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); end P1_fix(1,j)=(uw_sum(1,j)-ue_sum(1,j)/de1_sum(1,j);%對(duì)（i，j）點(diǎn)進(jìn)行速度修正% for i=3:1:a-2 U1_fix(i,j)=(dy/ap_u1(i,j)*P1_fix(1,j);%速度修正% end U1_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; if n>1000 break end A(a-4)*m+1:(a-4)*(m+1),2:b-1)=U1_old(3:a-2,2:b-1);%賦值，即m從1開始%endendE(1:a