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

CFD 实验报告二

姓名： 学号：

一、题目

求解Poisson 方程

y x y x c o s s i n 2222=??+??ψ

ψ，

10≤≤x ， 10≤≤y ， 0|0==x ψ，

y x ==1|ψ2

cos 1sin y

-

， ==0|y ψ2sin x -， x y ==1|ψ2

1

c o s s i n x -， 描出等值线：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 等迭代方法求解线性方程组，即可得到函数,i j ψ的在网

格点处的离散值。同时，计算域为Lx Ly ?的矩形区域，划分结构网格，均匀步长，设x 方向网格步长x ?，y 方向网格步长y ?，则x 方向网格点数为m=/1Lx x ?+， y 方向网格点数为n=/1Ly y ?+。

四、计算公式和步骤

4.1 精确解的计算

题中已知四个边界条件，可以通过现已有的解析求解不难求得精确解（解析

解）为1

sin cos 2

x y xy ψ=-?+。

4.2迭代求解一般过程

迭代法是求解离散代数方程组的主要方法。假如我们对题目中的PDE 给定一个离散格式，则对每一个确定的离散点（i ，j ），PDE 转化为FDE ，为

,1,,(,sin[(1)]cos[(1)])i j a b c d F i x i y ψψψ=-?-? ，其中,,~a b c d ψψ是离散格式中所有

与,i j ψ有关的点，函数1F 中所有元素都是线性叠加，即1F 是线性多元函数。所以，将边界上给定结果以及中间的离散格式得到的结果综合起来就是一个线性方程组如下：

22221,11,(1)1,11,11,11,1(1),1(1),(1)0sin11sin1cos 0sin1cos1n n n n n n n n a a a a ψψψψ++++++++??????

?????????????????????????

=??????????????????????????????????????

也就是AV B =，0

A ≠，令A N P =-,N A ≈ 故AV PV

B PV +=+,即NV B PV =+

1n n NV PV B +∴=+

其中N 可以分为对角阵、三对角阵、上下三角阵。 迭代式：

111n n V N PV N B +--=+,给出初值0V

或者1n n n n n NV PV B PV AV B AV +=+=++-,~n n R B AV =-残量

推出1n n n NV NV R +=+

收敛准则：1）n 足够大，残量趋于零（小于ε） 2）10()n n V V ε+-→< 4.3方程离散

上述泊松方程在(,)i j x y 的正五点差分格式为：

1,,1,,1,,1

2

2

22sin cos i j i j i j i j i j i j i j x y x

y

ψψψψψψ+-+--+-++

=??

将泊松方程在计算域网格点(,)i j x y 上进行差分得到差分方程：

22223,22,21,22,32,22,12222223,32,31,32,42,32,22322223,n 12,n 11,n 12,n 2,n 12,n 22,2(2)(2)sin cos 2,3

(2)(2)sin cos 2,1(2)(2)s i j y x x y x y i j y x x y x y i j n y x x y ψψψψψψψψψψψψψψψψψψ-----==?-++?-+=??==?-++?-+=??==-?-++?-+=??M M 2122224,23,22,23,33,23,13222224,33,32,33,43,33,233224,n 13,n 12,n 13,n 3,n in cos 3,2(2)(2)sin cos 3,3

(2)(2)sin cos 3,1(2)(2n x y i j y x x y x y i j y x x y x y i j n y x ψψψψψψψψψψψψψψψψψ----??

??

?

??

==?-++?-+=??==?-++?-+=??==-?-++?-M M 2213,n 231)sin cos n x y x y ψ---??

??

?

?+=???

M

2222m,2m 1,2m 2,2m 1,3m 1,2m 1,112

2222m,3m 1,3m 2,3m 1,4m 1,3m 1,2132m,n 1m 1,n 1m 2,n 11,2(2)(2)sin cos 1,3(2)(2)sin cos 1,1(2m m i m j y x x y x y i m j y x x y x y i m j n y ψψψψψψψψψψψψψψψ-----------------=-=?-++?-+=??=-=?-++?-+=??=-=-?-+M M 222m 1,n m 1,n 1m 1,n 211)(2)sin cos m n x x y x y ψψψ-------????

?

?+?-+=???

以下对边界条件进行离散：

1,m,k,1k,n 1:0,(1,2...)sin1cos :,(1),(1,2...)2sin 1:,(1),(1,2...)2sin cos1:,(1),(1,2...)

2

k k k k k k

k k k k i k n y i m y y k y k n x j x k x k m x j n x x k x k m ψψψψ===?

?

?==-=-?=???==-=-?=??

?==-=-?=??左边界右边界下边界上边界

综合差分方程及边界条件将方程整理成线性方程组K B Φ=的形式得到：

(2)(n 2)(2)(n 2)

m m G I I G I K I G I I G --?--??

??????=?

??

?????O O

O

其中，

2(2)(2)

11,=11m m I a a x -?-??

??????=????????

?O

其中22(2)(2)

2()

2(),=,;2()2()m m a+b b b a+b b G a x b y b a+b b b a+b -?--??

??-????=?=??

?

-?

?

??-??O O

O

其中 []

23122,23,2m 1,232,33,3m 1,312,n 13,n 1m 1,n 12,23,2m 1,22,3m 1,32,n 1m 1,n 1(m 2)(n 2)1

T

n n T

φφφφψψψφψψψφψψψψψψψψψψ---------------?Φ=??=????=????=??

??Φ=??

L

L L L

L L

L

；

[]2

31T

n B b b b -=L

222,11,2323,12121,1231,333313sin cos sin cos sin cos sin cos sin cos sin cos m m m ab x y a b ab x y a b ab x y a ab x y b ab x y b ab x y ψψψψψ-----??

??

-?

?=????-??-????

?

?=??????

M M

212,n 1,n 1313,n 221111,n sin cos sin cos ,=,;sin cos n n n m n m ab x y a b ab x y a b a x b y ab x y a ψψψψ---------????-??=?=?????-??

M 其中

4.4线性方程组迭代算法

对于线性方程组AX b =，可以利用以下迭代算法进行求解。

4.4.1 Jacobi 迭代

将系数矩阵分解A D L U =--，Jacobi 迭代格式为：

(1)()k k x Bx f +=+

其中1(L U)B D -=+，1f D b -=。

4.4.2超松弛迭代法(Successive Over-Relaxation)

将系数矩阵分解A D L U =--，SOR 迭代格式为：

(1)()k k x Bx f +=+

其中1()((1)D )B D wL w wU -=--+，1()f w D wL b -=-，w 为迭代松弛因子，当=1w 时，松弛迭代法就是Gauss-Seidel 迭代法；当1w >时被称为超松弛迭代。

4.4.3共轭梯度法(Conjugate Gradient)

共轭梯度法求解代数方程组AX b =的算法[2]为： （1）假定初场(0)X ；

（2）计算初余量(0)(0)r b AX =-； （3）令(0)(0)=p r ；

（4）对0,1,2k =L ，直到(k)r tolerence <（记(,)T a b a b =g ）

(k)(k)(k)(k)=(,)/(,)k r r Ap p α

(k 1)(k)(k)k X X p α+=+ (k 1)(k)(k)k r r Ap α+=-

=(k+1)(k+1)(k)(k)k (r ,r )/(r ,r )β

(k 1)

(k

1)k p r p β++=- ，结束。

五、结果比较与分析

由上面已经罗列的公式可知，精确求解的泊松方程为：

1

=-sin cos 2

x y xy ψ+，

利用MATLAB 编写程序，步长取x ?=y ?=0.01分别用Jacobi 迭代、超松弛迭代法、共轭梯度法进行计算，分别描出等值线：05.0=ψ, 0.2, 0.5, 0.75, 1，并与精确解对比，如下图所示：

图1 精确解的等值线

图2 Jacobi 迭代法的等值线

图3超松弛迭代法的等值线

图4共轭梯度法的等值线

综合来看，由图1、2、3、4可以看到，差分解与数值解基本相同。

此外，给出整个求解区域的数值解和精确解的云图，如下所示：

图5 数值解云图

图6精确解云图

定义数值解与精确解之间的误差为1

=*,*cal E nm

ψψψ-∑其中为精确解，计算得到：

-41

=

*8.026110(Jacobi )cal E nm

ψψ-=?∑迭代结果 根据题意以迭代误差小于610-为迭代收敛条件，将各个迭代算法求解整个问题的时间，求解线性方程组AX b =的迭代步数，计算误差列表如下：

其中：直接求解是MATLAB 中求解线性方程组的左除算法，对于AX b =，则\X A b =，其结果相当于-1X A b =；同时，SOR 方法中的松弛因子 1.75w 取。

综上所述：

1. 三种方法的数值解与精确解差别均很小，该算例验证了三种方法的正确性。

2. 根据题意采用的是均匀网格，从云图上我们可以看出，不同位置的等值线密度不同。因此，若采用非均匀网格，在等值线比较密的位置加密网格可以提高计算的精度。

3. 不同迭代算法效率不同，迭代步数最少的是CG 方法，迭代步数最多的是Jacobi 迭代，各个算法的迭代步数比较与理论相符；其中MATLAB 自带求解线性方程组的左除算法效率和精度均最高，其次是共轭梯度法，Jacobi 迭代，SOR 算法；

4. CG 的迭代步数虽然较少，但是精度有待提高，可以通过提高算法收敛标准提高精度；

5. 程序采用的是正方形网格，未研究两个方向网格数目不同以及渐近网格对计算结果的影响。可以进一步研究其对计算结果和计算效率的影响。

六、源程序及其主要符号说明

符号说明如下表所示：

源程序如下：

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=dx^2;b=dy^2;

I=a*eye(m-2);%生成矩阵I

I=sparse(I);

G=zeros(m-2,m-2);%生成矩阵G

G(1,1:2)=[-2*(a+b) b];

for i=2:m-3

G(i,i-1:i+1)=[b -2*(a+b) b];

end

G(m-2,m-3:m-2)=[b -2*(a+b)];

G=sparse(G);

M=(m-2)*(n-2);%生成矩阵B

B=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)=dx^2*dy^2*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;

end

end

K=cell(n-2);%生成分块系数矩阵K

K(:)={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};

end

K(n-2,n-3:n-2)={I G};

K=cell2mat(K);

disp('~~~~~~~~~~~~~~~~~~~SC16005041李文建~~~~~~~~~~~~~~~~~~~~~~~~'); disp(' ');

tic %开始计时并迭代求解线性方程组

% X=K\B; %左除算法

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)'};

end

ps=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.0f\n', kk);

disp('计算精度为：');

Q

contour(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);

end

3. 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<1

f=D\b;

X=B*X0+f;

i=1;

while norm(X-X0,inf)>=1e-6

X0=X;

X=B*X0+f;

i=i+1;

if i==1E6

disp('迭代步数太多，放弃计算！')

break;

end

end

kk=i; %输出迭代步数

else

disp('不收敛！')

X=[];

end

end

4. 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<1

f=w*((D-w*L)\b);

X=B*X0+f;

i=1;

while norm(X-X0,inf)>=1e-6

X0=X;

X=B*X0+f;

i=i+1;

if i==1E6

disp('迭代步数太多，放弃计算！')

break;

end

end

kk=i; %输出迭代步数

else

disp('不收敛！')

X=[];

end

5. 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;

end

X=x;

kk=k; %迭代次数

end

七、参考文献

[1] 吴清松. 计算热物理引论[M]. 中国科学技术大学出版社. 合肥.2009

[2] 孔凡铀. 用FFT法与SOR法求解二维泊松方程的比较[J]. 气象科学. 1988年第一期

[3] 杨金凤,邓居智,陈辉. 多重网格法在求解泊松方程中的应用进展[J]. 内蒙古石油化工. 2011年第24期

[4] 王忆锋, 唐利斌. 利用有限差分和MATLAB矩阵运算直接求解二维泊松方程[J]. 红外技术. 第32卷第4期. 2010.04