中科大计算流体力学CFD之大作业二(共14页).doc

《中科大计算流体力学CFD之大作业二(共14页).doc》由会员分享，可在线阅读，更多相关《中科大计算流体力学CFD之大作业二(共14页).doc（14页珍藏版）》请在淘文阁 - 分享文档赚钱的网站上搜索。

1、精选优质文档-倾情为你奉上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问题性质从题目给出的问题可以看出，该方程是一个二阶线性非齐次偏微分方程，并且给出

2、了四个边界条件，更具定解条件，所以该方程是可以求解的。3.2求解原则题中是关于x和y的二阶导数，因此可以用二阶中心差分离散，即用正五点格式离散泊松方程。将差分方程整理成以五对角矩阵为系数的线性方程组，用Jacobi或者CG或者SOR等迭代方法求解线性方程组，即可得到函数的在网格点处的离散值。同时，计算域为的矩形区域，划分结构网格，均匀步长，设方向网格步长，方向网格步长，则方向网格点数为m=， 方向网格点数为n=。四、计算公式和步骤4.1 精确解的计算题中已知四个边界条件，可以通过现已有的解析求解不难求得精确解（解析解）为。4.2迭代求解一般过程迭代法是求解离散代数方程组的主要方法。假如我们对题

3、目中的PDE给定一个离散格式，则对每一个确定的离散点（i ，j），PDE转化为FDE，为，其中是离散格式中所有与有关的点，函数中所有元素都是线性叠加，即是线性多元函数。所以，将边界上给定结果以及中间的离散格式得到的结果综合起来就是一个线性方程组如下：也就是，令,故,即其中N可以分为对角阵、三对角阵、上下三角阵。迭代式：,给出初值或者,残量推出收敛准则：1）n足够大，残量趋于零（小于） 2）4.3方程离散上述泊松方程在的正五点差分格式为：将泊松方程在计算域网格点上进行差分得到差分方程：以下对边界条件进行离散：综合差分方程及边界条件将方程整理成线性方程组的形式得到：其中，；4.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编写程序，步长取

5、=0.01分别用Jacobi迭代、超松弛迭代法、共轭梯度法进行计算，分别描出等值线：, 0.2, 0.5, 0.75, 1，并与精确解对比，如下图所示：图1 精确解的等值线图2 Jacobi迭代法的等值线图3超松弛迭代法的等值线图4共轭梯度法的等值线综合来看，由图1、2、3、4可以看到，差分解与数值解基本相同。此外，给出整个求解区域的数值解和精确解的云图，如下所示：图5 数值解云图图6精确解云图定义数值解与精确解之间的误差为，计算得到：根据题意以迭代误差小于为迭代收敛条件，将各个迭代算法求解整个问题的时间，求解线性方程组的迭代步数，计算误差列表如下：JacobiSORCG直接求解耗费时间(s)

6、0.785553.05530.04930.0321迭代步数74559961380精度8.0261E-045.6117E-052.6E-031.1457E-07其中：直接求解是MATLAB中求解线性方程组的左除算法，对于，则，其结果相当于；同时，SOR方法中的松弛因子。综上所述：1. 三种方法的数值解与精确解差别均很小，该算例验证了三种方法的正确性。2. 根据题意采用的是均匀网格，从云图上我们可以看出，不同位置的等值线密度不同。因此，若采用非均匀网格，在等值线比较密的位置加密网格可以提高计算的精度。3. 不同迭代算法效率不同，迭代步数最少的是CG方法，迭代步数最多的是Jacobi迭代，各个算法的

7、迭代步数比较与理论相符；其中MATLAB自带求解线性方程组的左除算法效率和精度均最高，其次是共轭梯度法，Jacobi迭代，SOR算法；4. CG的迭代步数虽然较少，但是精度有待提高，可以通过提高算法收敛标准提高精度；5. 程序采用的是正方形网格，未研究两个方向网格数目不同以及渐近网格对计算结果的影响。可以进一步研究其对计算结果和计算效率的影响。六、源程序及其主要符号说明符号说明如下表所示：lx,ly求解区间大小dx，dyx,y方向上的网格步长m，nx,y方向上的节点数K分块系数矩阵B非齐次项矩阵X解向量ps不考虑边界的解Ps考虑边界的解Ps_e精确解kk迭代次数Q计算精度源程序如下：1. 主程

8、序：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方向节点数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

9、=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

10、: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

11、(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) %画出数值解的等值线%

12、 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);%设

13、置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);

14、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)=

15、eps)|(k=N-1) break; end B=(norm(r_t)2)/(norm(r)2); d=r_t+B*d; r=r_t;endX=x;kk=k; %迭代次数end七、参考文献1 吴清松. 计算热物理引论M. 中国科学技术大学出版社. 合肥.20092 孔凡铀. 用FFT法与SOR法求解二维泊松方程的比较J. 气象科学. 1988年第一期3 杨金凤,邓居智,陈辉. 多重网格法在求解泊松方程中的应用进展J. 内蒙古石油化工. 2011年第24期4 王忆锋, 唐利斌. 利用有限差分和MATLAB矩阵运算直接求解二维泊松方程J. 红外技术. 第32卷第4期. 2010.04专心-专注-专业