中科大計算流體力學CFD之大作業(yè)二

CFD實驗報告二姓名： 學號：一、題目求解Poisson方程 ， ， ， ， ， ， ，描出等值線：, 0.2, 0.5, 0.75, 1. 要求所用方法：(1) Jacobi, G-S選一；SOR，線SOR，塊SOR選一。迭代法要求誤差； (2) CG方法，MG方法選一。二、報告要求1）簡述問題的性質、求解原則；2）列出全部計算公式和步驟；3）表列出程序中各主要符號和數組意義；4）計算結果與精確解比較5）結果分析（方案選擇、比較、討論、體會、建議等）；6）附源程序。三、問題簡要分析3.1問題性質從題目給出的問題可以看出，該方程是一個二階線性非齊次偏微分方程，并且給出了四個邊界條件，更具定解條件，所以該方程是可以求解的。3.2求解原則題中是關于x和y的二階導數，因此可以用二階中心差分離散，即用正五點格式離散泊松方程。將差分方程整理成以五對角矩陣為系數的線性方程組，用Jacobi或者CG或者SOR等迭代方法求解線性方程組，即可得到函數的在網格點處的離散值。同時，計算域為的矩形區(qū)域，劃分結構網格，均勻步長，設方向網格步長，方向網格步長，則方向網格點數為m=， 方向網格點數為n=。四、計算公式和步驟4.1 精確解的計算題中已知四個邊界條件，可以通過現(xiàn)已有的解析求解不難求得精確解（解析解）為。4.2迭代求解一般過程迭代法是求解離散代數方程組的主要方法。假如我們對題目中的PDE給定一個離散格式，則對每一個確定的離散點（i ，j），PDE轉化為FDE，為，其中是離散格式中所有與有關的點，函數中所有元素都是線性疊加，即是線性多元函數。所以，將邊界上給定結果以及中間的離散格式得到的結果綜合起來就是一個線性方程組如下：也就是，令,故,即其中N可以分為對角陣、三對角陣、上下三角陣。迭代式：,給出初值或者,殘量推出收斂準則：1）n足夠大，殘量趨于零（小于） 2）4.3方程離散上述泊松方程在的正五點差分格式為：將泊松方程在計算域網格點上進行差分得到差分方程：以下對邊界條件進行離散：綜合差分方程及邊界條件將方程整理成線性方程組的形式得到：其中，；4.4線性方程組迭代算法對于線性方程組，可以利用以下迭代算法進行求解。4.4.1 Jacobi迭代將系數矩陣分解，Jacobi迭代格式為：其中，。4.4.2超松弛迭代法(Successive Over-Relaxation)將系數矩陣分解，SOR迭代格式為：其中，為迭代松弛因子，當時，松弛迭代法就是Gauss-Seidel迭代法；當時被稱為超松弛迭代。4.4.3共軛梯度法(Conjugate Gradient)共軛梯度法求解代數方程組的算法2為：（1）假定初場；（2）計算初余量；（3）令；（4）對，直到（記） ，結束。五、結果比較與分析由上面已經羅列的公式可知，精確求解的泊松方程為：，利用MATLAB編寫程序，步長取=0.01分別用Jacobi迭代、超松弛迭代法、共軛梯度法進行計算，分別描出等值線：, 0.2, 0.5, 0.75, 1，并與精確解對比，如下圖所示：圖1 精確解的等值線圖2 Jacobi迭代法的等值線圖3超松弛迭代法的等值線圖4共軛梯度法的等值線綜合來看，由圖1、2、3、4可以看到，差分解與數值解基本相同。此外，給出整個求解區(qū)域的數值解和精確解的云圖，如下所示：圖5 數值解云圖圖6精確解云圖定義數值解與精確解之間的誤差為，計算得到：根據題意以迭代誤差小于為迭代收斂條件，將各個迭代算法求解整個問題的時間，求解線性方程組的迭代步數，計算誤差列表如下：JacobiSORCG直接求解耗費時間(s)0.785553.05530.04930.0321迭代步數74559961380精度8.0261E-045.6117E-052.6E-031.1457E-07其中：直接求解是MATLAB中求解線性方程組的左除算法，對于，則，其結果相當于；同時，SOR方法中的松弛因子。綜上所述：1. 三種方法的數值解與精確解差別均很小，該算例驗證了三種方法的正確性。2. 根據題意采用的是均勻網格，從云圖上我們可以看出，不同位置的等值線密度不同。因此，若采用非均勻網格，在等值線比較密的位置加密網格可以提高計算的精度。3. 不同迭代算法效率不同，迭代步數最少的是CG方法，迭代步數最多的是Jacobi迭代，各個算法的迭代步數比較與理論相符；其中MATLAB自帶求解線性方程組的左除算法效率和精度均最高，其次是共軛梯度法，Jacobi迭代，SOR算法；4. CG的迭代步數雖然較少，但是精度有待提高，可以通過提高算法收斂標準提高精度；5. 程序采用的是正方形網格，未研究兩個方向網格數目不同以及漸近網格對計算結果的影響?？梢赃M一步研究其對計算結果和計算效率的影響。六、源程序及其主要符號說明符號說明如下表所示：lx,ly求解區(qū)間大小dx，dyx,y方向上的網格步長m，nx,y方向上的節(jié)點數K分塊系數矩陣B非齊次項矩陣X解向量ps不考慮邊界的解Ps考慮邊界的解Ps_e精確解kk迭代次數Q計算精度源程序如下：1. 主程序：clear %清空global kk;kk=0;dx=0.01;dy=0.01; %設置網格步長lx=1;ly=1; %設置計算域大小m=lx/dx+1;n=ly/dy+1; %計算x、y方向節(jié)點數a=dx2;b=dy2;I=a*eye(m-2);%生成矩陣II=sparse(I);G=zeros(m-2,m-2);%生成矩陣GG(1,1:2)=-2*(a+b) b;for i=2:m-3 G(i,i-1:i+1)=b -2*(a+b) b;endG(m-2,m-3:m-2)=b -2*(a+b);G=sparse(G);M=(m-2)*(n-2);%生成矩陣BB=zeros(M,1);k=1;for i=2:m-1 for j=2:n-1 Xi=(i-1)*dx; Yj=(j-1)*dy; B(k)=dx2*dy2*sin(Xi)*cos(Yj); if j=2 B(k)=B(k)+a*sin(Xi)/2; end if i=m-1 B(k)=B(k)-b*(Yj-sin(1)*cos(Yj)/2); end if j=n-1 B(k)=B(k)-a*(Xi-sin(Xi)*cos(1)/2); end k=k+1; endendK=cell(n-2);%生成分塊系數矩陣KK(:)=zeros(m-2,m-2);K(1,1:2)=G I;for i=2:n-3 K(i,i-1:i+1)=I G I;endK(n-2,n-3:n-2)=I G;K=cell2mat(K);disp(SC李文建);disp( );tic %開始計時并迭代求解線性方程組% X=KB; %左除算法X=Jacobi(K,B);% X=SOR(K,B);% X=CG(K,B);ps=cell(n-2,1);ps(:)=zeros(1,m-2);for k=1:n-2 ps(n-2-k+1)=X(m-2)*(k-1)+1:(m-2)*k);endps=cell2mat(ps);Ps=zeros(n,m);Ps(2:n-1,2:m-1)=fliplr(rot90(ps);Ps(:,1)=zeros(n,1);%加入邊界值Ps(:,m)=n-1:-1:0*dy-sin(1)*cos(n-1:-1:0*dy)/2;Ps(n,:)=-sin(dx*0:m-1)/2;Ps(1,:)=(dx*0:m-1)-sin(dx*0:m-1)*cos(1)/2;toc %結束計時Ps_e=fun(m,n,dx,dy); %精確解Q=sum(sum(abs(Ps-Ps_e)/(m*n); %平均誤差大小fprintf(迭代次數為： %8.0fn, kk);disp(計算精度為：);Qcontour(flipud(Ps),0.05,0.2 0.5 0.75 1,ShowText,on) %畫出數值解的等值線% contour(flipud(Ps_e),0.05,0.2 0.5 0.75 1,ShowText,on) %畫出精確解的等值線% contourf(flipud(Ps_e) %畫出精確解分布云圖% contourf(flipud(Ps) %畫出數值解分布云圖2. 求精確解函數function Ps_e=fun(m,n,dx,dy)x=0:m-1*dx;y=0:n-1*dy;Ps_e=-sin(x)*cos(y)/2+x*y;Ps_e=rot90(Ps_e,1);end3. Jacobi迭代函數function X=Jacobi(A,b)global kk;D=diag(diag(A);%設置Jacobi迭代法需要的矩陣L=D-tril(A);U=D-triu(A);N=length(b);X0=zeros(N,1);B=D(L+U);R=0.75;if R=1e-6 X0=X; X=B*X0+f; i=i+1; if i=1E6 disp(迭代步數太多，放棄計算！) break; end end kk=i; %輸出迭代步數else disp(不收斂！) X=;endend4. SOR算法函數function X=SOR(A,b) global kk;D=diag(diag(A);L=D-tril(A);U=D-triu(A);N=length(b);X0=zeros(N,1);w=1.75;B=(D-w*L)(1-w)*D+w*U);R=0.75; if R=1e-6 X0=X; X=B*X0+f; i=i+1; if i=1E6 disp(迭代步數太多，放棄計算！) break; end end kk=i; %輸出迭代步數else disp(不收斂！) X=;end5. CG方法函數function X=CG(A,b) global kk;N=length(b);x=zeros(N,1);eps=1e-6;r=b-A*x;d=r;for k=0:N-1 a=(norm(r)2)/(d*A*d); x=x+a*d; r_t=b-A*x; if (norm(r_t)=eps)|(k=N-1) break; end B=(norm(r_t)2)/(norm(r)2); d=r_t+B*d; r=r_t;en