哈工大两相流作业

1、精选优质文档-倾情为你奉上2014 年 秋 季学期研究生课程考核（课程考核报告）考核科目: 多相流 学生所在院（系）: 能源科学与工程学院学生所在学科: 动力工程及工程热物理学 生 姓 名: xxx学 号: 14S0020xx学 生 类 别: 学术型考核结果阅卷人研究生课程：多相流考试和作业1流体携带颗粒的流动过程。假设流体和颗粒具有相同的温度，两者之间无质量交换，颗粒在流体携带下通过一个垂直管道，见图所示。计算条件：计算管长为5.0m，管直径为50mm。颗粒直径取为学号数的最后2位数，mm（例如，即颗粒直径为55mm）。颗粒密度为学号数的最后4位数，kg/m3。气体密度为学号数的最后3位数，

2、kg/m3。进口气体速度为：, 其中，ulz为入口轴向速度分量，ulzo为管中心轴向速度，取为5.0m/s。n取学号数的最后1位数（当0时，取学号的最后第2位数）。ulr为入口径向速度分量。并且假设管内气体速度分布与入口具有相同的速度分布。入口处颗粒的初始位置：n取学号的最后第2位数，若大于50，取为学号的最后第1位数。入口处颗粒轴向和径向速度分量为，其中，n取学号数的最后1位数（由于r=0和1.0为壁面，学号尾数为0和50时，取学号的最后第2位数）。Uszo取为2.0m/s。计算中所需要的其他参数自行确定。要求：1. 给出该颗粒运动速度的变化。 2. 给出该颗粒的运动轨迹（颗粒到达壁面或者出

3、口视为颗粒运动结束）。 3. 提供计算的编程。 4. 提供纸质版。初始数据和条件本人学号14S，n=6,所以流体的轴向速度分布为： ，其中，D=0.05m；径向速度：； 颗粒的初始位置：；入口处轴向速度：，其中；入口径向速度：； 物性参数为：，；气体粘度取常温下空气的粘度：。解题思路和步骤直角坐标系下，颗粒相速度满足如下偏微分方程：dud,xdt=(uc,x-ud,x)rpdud,ydt=(uc,y-ud,y)rpdud,zdt=(uc,z-ud,z)rp+g所以当给定初始速度、位移和合适的时间步长后，可对其后的速度和位移进行求解。1固定网格法将整个计算区域划分成均匀的计算网格，以单一网格作为

4、基本计算区域，确定时间步长，计算颗粒的运动速度，判断颗粒的位置。不断缩小时间步长，直至颗粒落到网格的边界上，进入下一网格计算。具体步骤如下：1. 划分网格，因为本题目径向与轴向尺寸差异较大，且径向与轴向速度也相差较大，为了保证计算精度和计算速度，采用径向宽度和轴向高度不相等的长方形网格，径向宽度dx=10-4m，轴向高度dz=2X10-3m；2. 初选时间步长： t1=minx1ux,z1uzx和z分别为颗粒到网格边界的距离;3. 利用四阶龙哥库塔法求解微分方程，求得颗粒的速度ux,t和uz,t;4. 计算颗粒相的位置：xt=x0+0.5×(ux,0+ux,t) zt=z0+0.5&

5、#215;(uz,0+uz,t)；5. 判断颗粒的位置xt和zt是否落在网格边界上（当颗粒到网格边界的距离小于10-7m时即认为已经到达网格边界上），如果落到网格边界上或已出网格，此步计算结束，进入下一网格进行计算；如果落到该计算网格内部，则从新选择时间步长： t2=minx2(ux+ux，t)/2,z2(uz+uz,t)/2之后重复上述2-5步骤即可。2.移动网格法基本思路是不划定位置确定的网格，颗粒每经过一个时间步长t后，再以颗粒现所在的位置为原点重新建立网格（网格大小始终保持一样），确定下一步的时间步长，进而计算颗粒的速度和位移。具体步骤如下：1. 选定网格大小，和固定网格法一样，径向宽

6、度dx=10-4m，轴向高度dz=2X10-3m；2. 选定时间步长：t=mindxux,dzuz；3. 利用四阶龙哥库塔法求解微分方程，求得颗粒的速度ux,t和uz,t;4. 计算颗粒相的位置：xt=x0+0.5×(ux,0+ux,t) zt=z0+0.5×(uz,0+uz,t)；5. 以颗粒现所在的位置为原点重新建立网格,重复2-5步骤即可。计算结果1. 颗粒相运动速度的变化1. 固定网格法上方曲线为轴向速度，下方为径向速度2.移动网格法上方轴向速度，下方径向速度2.运动轨迹1.固定网格法2. 移动网格法附录-程序代码1.0固定网格法%初始数据%r0=0.05/6;%颗

7、粒初始径向位置,单位mz0=0;%颗粒初始轴向位置dr=0.0001;%网格径向长度，单位mdz=20*dr;%网格轴向长度，单位md=6.6e-5;%颗粒的直径，单位mpc=66;%气体密度pd=2066;%颗粒密度uc=17.9e-6;%气体的动力粘度ucr=0;%气体径向速度，为0trp=d2*pc/(18*uc);%颗粒的松弛时间udz=2*(1-(r0/0.05)6);%颗粒初始轴向速度udr=udz/6;%颗粒初始径向速度rt=r0;%初始化径向位移zt=z0;%初始化轴向位移t=0;%初始化时间x(1)=0;%网格起点y(1)=0;%网格起点%建立网格，储存网格内气体轴向速度%f

8、or i=2:(0.05/dr+1) u(i)=0.5*(5*(1-(x(i-1)/0.05)(1/6)+5*(1-(x(i-1)+dr)/0.05)(1/6); x(i)=x(i-1)+dr;endfor k=2:(5/dz+1) y(k)=y(k-1)+dz;end%循环主体%j=84;l=2;m=1;while and(and(rt>=0,rt<=0.05),and(zt>=0,zt<=5),%判断是否出垂直管道 ucz=u(j);%确定该网格内气体的轴向速度值 while and(y(l)-zt)>=1e-7,x(j)-rt>=1e-7),%满足条件

9、则循环 dt=min(abs(x(j)-rt)/udr),abs(y(l)-zt)/udz);%选取最小时间步 %四阶Rongue-Kutta% kz1=(ucz-udz)/trp-(pd-d)/pd*9.8; kz2=(ucz-(udz+kz1*dt/2)/trp-(pd-d)/pd*9.8; kz3=(ucz-(udz+kz2*dt/2)/trp-(pd-d)/pd*9.8; kz4=(ucz-(udz+kz3*dt)/trp-(pd-d)/pd*9.8; kr1=(ucr-udr)/trp; kr2=(ucr-(udr+kr1*dt/2)/trp; kr3=(ucr-(udr+kr2*d

10、t/2)/trp; kr4=(ucr-(udr+kr3*dt)/trp; udzt=udz+dt/6*(kz1+2*kz2+2*kz3+kz4);%计算dt时刻后的颗粒轴向速度 udrt=udr+dt/6*(kr1+2*kr2+2*kr3+kr4);%计算dt时刻后的颗粒径向速度 zt=zt+0.5*(udz+udzt)*dt;%计算dt时刻后的颗粒轴向位置 rt=rt+0.5*(udr+udrt)*dt;%计算dt时刻后的颗粒径向位置 t=t+dt; %输出数据excel% a2(m,1)=rt; a2(m,2)=zt; a2(m,3)=udrt; a2(m,4)=udzt; a2(m,5)

11、=t; m=m+1; %更新初试速度% udr=0.5*（udrt+udr）; udz=0.5* (udzt+udr); fprintf('udz=%dn',udz); fprintf('udr=%dn',udr); fprintf('zt=%dn',zt); fprintf('rt=%dn',rt); fprintf('t=%dn',t); end %判断下一网格位置% if (x(j)-rt<=1e-7 && y(l)-zt>1e-7)|rt>=x(j) j=j+1; else

12、if (x(j)-rt>=1e-7 && y(l)-zt<1e-7)|zt>=y(l) l=l+1; else j=j+1;l=l+1; end end3. 移动网格法%初始数据%r0=0.05/6;%颗粒初始径向位置,单位mz0=0;%颗粒初始轴向位置dr=0.0001;%网格径向长度，单位mdz=20*dr;%网格轴向长度，单位md=6.6e-5;%颗粒的直径，单位mpc=66;%气体密度pd=2066;%颗粒密度uc=17.9e-6;%气体的动力粘度ucr=0;trp=d2*pc/(18*uc);%颗粒的松弛时间udz=2*(1-(r0/0.05)6);

13、%颗粒初始轴向速度udr=udz/6;%颗粒初始径向速度rt=r0;%初始化径向位移zt=z0;%初始化轴向位移t=0;%初始化时间m=1%程序主体%while and(and(rt>=0,rt<=0.05),and(zt>=0,zt<=5),%判断是否出垂直管道 dt=min(dr/udr,dz/udz);%选取最小时间步 ucz=0.5*(5*(1-(rt/0.05)(1/6)+5*(1-(rt+dr)/0.05)(1/6);%确定该网格内气体的轴向速度值 %四阶Rongue-Kutta% kz1=(ucz-udz)/trp-(pd-d)/pd*9.8; kz2=(

14、ucz-(udz+kz1*dt/2)/trp-(pd-d)/pd*9.8; kz3=(ucz-(udz+kz2*dt/2)/trp-(pd-d)/pd*9.8; kz4=(ucz-(udz+kz3*dt)/trp-(pd-d)/pd*9.8; kr1=(ucr-udr)/trp; kr2=(ucr-(udr+kr1*dt/2)/trp; kr3=(ucr-(udr+kr2*dt/2)/trp; kr4=(ucr-(udr+kr3*dt)/trp; udzt=udz+dt/6*(kz1+2*kz2+2*kz3+kz4);%计算dt时刻后的颗粒轴向速度 udrt=udr+dt/6*(kr1+2*kr2+2*kr3+kr4);%计算dt时刻后的颗粒径向速度 zt=zt+0.5*(udz+udzt)*dt;%计算dt时刻后的颗粒轴向位置 rt=rt+0.5*(udr+udrt)*dt;%计算dt时刻后的颗粒径向位置 t=t+dt; %输出表格% a1(m,1)=rt; a1(m,2)=zt; a1(m,3)=udrt; a1(m,4)=udzt; a1(m,5)=t; m=m+1;