方腔环流的流场计算

Page 5 of 10
中山大学工学院、理论与应用力学刘广编制
实验编号及题目：实验六
方腔环流
《计算流体力学试验》课程实验报告纸
院系：工学院
E0=0 !循环中进行判断，但是这种方法不好，可以采用-1，0，1标记法，循环区间为大区间，计算之前判断，-1 为外点，CYCLE掉，1为内点，用迭代公式，0为边界点，值不变 DO I=2,X,1 DO J=2,X,1 M(I,J)=(1-W)*TEMPM(I,J)+(W/4)*(M(I-1,J)+TEMPM(I+1,J)+M(I,J1)+TEMPM(I,J+1)+H*H*TEMPM1(I,J)) M1(I,J)=(W/4)*(TEMPM1(I+1,J)+M1(I-1,J)+TEMPM1(I,J+1)+M1(I,J-1)+(RE/4)*((TEMPM(I,J+1)M(I,J-1))*(TEMPM1(I+1,J)-M1(I-1,J))-(TEMPM(I+1,J)-M(I-1,J))*(TEMPM1(I,J+1)-M1(I,J-1))))+(1W)*TEMPM1(I,J) E1=ABS(TEMPM(I,J)-M(I,J)) E2=ABS(TEMPM1(I,J)-M1(I,J)) IF(MAX(E1,E2)>E0)THEN E0=MAX(E1,E2) END IF END DO END DO DO J=1,X+1!左边 M1(1,J)=-2*(TEMPM(2,J)-TEMPM(1,J))/(H*H) E2=ABS(TEMPM1(1,J)-M1(1,J)) E1=ABS(TEMPM(1,J)-M(1,J)) IF(MAX(E1,E2)>E0)THEN E0=MAX(E1,E2) END IF END DO DO J=1,X+1!右边 M1(X+1,J)=-2*(TEMPM(X,J)-TEMPM(X+1,J))/(H*H) E2=ABS(TEMPM1(X+1,J)-M1(X+1,J)) E1=ABS(TEMPM(X+1,J)-M(X+1,J)) IF(MAX(E1,E2)>E0)THEN E0=MAX(E1,E2) END IF END DO DO I=2,X!底边 M1(I,1)=-2*(TEMPM(I,2)-TEMPM(I,1))/(H*H)
中山大学工学院、理论与应用力学刘广编制
实验编号及题目：实验六
方腔环流
《计算流体力学试验》课程实验报告纸
院系：工学院 姓名：刘广 学号：11309018 日期：2014 年 05 月 23 号
方腔环流示意图
我们取网格做如下说明，取正方形网格，网格数量由用户输入，但是建议不超过 200，因 为超过 200 计算量会急剧增大，雷诺数也是由用户输入，建议不超过 5000，因为超过 5000 流 场已经呈现出非线性，开始出现偏差，雷诺数超过 10000 甚至根本不能算出结果，这是因为强 非线性情况下描述流场行为的方程已经不能做如上简化， 同样的， 在计算过程中我们采用超松 弛迭代法，下面对离散形式做一下说明： 我们用二阶精度的差商代替微商做以下说明，得
( i 1, j 2 i , j i 1, j ) / h 2 ( i , j 1 2 i , j i , j 1 ) / h 2 i , j ( i , j 1 i , j 1 )(
2 i 1, j
i 1, j )
4h 4h 2 1 i 1, j 2 i , j i 1, j i , j 1 2 i , j i , j 1 ( ) Re h2 h2
松弛因子 代入， i , j
n 1 2 n n n n ( in in, j , j h i 1, j i 1, j i , j 1 i , j 1 ) / 4 (1 )

( i 1, j i 1, j )(
i , j 1
院系：工学院 实验人姓名：刘广 日期：2014 年 05 月 23 号 专业：理论与应用力学 参加人姓名：刘广 温度：18° C 学号：11309018 年级：2011 级
实验六：方腔环流 一、 实验目的 1、在正方形水腔当中布满液体，水腔三方为墙壁，上方为盖板，盖板以速度 1 向右匀速 运动，当水腔内部液体达到平衡时，计算水箱内的流场。 2、计算全区域的流函数分布，并绘制流场图像。 二：实验仪器与设备: ① 装有 Fortran 的计算机 三、实验原理 1.方程与边界条件 用流函数涡量法， 和 满足下列无量纲定常方程与边界条件 1台

Re n n [( in, j 1 in, j 1 )( in in1, j in1, j )( in 1, j i 1, j ) ( , j 1 i , j 1 )]} 4 (1 ) in ,j
边界值的点需要特殊处理，然后由内点差分格式顺序扫描，一次次迭代直到满足用户输入的 精度为止。边界条件如下：
姓名：刘广
学号：11309018
日期：2014 年 05 月 23 号
Page 6 of 10
中山大学工学院、理论与应用力学刘广编制
实验编号及题目：实验六
方腔环流
《计算流体力学试验》课程实验报告纸
院系：工学院
E2=ABS(TEMPM1(I,1)-M1(I,1)) E1=ABS(TEMPM(I,1)-M(I,1)) IF(MAX(E1,E2)>E0)THEN E0=MAX(E1,E2) END IF END DO DO I=2,X M1(I,X+1)=-2*(TEMPM(I,X)-TEMPM(I,X+1)+H)/(H*H) E2=ABS(TEMPM1(I,X+1)-M1(I,X+1)) E1=ABS(TEMPM(I,X+1)-M(I,X+1)) IF(MAX(E1,E2)>E0)THEN E0=MAX(E1,E2) END IF END DO !层之间交换 TEMPM(:,:)=M(:,:) TEMPM1(:,:)=M1(:,:) END DO DO I=1,X+1,1 DO J=1,X+1,1 WRITE(8,*)I,',',J,',',M(I,J) END DO END DO DO I=1,X+1,1 DO J=1,X+1,1 WRITE(10,*)I,',',J,',',M1(I,J) END DO END DO END SUBROUTINE
PBiblioteka Baiduge 7 of 10
中山大学工学院、理论与应用力学刘广编制
实验编号及题目：实验六
方腔环流
《计算流体力学试验》课程实验报告纸
院系：工学院 姓名：刘广 学号：11309018 日期：2014 年 05 月 23 号
我们写如下MATLAB代码：b=importdata('OUTPUT.txt');
SUBROUTINE StreamFunctionVorticity() INTEGER::I,J,X !迭代精度，松弛因子,上层误差，最终误差 REAL::FACTOR,W,E1,E2,E0,RE,H !声明动态矩阵，用于存储计算结果，维数为I，J，M为流函数，TEMPM为涡量，边界两者相等 REAL,ALLOCATABLE::M(:,:),TEMPM(:,:),M1(:,:),TEMPM1(:,:) !打开通道号8，输出OUTPUT_C.CSV OPEN(UNIT=8,FILE='OUTPUTLiu.txt') OPEN(UNIT=10,FILE='OUTPUTWo.txt') WRITE(*,*)'请输入横向网格数（纵向相同；注意，大于200计算量会非常大）' READ(*,*)X !确定M大小，开始在网格第一层赋初值，其中循环不能从负数开始，X网格右移 ALLOCATE(M(X+1,X+1),M1(X+1,X+1)) ALLOCATE(TEMPM(X+1,X+1),TEMPM1(X+1,X+1)) WRITE(*,*)'请输入雷诺数Re（10，100或1000；注意，雷诺数不能超过1万）' READ(*,*)RE WRITE(*,*)'请输入松弛因子(RE=10为1，100为0.7，1000为0.5)' READ(*,*)W H=0.001 M(:,:)=0 TEMPM(:,:)=0 M1(:,:)=0 TEMPM1(:,:)=0 WRITE(*,*)'请输入迭代精度Factor(小于0.001计算量会急剧增大)' READ(*,*)FACTOR E0=FACTOR !开始计算内点，采用超松弛因子法计算，松弛因子取1.7269，迭代精度factor用户输入 DO WHILE(E0>=FACTOR) !误差规零
sn 1 0, 在四边
2( sn * sn ) , 在两侧和底边 h2 2( sn * sn h) n 1 s , 在顶边 h2
n 1 s
三：实验步骤: 学会对MS-FORTRAN的基本操作。 用Fortran编写程序计算方腔环流流场计算源代码，我们可 以写如下Fortran90程序： 主程序如下：
Page 4 of 10
中山大学工学院、理论与应用力学刘广编制
实验编号及题目：实验六
方腔环流
《计算流体力学试验》课程实验报告纸
院系：工学院
END IF WRITE(*,*)'请在目录下查看结果' PAUSE STOP END
姓名：刘广
学号：11309018
日期：2014 年 05 月 23 号
以下为流函数涡量法的子程序，
! Author : GUANG_LIU * owenyaa@gmail.com * ! Created Time : X+14-04-22 19:50 ! Last Revised : GUANG_LIU ,X+14-04-22 ! Remark : 以下实验结果中均为网格19×19，步长0.001，松弛因子0.5，迭代精度0.000001。 !Re=100 PROGRAM MAIN IMPLICIT NONE INTEGER::I WRITE(*,*)'请选择计算方法,流函数涡量法输入1,速度压力法输入2' READ(*,*)I IF (I==1)THEN CALL StreamFunctionVorticity() ELSE IF (I==2)THEN CALL PressureVelocity() Re=1000
中山大学工学院计算流体力学实验报告
实验名称：方腔环流的流场计算 姓 名： 刘广 刘广 11309018 詹杰民
参与组员： 学 号：
任课教师：
学科专业：工学院 理论与应用力学
中山大学 2014 年 05 月 23 日
中山大学工学院、理论与应用力学刘广编制
实验编号及题目：实验六
方腔环流
《计算流体力学试验》课程实验报告纸
2 2 x 2 x 2 1 2 y x x y Re
（1）
（2）
0 ， 在腔体四边
v 0 ，在 AD 和 BC 上 x
u 0 ，在 AB 和 DC 上 y
其中，Re 为雷诺数。
Page 2 of 10
姓名：刘广
学号：11309018
日期：2014 年 05 月 23 号
程序执行时，我输入的是10E-3作为计算精度，对于程序的说明程序注释已经说的很清楚， 这里不再赘述。最终会在目录下生成OUTPUT.txt 。 四：实验结果 最终会在目录下生成OUTPUT.txt 。然后我们可以用Matlab软件对计算出来的数据进一步整 理，画成云图和画出流线，我们做一下说明。
i , j 1 )
Page 3 of 10
中山大学工学院、理论与应用力学刘广编制
实验编号及题目：实验六
方腔环流
《计算流体力学试验》课程实验报告纸
院系：工学院
1 in ,j
姓名：刘广
学号：11309018
日期：2014 年 05 月 23 号

4
n n n { in 1, j i 1, j i , j 1 i , j 1