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

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

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

3、編寫MATLAB程序對整個流場進行計算，另外借助fluent計算流體力學軟件數值模擬，通過對比分析計算結果，得出整個流場的流速矢量圖和速度云圖。上圖為兩種計算方法的計算結果，發(fā)現(xiàn)兩者速度分布趨勢一致，編程計算結果數值偏大，原因在于計算采用一階迎風格式，節(jié)點數值趨近于其上游節(jié)點值。下圖為速度云圖，入口段效應的影響導致流動未達到充分發(fā)展。 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é)點%邊界條件的設置、初始值的設置%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%整個流場迭代% if m>4 break end m=m+1;for j=2:1:b-1; if j>2 break end n=1;%迭代次數% 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%迭代，算兩列，其他為已知%對（i,j）點求uxin% for i=3:1:a-2 Fe1(i,j)=den*(U1_old(i,j)+U1_old(i,j+1)*dy/2;%對流強度% 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;%離散方程對應系數% 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)點的源項% end A1_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(

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);%對（i，j）點進行速度修正% 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