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

CFD实验报告二姓名：学号：一、题目求解Poisson方程yxyxcossin2222，10x，10y，0|0x，yx1|2cos1siny，0|y2sinx，xy1|21cossinx，描出等值线：05.0,0.2,0.5,0.75,1.要求所用方法：(1)Jacobi,G-S选一；SOR，线SOR，块SOR选一。迭代法要求误差610；(2)CG方法，MG方法选一。二、报告要求1）简述问题的性质、求解原则；2）列出全部计算公式和步骤；3）表列出程序中各主要符号和数组意义；4）计算结果与精确解比较5）结果分析（方案选择、比较、讨论、体会、建议等）；6）附源程序。三、问题简要分析3.1问题性质从题目给出的问题可以看出，该方程是一个二阶线性非齐次偏微分方程，并且给出了四个边界条件，更具定解条件，所以该方程是可以求解的。3.2求解原则题中是关于x和y的二阶导数，因此可以用二阶中心差分离散，即用正五点格式离散泊松方程。将差分方程整理成以五对角矩阵为系数的线性方程组，用Jacobi或者CG或者SOR等迭代方法求解线性方程组，即可得到函数,ij的在网格点处的离散值。同时，计算域为LxLy的矩形区域，划分结构网格，均匀步长，设x方向网格步长x，y方向网格步长y，则x方向网格点数为m=/1Lxx，y方向网格点数为n=/1Lyy。四、计算公式和步骤4.1精确解的计算题中已知四个边界条件，可以通过现已有的解析求解不难求得精确解（解析解）为1sincos2xyxy。4.2迭代求解一般过程迭代法是求解离散代数方程组的主要方法。假如我们对题目中的PDE给定一个离散格式，则对每一个确定的离散点（i，j），PDE转化为FDE，为,1,,(,sin[(1)]cos[(1)])ijabcdFixiy，其中,,~abcd是离散格式中所有与,ij有关的点，函数1F中所有元素都是线性叠加，即1F是线性多元函数。所以，将边界上给定结果以及中间的离散格式得到的结果综合起来就是一个线性方程组如下：22221,11,(1)1,11,11,11,1(1),1(1),(1)0sin11sin1cos0sin1cos1nnnnnnnnaaaa也就是AVB，0A，令ANP,NA故AVPVBPV,即NVBPV1nnNVPVB其中N可以分为对角阵、三对角阵、上下三角阵。迭代式：111nnVNPVNB,给出初值0V或者1nnnnnNVPVBPVAVBAV,~nnRBAV残量推出1nnnNVNVR收敛准则：1）n足够大，残量趋于零（小于）2）10()nnVV4.3方程离散上述泊松方程在(,)ijxy的正五点差分格式为：1,,1,,1,,12222sincosijijijijijijijxyxy将泊松方程在计算域网格点(,)ijxy上进行差分得到差分方程：22223,22,21,22,32,22,12222223,32,31,32,42,32,22322223,n12,n11,n12,n2,n12,n22,2(2)(2)sincos2,3(2)(2)sincos2,1(2)(2)sijyxxyxyijyxxyxyijnyxxyMM2122224,23,22,23,33,23,13222224,33,32,33,43,33,233224,n13,n12,n13,n3,nincos3,2(2)(2)sincos3,3(2)(2)sincos3,1(2)(2nxyijyxxyxyijyxxyxyijnyxMM2213,n231)sincosnxyxyM2222m,2m1,2m2,2m1,3m1,2m1,1122222m,3m1,3m2,3m1,4m1,3m1,2132m,n1m1,n1m2,n11,2(2)(2)sincos1,3(2)(2)sincos1,1(2mmimjyxxyxyimjyxxyxyimjnyMM222m1,nm1,n1m1,n211)(2)sincosmnxxyxy以下对边界条件进行离散：1,m,k,1k,n1:0,(1,2...)sin1cos:,(1),(1,2...)2sin1:,(1),(1,2...)2sincos1:,(1),(1,2...)2kkkkkkkkkkiknyimyykyknxjxkxkmxjnxxkxkm左边界右边界下边界上边界综合差分方程及边界条件将方程整理成线性方程组KB的形式得到：(2)(n2)(2)(n2)mmGIIGIKIGIIGOOO其中，2(2)(2)11,=11mmIaaxO其中22(2)(2)2()2(),=,;2()2()mma+bbba+bbGaxbyba+bbba+bOOO其中23122,23,2m1,232,33,3m1,312,n13,n1m1,n12,23,2m1,22,3m1,32,n1m1,n1(m2)(n2)1TnnTLLLLLLL；231TnBbbbL222,11,2323,12121,1231,333313sincossincossincossincossincossincosmmmabxyababxyababxyaabxybabxybabxyMM212,n1,n1313,n221111,nsincossincos,=,;sincosnnnmnmabxyababxyabaxbyabxyaM其中4.4线性方程组迭代算法对于线性方程组AXb，可以利用以下迭代算法进行求解。4.4.1Jacobi迭代将系数矩阵分解ADLU，Jacobi迭代格式为：(1)()kkxBxf其中1(LU)BD，1fDb。4.4.2超松弛迭代法(SuccessiveOver-Relaxation)将系数矩阵分解ADLU，SOR迭代格式为：(1)()kkxBxf其中1()((1)D)BDwLwwU，1()fwDwLb，w为迭代松弛因子，当=1w时，松弛迭代法就是Gauss-Seidel迭代法；当1w时被称为超松弛迭代。4.4.3共轭梯度法(ConjugateGradient)共轭梯度法求解代数方程组AXb的算法[2]为：（1）假定初场(0)X；（2）计算初余量(0)(0)rbAX；（3）令(0)(0)=pr；（4）对0,1,2kL，直到(k)rtolerence（记(,)Tababg）(k)(k)(k)(k)=(,)/(,)krrApp(k1)(k)(k)kXXp(k1)(k)(k)krrAp=(k+1)(k+1)(k)(k)k(r,r)/(r,r)(k1)(k1)(k)kprp，结束。五、结果比较与分析由上面已经罗列的公式可知，精确求解的泊松方程为：1=-sincos2xyxy，利用MATLAB编写程序，步长取x=y=0.01分别用Jacobi迭代、超松弛迭代法、共轭梯度法进行计算，分别描出等值线：05.0,0.2,0.5,0.75,1，并与精确解对比，如下图所示：图1精确解的等值线图2Jacobi迭代法的等值线图3超松弛迭代法的等值线图4共轭梯度法的等值线综合来看，由图1、2、3、4可以看到，差分解与数值解基本相同。此外，给出整个求解区域的数值解和精确解的云图，如下所示：图5数值解云图图6精确解云图定义数值解与精确解之间的误差为1=*,*calEnm其中为精确解，计算得到：-41=*8.026110(Jacobi)calEnm迭代结果根据题意以迭代误差小于610为迭代收敛条件，将各个迭代算法求解整个问题的时间，求解线性方程组AXb的迭代步数，计算误差列表如下：JacobiSORCG直接求解耗费时间(s)0.785553.05530.04930.0321迭代步数74559961380精度8.0261E-045.6117E-052.6E-031.1457E-07其中：直接求解是MATLAB中求解线性方程组的左除算法，对于AXb，则\XAb，其结果相当于-1XAb；同时，SOR方法中的松弛因子1.75w取。综上所述：1.三种方法的数值解与精确解差别均很小，该算例验证了三种方法的正确性。2.根据题意采用的是均匀网格，从云图上我们可以看出，不同位置的等值线密度不同。因此，若采用非均匀网格，在等值线比较密的位置加密网格可以提高计算的精度。3.不同迭代算法效率不同，迭代步数最少的是CG方法，迭代步数最多的是Jacobi迭代，各个算法的迭代步数比较与理论相符；其中MATLAB自带求解线性方程组的左除算法效率和精度均最高，其次是共轭梯度法，Jacobi迭代，SOR算法；4.CG的迭代步数虽然较少，但是精度有待提高，可以通过提高算法收敛标准提高精度；5.程序采用的是正方形网格，未研究两个方向网格数目不同以及渐近网格对计算结果的影响。可以进一步研究其对计算结果和计算效率的影响。六、源程序及其主要符号说明符号说明如下表所示：lx,ly求解区间大小dx，dyx,y方向上的网格步长m，nx,y方向上的节点数K分块系数矩阵B非齐次项矩阵X解向量ps不考虑边界的解Ps考虑边界的解Ps_e精确解kk迭代次数Q计算精度源程序如下：1.主程序：clear%清空globalkk;kk=0;dx=0.01;dy=0.01;%设置网格步长lx=1;ly=1;%设置计算域大小m=lx/dx+1;n=ly/dy+1;%计算x、y方向节点数a=dx^2;b=dy^2;I=a*eye(m-2);%生成矩阵II=sparse(I);G=zeros(m-2,m-2);%生成矩阵GG(1,1:2)=[-2*(a+b)b];fori=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;fori=2:m-1forj=2:n-1Xi=(i-1)*dx;Yj=(j-1)*dy;B(k)=dx^2*dy^2*sin(Xi)*cos(Yj);ifj==2B(k)=B(k)+a*sin(Xi)/2;endifi==m-1B(k)=B(k)-b*(Yj-sin(1)*cos(Yj)/2);endifj==n-1B(k)=B(k)-a*(Xi-sin(Xi)*cos(1)/2);endk=k+1;endendK=cell(n-2);%生成分块系数矩阵KK(:)={zeros(m-2,m-2)}