计算流体力学上机实验报告

《计算流体力学》

上机实验报告

班级：

姓名：

学号：

北京航空航天大学流体力学研究所

上机实验名称

两平行平板间不可压缩流体绕物体的平面无旋流动

一、实验目的

通过具体算例，熟悉和掌握使用CFD方法获取给定流场的流动参数。

二、实验内容、方法及步骤

1.流动问题描述

如下图所示，考虑在平行放置的两平板之间流过的理想不可压缩流体绕方形物体的平面无旋流动。

2.求解区域

H=；绕流物体是边长为2的正方形，设两平行平板之间的距离为6

L=。根据流动的对称性，可取流上游来流入口位置与物体中心的距离为3

动区域的四分之一作为求解区域W，如下图所示。

3. 控制方程

对于不可压缩流体的平面无旋流动，流函数 y 在区域 W 内满足Laplace 方程

22

2

20x

y

y

y

抖+=抖

4. 边界条件

（1）OABC 是一条流线，规定

0OA BC

y

= ；

（2）对 OE 上任意一点 ()0,P y ，有

P y y = ；

（3） ED 也是一条流线，所以

2ED

H

y =

； （4）根据对称性，在 CD 上有

0CD

x

y ?=? 。

5. 定解问题

对于这里考虑的流动，可用下述定解问题来描述

22

2

20 , 0 , , , 2 0 , x

y

OA BC y OE H ED CD x

y

y

y y y y 抖+=W 抖===?=?在 内在上在上在上

在 上

6. 求解区域的离散化 － 计算网格

将单位长度等分成n 份，记1

h n

=

，于是求解区域W 沿x 方向可划分成M L n =?个网格，用0,1,2,3,,j M =L 来标记；沿y 方

向可划分成2

H

N n =?个网格，用0,1,2,3,,k N =L 来标记。 这些网格点可分成三类：

（1）当 ()01j L n <<- 且 0k N << ，

或者当 ()1L n j M -<< 且 n k N << 时， 网格点落在流场 W 内部，称为内点。这些网格点上的流函数需通过求解方程组来计算； （2）当 ()1L n j M -

n ? 时，

网格点落在正方形物体内部，网格点上不存在流场，无需计算；

（3）其余的网格点落在流场 W 的边界上，称为边界点。这些

网格点上的流函数直接由边界条件给定，也无需计算。

7. 定解问题的离散化 － 差分格式

Laplace 方程 22

2

20x

y

y

y

抖+=抖 的差分近似为 1,,1,,1,,1

2

2

220j k j k j k

j k j k j k h

h

y y y y y y +-+--+-++

=

边界条件

0x

?=?y

的差分近似为 ,1,0j k j k

h

y y --=

8. 内点上流函数的计算－ 迭代算法

在实际的计算中，内点的数量非常多，计算流函数需要求解大型的代数方程组。将Laplace 方程的差分近似改写成

(),1,1,,1,11

4

j k j k j k j k j k y y y y y +-+-=

+++ 使用高斯-赛德尔迭代算法计算：

()

(1)

()(1)()(1),1,1,,1,1

14

m m m m m j k

j k j k j k j k y

y y y y ++++-+-=+++ （ 0,1,2,3,m =L ） 给定迭代的初始近似 (0),1j k

y = 和迭代的收敛精度 δ15e =- 。

对 0,1,2,3,m =L ：

在流场内的每一个内点上用迭代公式

()

(1)()(1)()(1),1,1,,1,1

14

m m m m m j k

j k j k j k j k y

y y y y ++++-+-=+++ 进行计算。等式右边如果有边界点，可直接将边界条件指定的流函数值代入。

在边界 CD 上的每一个网格点（不包括端点C 、D ）处，

用公式 (1)(1)

,1,m m j k j k

++-=y y 计算。 对所有内点，计算 (1)()

,,,max m m j k j k

j k

+-y y 。 若 (1)()

,,,max m

m j k j k

j k

δ+-

否则将m加1，返回步骤3

9.计算源程序（matlab 编制）

function [liuxian,vx,vy]=jisuanliuti(x,y)

%%%%%%%%%%%%%%%%%%this is for Computational Fluid Dynamics;

%%%%%%%%%%%%%%%%使用方式：

[liuxian,vx,vy]=jisuanliuti(x,y)；

%%%%%%%%%%%%%%%输入x:x轴方向单位长度上划分的网格（x>=3）

%%%%%%%%%%%%%%%输入y:y方向上单位长度上划分的网格(y>=3)

%%%%%%%%%%%%%%输出liuxian:输出各个网格点的流函数值

%%%%%%%%%%%%%%输出vx：输出x方向的速度u

%%%%%%%%%%%%%%输出vy：输出y方向的速度v

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% if (x<3)||(y<3)

error('number of grids can not be divided too little');

end

%%%%%%%%%%%%%%%%参数定义

%%%%%%%%%%%%%%%%%%%%%%%%%%%%

n=x;%x方向单位长度划分n个

h=y;%y方向单位长度划分h个

H=3;%外边界的流函数值

fai=ones(3*n+1,3*h+1);%流函数%%%%%%%%%%%%%%边界条件%%%%%%%%%%%%%% fai(:,1)=0;

fai(2*n+1:3*n+1,1:h+1)=0;

fai(1,:)=((1:3*h+1)-1)*3/(3*h);

fai(:,end)=H;

fai(end,:)=fai(end-1,:); %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%计算流函数，用的是高斯-赛德尔迭代%%%%%%%%%%%%%%%%%%%%%%%%%%%%

fai1=fai;

for i=1:1000

for j=2:1:3*n

for k=2:1:3*h

t=(j<2*n+1)||((j>=2*n+1)&&(k>=h+2));

if t==1

fai1(j,k)=1/4*(fai1(j+1,k)+fai1(j-1,k)+fai1(j,k+1)+fai1(j,k-1));

end

end

end

c=max(max(abs(fai1-fai)));

if c>1e-5

fai=fai1;

else break;

end

fai1(end,:)=fai1(end-1,:);

end

liuxian=fai'; %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%画等高线图%%%%%%%%%%%%%%%%% figure;

contour(0:1/n:3,0:1/h:3,liuxian,8,'LineWidth',2);

title('流线图','fontsize',14); %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%改变流函数矩阵格式%%%%%%%%%%%% for i=1:length(liuxian(:,1))

liuxian1(i,:)=liuxian(length(liuxian(:,1))+1-i,:);

end

liuxian=liuxian1; %%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%速度的计算%%%%%%%%%%%%%% %%%%%%%%%%%x方向速度%%%%%%%%%%%%%%%% vx=(liuxian(1:end-1,:)-liuxian(2:end,:))*h;

%%%%%%%%%%%y 方向速度%%%%%%%%%%%%%%% vy=-(liuxian(:,2:end)-liuxian(:,1:end-1))*n; end

三、 实验结果

x

y

流线图

边界CD 上的速度u 的值如下表：

U 1.6765 1.656 1.615 1.558 1.502 Y 1 1.125 1.375 1.625 1.875 U 1.456 1.423 1.401 1.39 1.3845 Y 2.125

2.375

2.625

2.875

3

四、 讨论

1. x,y 方向划分不同的网格计算精度有所不同，划分网格越多，精度越高。下图是分别按（3,3）和（5,5）划分网格时的流线图。

流线图

流线图

流线图

Figure 1网格划分为（5，5）

从两幅图处可以清晰地看到，网格数较多的得到的流线更加光滑。 2. 网格形状的影响

分别考虑了（3,5）和（5,3）两种网格划分方式，对比如下

可见那个方向划分的网格数较多，那个方向的计算更加精细。

流线图

Figure 3网格划分为（3，5）

Figure 4网格划分为（5，3）

Figure 2网格划分为（3，3）