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

《计算流体力学》上机实验报告班级：姓名：学号：北京航空航天大学流体力学研究所1上机实验名称两平行平板间不可压缩流体绕物体的平面无旋流动一、实验目的通过具体算例，熟悉和掌握使用CFD方法获取给定流场的流动参数。二、实验内容、方法及步骤1.流动问题描述如下图所示，考虑在平行放置的两平板之间流过的理想不可压缩流体绕方形物体的平面无旋流动。2.求解区域设两平行平板之间的距离为6H=；绕流物体是边长为2的正方形，上游来流入口位置与物体中心的距离为3L=。根据流动的对称性，可取流动区域的四分之一作为求解区域W，如下图所示。23.控制方程对于不可压缩流体的平面无旋流动，流函数y在区域W内满足Laplace方程22220xyyy抖+=抖4.边界条件（1）OABC是一条流线，规定0OABCy=；（2）对OE上任意一点()0,Py，有Pyy=；3（3）ED也是一条流线，所以2EDHy=；（4）根据对称性，在CD上有0CDxy¶=¶。5.定解问题对于这里考虑的流动，可用下述定解问题来描述22220,0,,,20,xyOABCyOEHEDCDxyyyyyy抖+=W抖===¶=¶在内在上在上在上在上6.求解区域的离散化－计算网格将单位长度等分成n份，记1hn=，于是求解区域W沿x方向可划分成MLn=?个网格，用0,1,2,3,,jM=L来标记；沿y方向可划分成2HNn=?个网格，用0,1,2,3,,kN=L来标记。这些网格点可分成三类：4（1）当()01jLn-且0kN，或者当()1LnjM-且nkN时，网格点落在流场W内部，称为内点。这些网格点上的流函数需通过求解方程组来计算；（2）当()1LnjM-?且0kn?时，网格点落在正方形物体内部，网格点上不存在流场，无需计算；（3）其余的网格点落在流场W的边界上，称为边界点。这些网格点上的流函数直接由边界条件给定，也无需计算。7.定解问题的离散化－差分格式Laplace方程22220xyyy抖+=抖的差分近似为1,,1,,1,,122220jkjkjkjkjkjkhhyyyyyy+-+--+-++=边界条件0x¶=¶y的差分近似为,1,0jkjkhyy--=58.内点上流函数的计算－迭代算法在实际的计算中，内点的数量非常多，计算流函数需要求解大型的代数方程组。将Laplace方程的差分近似改写成(),1,1,,1,114jkjkjkjkjkyyyyy+-+-=+++使用高斯-赛德尔迭代算法计算：()(1)()(1)()(1),1,1,,1,114mmmmmjkjkjkjkjkyyyyy++++-+-=+++（0,1,2,3,m=L）给定迭代的初始近似(0),1jky=和迭代的收敛精度15e=-。对0,1,2,3,m=L：在流场内的每一个内点上用迭代公式()(1)()(1)()(1),1,1,,1,114mmmmmjkjkjkjkjkyyyyy++++-+-=+++进行计算。等式右边如果有边界点，可直接将边界条件指定的流函数值代入。在边界CD上的每一个网格点（不包括端点C、D）处，用公式(1)(1),1,mmjkjk++-=yy计算。对所有内点，计算(1)(),,,maxmmjkjkjk+-yy。若(1)(),,,maxmmjkjkjk+-yy，迭代收敛。6否则将m加1，返回步骤39.计算源程序（matlab编制）function[liuxian,vx,vy]=jisuanliuti(x,y)%%%%%%%%%%%%%%%%%%thisisforComputationalFluidDynamics;%%%%%%%%%%%%%%%%使用方式：[liuxian,vx,vy]=jisuanliuti(x,y)；%%%%%%%%%%%%%%%输入x:x轴方向单位长度上划分的网格（x=3）%%%%%%%%%%%%%%%输入y:y方向上单位长度上划分的网格(y=3)%%%%%%%%%%%%%%输出liuxian:输出各个网格点的流函数值%%%%%%%%%%%%%%输出vx：输出x方向的速度u%%%%%%%%%%%%%%输出vy：输出y方向的速度v%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%if(x3)||(y3)error('numberofgridscannotbedividedtoolittle');end%%%%%%%%%%%%%%%%参数定义%%%%%%%%%%%%%%%%%%%%%%%%%%%%7n=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;fori=1:1000forj=2:1:3*nfork=2:1:3*ht=(j2*n+1)||((j=2*n+1)&&(k=h+2));ift==1fai1(j,k)=1/4*(fai1(j+1,k)+fai1(j-1,k)+fai1(j,k+1)+fai1(j,k-1));endend8endc=max(max(abs(fai1-fai)));ifc1e-5fai=fai1;elsebreak;endfai1(end,:)=fai1(end-1,:);endliuxian=fai';%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%画等高线图%%%%%%%%%%%%%%%%%figure;contour(0:1/n:3,0:1/h:3,liuxian,8,'LineWidth',2);title('流线图','fontsize',14);%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%改变流函数矩阵格式%%%%%%%%%%%%fori=1:length(liuxian(:,1))liuxian1(i,:)=liuxian(length(liuxian(:,1))+1-i,:);endliuxian=liuxian1;%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%速度的计算%%%%%%%%%%%%%%%%%%%%%%%%%x方向速度%%%%%%%%%%%%%%%%vx=(liuxian(1:end-1,:)-liuxian(2:end,:))*h;9%%%%%%%%%%%y方向速度%%%%%%%%%%%%%%%vy=-(liuxian(:,2:end)-liuxian(:,1:end-1))*n;end三、实验结果xy流线图00.511.522.5300.511.522.53边界CD上的速度u的值如下表：U1.67651.6561.6151.5581.502Y11.1251.3751.6251.875U1.4561.4231.4011.391.3845Y2.1252.3752.6252.8753四、讨论1.x,y方向划分不同的网格计算精度有所不同，划分网格越多，精度越高。下图是分别按（3,3）和（5,5）划分网格时的流线图。10流线图00.511.522.5300.511.522.53流线图00.511.522.5300.511.522.53流线图00.511.522.5300.511.522.53Figure1网格划分为（5，5）从两幅图处可以清晰地看到，网格数较多的得到的流线更加光滑。2.网格形状的影响分别考虑了（3,5）和（5,3）两种网格划分方式，对比如下可见那个方向划分的网格数较多，那个方向的计算更加精细。流线图00.511.522.5300.511.522.53Figure3网格划分为（3，5）Figure4网格划分为（5，3）Figure2网格划分为（3，3）