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

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