方腔环流的流场计算

中山大学工学院计算流体力学实验报告实验名称：方腔环流的流场计算姓名：刘广参与组员：刘广学号：11309018任课教师：詹杰民学科专业：工学院理论与应用力学中山大学2014年05月23日中山大学工学院、理论与应用力学刘广编制实验编号及题目：实验六方腔环流《计算流体力学试验》课程实验报告纸院系：工学院专业：理论与应用力学年级：2011级实验人姓名：刘广参加人姓名：刘广日期：2014年05月23号温度：18°C学号：11309018实验六：方腔环流一、实验目的1、在正方形水腔当中布满液体，水腔三方为墙壁，上方为盖板，盖板以速度1向右匀速运动，当水腔内部液体达到平衡时，计算水箱内的流场。2、计算全区域的流函数分布，并绘制流场图像。二：实验仪器与设备:①装有Fortran的计算机1台三、实验原理1.方程与边界条件用流函数涡量法，和满足下列无量纲定常方程与边界条件2222xx（1）2Re1yxxy（2）0，在腔体四边0vx，在AD和BC上0uy，在AB和DC上其中，Re为雷诺数。Page2of10中山大学工学院、理论与应用力学刘广编制实验编号及题目：实验六方腔环流《计算流体力学试验》课程实验报告纸院系：工学院姓名：刘广学号：11309018日期：2014年05月23号方腔环流示意图我们取网格做如下说明，取正方形网格，网格数量由用户输入，但是建议不超过200，因为超过200计算量会急剧增大，雷诺数也是由用户输入，建议不超过5000，因为超过5000流场已经呈现出非线性，开始出现偏差，雷诺数超过10000甚至根本不能算出结果，这是因为强非线性情况下描述流场行为的方程已经不能做如上简化，同样的，在计算过程中我们采用超松弛迭代法，下面对离散形式做一下说明：我们用二阶精度的差商代替微商做以下说明，得221,,1,,1,,1,,1,11,1,1,1,,1,1221,,1,,1,,122(2)/(2)/()()()()44221()Reijijijijijijijijijijijijijijijijijijijijijhhhhhh松弛因子代入，12,,1,1,,1,1,()/4(1)nnnnnnnijijijijijijijhPage3of10中山大学工学院、理论与应用力学刘广编制实验编号及题目：实验六方腔环流《计算流体力学试验》课程实验报告纸院系：工学院姓名：刘广学号：11309018日期：2014年05月23号1,1,1,,1,1,1,11,1,1,1,,1,1,{4Re[()()()()]}4(1)nnnnnijijijijijnnnnnnnnijijijijijijijijnij边界值的点需要特殊处理，然后由内点差分格式顺序扫描，一次次迭代直到满足用户输入的精度为止。边界条件如下：112120,2(*),2(*),nsnnnsssnnnssshhh在四边在两侧和底边在顶边三：实验步骤:学会对MS-FORTRAN的基本操作。用Fortran编写程序计算方腔环流流场计算源代码，我们可以写如下Fortran90程序：主程序如下：!Author:GUANG_LIU*[e-mail]*!CreatedTime:X+14-04-2219:50!LastRevised:GUANG_LIU,X+14-04-22!Remark:以下实验结果中均为网格19×19，步长0.001，松弛因子0.5，迭代精度0.000001。!Re=100Re=1000PROGRAMMAINIMPLICITNONEINTEGER::IWRITE(*,*)'请选择计算方法,流函数涡量法输入1,速度压力法输入2'READ(*,*)IIF(I==1)THENCALLStreamFunctionVorticity()ELSEIF(I==2)THENCALLPressureVelocity()Page4of10中山大学工学院、理论与应用力学刘广编制实验编号及题目：实验六方腔环流《计算流体力学试验》课程实验报告纸院系：工学院姓名：刘广学号：11309018日期：2014年05月23号ENDIFWRITE(*,*)'请在目录下查看结果'PAUSESTOPEND以下为流函数涡量法的子程序，SUBROUTINEStreamFunctionVorticity()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.CSVOPEN(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(*,*)REWRITE(*,*)'请输入松弛因子(RE=10为1，100为0.7，1000为0.5)'READ(*,*)WH=0.001M(:,:)=0TEMPM(:,:)=0M1(:,:)=0TEMPM1(:,:)=0WRITE(*,*)'请输入迭代精度Factor(小于0.001计算量会急剧增大)'READ(*,*)FACTORE0=FACTOR!开始计算内点，采用超松弛因子法计算，松弛因子取1.7269，迭代精度factor用户输入DOWHILE(E0=FACTOR)!误差规零Page5of10中山大学工学院、理论与应用力学刘广编制实验编号及题目：实验六方腔环流《计算流体力学试验》课程实验报告纸院系：工学院姓名：刘广学号：11309018日期：2014年05月23号E0=0!循环中进行判断，但是这种方法不好，可以采用-1，0，1标记法，循环区间为大区间，计算之前判断，-1为外点，CYCLE掉，1为内点，用迭代公式，0为边界点，值不变DOI=2,X,1DOJ=2,X,1M(I,J)=(1-W)*TEMPM(I,J)+(W/4)*(M(I-1,J)+TEMPM(I+1,J)+M(I,J-1)+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))))+(1-W)*TEMPM1(I,J)E1=ABS(TEMPM(I,J)-M(I,J))E2=ABS(TEMPM1(I,J)-M1(I,J))IF(MAX(E1,E2)E0)THENE0=MAX(E1,E2)ENDIFENDDOENDDODOJ=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)THENE0=MAX(E1,E2)ENDIFENDDODOJ=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)THENE0=MAX(E1,E2)ENDIFENDDODOI=2,X!底边M1(I,1)=-2*(TEMPM(I,2)-TEMPM(I,1))/(H*H)Page6of10中山大学工学院、理论与应用力学刘广编制实验编号及题目：实验六方腔环流《计算流体力学试验》课程实验报告纸院系：工学院姓名：刘广学号：11309018日期：2014年05月23号E2=ABS(TEMPM1(I,1)-M1(I,1))E1=ABS(TEMPM(I,1)-M(I,1))IF(MAX(E1,E2)E0)THENE0=MAX(E1,E2)ENDIFENDDODOI=2,XM1(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)THENE0=MAX(E1,E2)ENDIFENDDO!层之间交换TEMPM(:,:)=M(:,:)TEMPM1(:,:)=M1(:,:)ENDDODOI=1,X+1,1DOJ=1,X+1,1WRITE(8,*)I,',',J,',',M(I,J)ENDDOENDDODOI=1,X+1,1DOJ=1,X+1,1WRITE(10,*)I,',',J,',',M1(I,J)ENDDOENDDOENDSUBROUTINE程序执行时，我输入的是10E-3作为计算精度，对于程序的说明程序注释已经说的很清楚，这里不再赘述。最终会在目录下生成OUTPUT.txt。四：实验结果最终会在目录下生成OUTPUT.txt。然后我们可以用Matlab软件对计算出来的数据进一步整理，画成云图和画出流线，我们做一下说明。Page7of10中山大学工学院、理论与应用力学刘广编制实验编号及题目：实验六方腔环流《计算流体力学试验》课程实验报告纸院系：工学院姓名：刘广学号：11309018日期：2014年05月23号我们写如下MATLAB代码：b=importdata('OUTPUT.txt');x=b(:,1);y=b(:,2);z=b(:,3);%对X、Y轴做网格划分xtemp=linspace(min((x-1)*0.25),max((x-1)*0.25),80);ytemp=linspace(min((y-1)*0.25),max((y-1)*0.25),50);[X,Y]=meshgrid(xtemp,ytemp);Z=griddata((x-1)*0.25,(y-1)*0.25,z,X,Y);%画出图像，设定视角，设定座标区间%holdoncontour3(X,Y,Z,500)xlabel('x轴');ylabel('y轴')zlabel('二维流场计算')执行后生成对应云图如下图1.以下我们均取网格为100*100，雷诺数依次为100，1000，2000，5000.迭代精度为0.001；画出对应图像如下图：图1雷诺数为100Page8of10中山大学工学院、理论与应用力学刘广编制实验编号及题目：实验六方腔环流《计算流体力学试验》课程实验报告纸院系：工学院姓名：刘广学号：11309018日期：2014年05月23号图2雷诺数为1000图3雷诺数为2000Page9of10中山大学工学院、理论与应用力学刘广编制实验编号及题目：实验六方腔环流《计算流体力学试验》课程实验报告纸院系：工学院姓名：刘广学号：11309018日期：2014年0