中科大计算流体力学CFD之大作业二

1、CFD 实验报告二学号:、题目求解 Poisson 方程2 2Y 2sinxcosy， x y0 x 1，0 y 1，|x 00，|x 1ysin 1cosy2 ，sinx|y 02，|y 1xsin xcos12 ，描出等值线：0.05,0.2, 0.5, 0.75, 1.要求所用方法：Jacobi, G-S 选一；SOR,线 SOR,块 SOR 选一。迭代法 要求误差 106；（2） CG 方法，MG 方法选一。二、报告要求1）简述问题的性质、求解原则；2）列出全部计算公式和步骤；3）表列出程序中各主要符号和数组意义；4）计算结果与精确解比较5）结果分析（方案选择、比较、讨论、体会、建议等

2、）；6）附源程序。三、问题简要分析 3.1 冋题性质从题目给出的问题可以看出，该方程是一个二阶线性非齐次偏微分方程，并 且给出了四个边界条件，更具定解条件，所以该方程是可以求解的。3.2 求解原则题中是关于 x 和 y 的二阶导数，因此可以用二阶中心差分离散，即用正五点格式离散泊松方程。将差分方程整理成以五对角矩阵为系数的线性方程组，用Jacobi 或者 CG 或者 SOR 等迭代方法求解线性方程组，即可得到函数 的在网格点处的离散值。同时，计算域为 Lx Ly 的矩形区域，划分结构网格，均匀步 长，设 x 方向网格步长 x，y方向网格步长 y,则 x 方向网格点数为m=Lx/ x 1，y方向

3、网格点数为 n= Ly / y 1。四、计算公式和步骤4.1 精确解的计算 题中已知四个边界条件，可以通过现已有的解析求解不难求得精确解(解析、 1解)为sinx cosy xy。4.2 迭代求解一般过程 迭代法是求解离散代数方程组的主要方法。假如我们对题目中的 PDE 给定一个离散格式，则对每一个确定的离散点(i，j)，PDE 转化为 FDE，为i,jF1(a,b| c,d,si n(i 1) XCOS(i 1) y)，其中a,bc,d是离散格式中所有与i,j有关的点，函数 F1中所有元素都是线性叠加，即 F1是线性多元函数。所以, 将边界上给定结果以及中间的离散格式得到的结果综合起来就是一

4、个线性方程组如下：匀 M III III III 川 51,10hhih4IfI 4p*1,n 1sin1 1RiI 4HIibhih4*I I*4HIIb4ih4PIIP4% 1)2,1I | | Id IDSn1,n1sin 1cOs1也就是 AV B，A 0，令 A N P，|N| IIAI故 AV PV B PV 即 NV B PVNVn1PVnB其中 N 可以分为对角阵、三对角阵、上下三角阵。迭代式：Vn 1N1PVnN1B ,给出初值 V0或者 NVn1PVnB PVnAVnB AVn,RnB AVn 残量推出 NVn1NVnRn收敛准则：1） n 足够大，残量趋于零（小于 ）2）

5、 Vn 1Vn0（）4.3 方程离散上述泊松方程在（xyj）的正五点差分格式为:i 1,j2i,ji 1,ji,j 12i,ji,j 122Xy将泊松方程在计算域网格点sin xicos yj2 y sin x2cosy22 y sinx2cosy32 2 x y sin x2cos yn2 y sin X3cosy22 y sin x3cosy32 2 .x y sin x3cosyn: 2 .y sin Xm 1cos y2: 2 .y sin Xm icos *2 2x y sin Xm icosy1,2n)(k 1) y,(k 1,2n)1) x,( k 1,2.m)(k 1) x,(

6、k 1,2.m)（Xi, Yj）上进行差分得到差分方程:2,j2y2(3,222,21,2)x2(2,322,22,1)2X2,j3y2(3,322,31,3)x2(2,4212,32,2)2X2,j n 1y2(3,n 122,n 11,n 1)x2(2,n22,n 12,n2)3,j2y2(4,223,22,2)x2(3,323,23,1)2X3,j3y2(4,323,32,3)x2(3,423,33,2)2X13,jn 1y2(4,n 123,n 12,n 1)x2(3,n123,n 13,n2)i m 1,j2厶/y(m,22m 1,2m 2,2)i m 1,j3ay (m,32m 1

7、,3m 2,3)i m 1, j n 12y (m,n 12-m 1,n 1m 2,n1)以下对边界条件进行离散:左边界 i1:1,k0,(k右边界 i m:m, kyksin 1cos2yk,yk下边界 j1:k,1sin xk2,X:k(k上边界 j n :k,nXksin xkco.s1,Xk2x2(x2(m 1,3m 1,41,21,3m 1,1)m 1,2)m 1,n1,nm 1,n2)综合差分方程及边界条件将方程整理成线性方程组 KB 的形式得到:11G(m 2)(n 2) (m 2)(n2)其中,其中 a=x2b2b3bn2(a+b )b1(m 2)b2(a+b )t2,n2,2

8、b2b31,23,2HI(m2)2(a+b )b2(a+b )(m 2) (m2),其中 a=2 2x ,b y ;3,n 1HIbn 1ab sin x2cos y2mHI1,3aabsin x3cosy2absin xm 1cos y2ab sin X2cosybabsin X3cosy31,22,1a1,3m 1,n2,3III1,32,n 1卅Tm 1,n 1(m2)(n 2) 1b1,23,1m 1,1absin xm 1cosy3abs inx2cosyn 1absin沁cosyn 12,na3,n1,n 1,其中 a= x2,b y2;abs in xm 1cosyn 1m 1,

9、n114.4 线性方程组迭代算法对于线性方程组 AX b，可以利用以下迭代算法进行求解4.4.1Jacobi 迭代将系数矩阵分解 A D L U，Jacobi 迭代格式为：x(k 1)Bx(k)f其中 B D1(L U)， f D1b。4.4.2 超松弛迭代法(Successive Over-Relaxation)将系数矩阵分解 A D L U，SOR 迭代格式为:(k 1)(k)x Bx f其中 B (D wL)1(1 w)D wU )， f w(D wL)1b，w 为迭代松弛因子,当 w=1 时，松弛迭代法就是 Gauss-Seide 迭代法；当 w 1 时被称为超松弛迭代。4.4.3 共

10、轭梯度法(Conjugate Gradient)共轭梯度法求解代数方程组 AX b 的算法2为:(1 )假定初场 X(0)；(2)计算初余量 r(0)b AX(0)；令 p(0)=r(0)；(4)对 k 0,1,2 川，直到 r(k)tolerenee (记(a,b) aT|b)k=(r(k),r(k)/(Ap(k), p(k)(k 1)r(k)rkAp(k)(k+1)(k+1)(k)(k)、k=(r,r)/(r,r )(k 1)(k 1)(k)PrkP，结束X(k 1)X(k)kP(k)五、结果比较与分析1由上面已经罗列的公式可知，精确求解的泊松方程为： =-?sin xcosy xy，利用

11、 MATLAB 编写程序，步长取 x= y=0.01 分别用 Jacobi 迭代、超松弛迭代法、共轭梯度法进行计算，分别描出等值线：0.05, 0.2, 0.5, 0.75,1，并与精确解对比，如下图所示：图 1 精确解的等值线图 2 Jacobi 迭代法的等值线图 3 超松弛迭代法的等值线图 4 共轭梯度法的等值线综合来看，由图 1、2、3、4 可以看到，差分解与数值解基本相同此外，给出整个求解区域的数值解和精确解的云图，如下所示:图 5 数值解云图图 6 精确解云图1定义数值解与精确解之间的误差为E = cal*,其中*为精确解，计nm1E=一cai*8.0261 10-4（Jacobi

12、迭代结果）nm根据题意以迭代误差小于 106为迭代收敛条件，将各个迭代算法求解整个问题的时间，求解线性方程组 AX b 的迭代步数，计算误差列表如下：JacobiSORCG直接求解耗费时间（S）0.785553.05530.04930.0321迭代步数74559961380精度8.0261E-045.6117E-052.6E-031.1457E-07其中：直接求解是 MATLAB 中求解线性方程组的左除算法，对于 AX b , 则 X A b，其结果相当于 X A-1b ;同时，SOR 方法中的松弛因子 w 取 1.75。综上所述：1. 三种方法的数值解与精确解差别均很小，该算例验证了三种方法

13、的正确性。2. 根据题意采用的是均匀网格，从云图上我们可以看出，不同位置的等值线 密度不同。因此，若采用非均匀网格，在等值线比较密的位置加密网格可以提 高计算的精度。3. 不同迭代算法效率不同，迭代步数最少的是 CG 方法，迭代步数最多的是Jacobi 迭代，各个算法的迭代步数比较与理论相符；其中 MATLAB 自带求解线性方程组的左除算法效率和精度均最高，其次是共轭梯度法，Jacobi 迭代，SOR算法；4. CG 的迭代步数虽然较少，但是精度有待提高，可以通过提高算法收敛标准 提高精度；5. 程序采用的是正方形网格， 未研究两个方向网格数目不同以及渐近网格对 计算结果的影响。可以进一步研究

14、其对计算结果和计算效率的影响。六、源程序及其主要符号说明符号说明如下表所示：lx,ly求解区间大小dx， dyx,y 方向上的网格步长m nx,y 方向上的节点数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 方向节点数a=dxA2;b=dyA2;I=a*eye(m-2);% 生成矩阵 II=sparse(I);G

15、=zeros(m-2,m-2);%生成矩阵 GG(1,1:2)=-2*(a+b) b;for i=2:m-3G(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-1for j=2:n-1Xi=(i-1)*dx;Yj=(j-1)*dy;B(k)=dxA2*dyA2*sin(Xi)*cos(Yj);if j=2B(k)=B(k)+a*sin(Xi)/2;endif i=m-1 B(k)=B(k)-b*(Yj-sin(1)*

16、cos(Yj)/2);endif j=n-1B(k)=B(k)-a*(Xi-sin(Xi)*cos(1)/2); endk=k+1;endendK=cell(n-2);%生成分块系数矩阵 K K(:)=zeros(m-2,m-2);K(1,1:2)=G I;for i=2:n-3K(i,i-1:i+1)=I G I;endK(n-2,n-3:n-2)=I G;K=cell2mat(K);disp( );tic %开始计时并迭代求解线性方程组% X=KB;% 左除算法X=Jacobi(K,B);% X=SOR(K,B);% X=CG(K,B);ps=cell(n-2,1);ps(:)=zeros

17、(1,m-2);for k=1:n-2ps(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(

18、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,

19、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=1E6disp(迭代步数太多，放弃计算！) break;end end kk=i; %输出迭代步数 elsedisp(不收敛！)X=;endend4.SOR 算法函数 function X=SOR(A,b) global kk;D=diag(dia