哈工大两相流作业

2014年秋季学期研究生课程考核（课程考核报告）考核科目:多相流学生所在院（系）:能源科学与工程学院学生所在学科:动力工程及工程热物理学生姓名:xxx学号:14S0020xx学生类别:学术型考核结果阅卷人研究生课程：多相流考试和作业1．流体携带颗粒的流动过程。假设流体和颗粒具有相同的温度，两者之间无质量交换，颗粒在流体携带下通过一个垂直管道，见图所示。计算条件：计算管长为5.0m，管直径为50mm。颗粒直径取为学号数的最后2位数，m（例如20141055，即颗粒直径为55m）。颗粒密度为学号数的最后4位数，kg/m3。气体密度为学号数的最后3位数，kg/m3。进口气体速度为：l,z1/nlzu[1(r/D)]u,o,l,ru0其中，ulz为入口轴向速度分量，ulzo为管中心轴向速度，取为5.0m/s。n取学号数的最后1位数（当0时，取学号的最后第2位数）。ulr为入口径向速度分量。并且假设管内气体速度分布与入口具有相同的速度分布。入口处颗粒的初始位置：r(n/D)n取学号的最后第2位数，若大于50，取为学号的最后第1位数。入口处颗粒轴向和径向速度分量为nszszu[1(r/D)]uo，srszuu/n其中，n取学号数的最后1位数（由于r=0和1.0为壁面，学号尾数为0和50时，取学号的最后第2位数）。Uszo取为2.0m/s。计算中所需要的其他参数自行确定。要求：1.给出该颗粒运动速度的变化。2.给出该颗粒的运动轨迹（颗粒到达壁面或者出口视为颗粒运动结束）。3.提供计算的编程。4.提供纸质版。初始数据和条件本人学号14S002066，n=6,所以流体的轴向速度分布为：l,z1/6lzu[1(r/D)]u,o，其中l,zu5.0m/s，D=0.05m；径向速度：l,ru0；颗粒的初始位置：rD/6；入口处轴向速度：6szszu[1(r/D)]uo，其中sz,0u2.0m/s；入口径向速度：srszuu/6；物性参数为：3s2066kg/m，3c066kg/m，sd66m；气体粘度取常温下空气的粘度：51.789410kg/(ms)。解题思路和步骤直角坐标系下，颗粒相速度满足如下偏微分方程：𝑑𝑢𝑑,𝑥𝑑𝑡=(𝑢𝑐,𝑥−𝑢𝑑,𝑥)τ𝑟𝑝𝑑𝑢𝑑,𝑦𝑑𝑡=(𝑢𝑐,𝑦−𝑢𝑑,𝑦)τ𝑟𝑝𝑑𝑢𝑑,𝑧𝑑𝑡=(𝑢𝑐,𝑧−𝑢𝑑,𝑧)τ𝑟𝑝+𝑔所以当给定初始速度、位移和合适的时间步长后，可对其后的速度和位移进行求解。1固定网格法将整个计算区域划分成均匀的计算网格，以单一网格作为基本计算区域，确定时间步长，计算颗粒的运动速度，判断颗粒的位置。不断缩小时间步长，直至颗粒落到网格的边界上，进入下一网格计算。具体步骤如下：1.划分网格，因为本题目径向与轴向尺寸差异较大，且径向与轴向速度也相差较大，为了保证计算精度和计算速度，采用径向宽度和轴向高度不相等的长方形网格，径向宽度dx=10-4m，轴向高度dz=2X10-3m；2.初选时间步长：∆𝑡(1)=min{∆𝑥(1)𝑢𝑥,∆𝑧(1)𝑢𝑧}∆x和∆z分别为颗粒到网格边界的距离;3.利用四阶龙哥库塔法求解微分方程，求得颗粒的速度𝑢𝑥,t和𝑢𝑧,t;4.计算颗粒相的位置：𝑥𝑡=𝑥0+0.5×(𝑢𝑥,0+𝑢𝑥,𝑡)𝑧𝑡=𝑧0+0.5×(𝑢𝑧,0+𝑢𝑧,𝑡)；5.判断颗粒的位置xt和zt是否落在网格边界上（当颗粒到网格边界的距离小于10-7m时即认为已经到达网格边界上），如果落到网格边界上或已出网格，此步计算结束，进入下一网格进行计算；如果落到该计算网格内部，则从新选择时间步长：∆𝑡(2)=min{∆𝑥(2)(𝑢𝑥+𝑢x，t)/2,∆𝑧(2)(𝑢𝑧+𝑢𝑧,𝑡)/2}之后重复上述2-5步骤即可。2.移动网格法基本思路是不划定位置确定的网格，颗粒每经过一个时间步长∆t后，再以颗粒现所在的位置为原点重新建立网格（网格大小始终保持一样），确定下一步的时间步长，进而计算颗粒的速度和位移。具体步骤如下：1.选定网格大小，和固定网格法一样，径向宽度dx=10-4m，轴向高度dz=2X10-3m；2.选定时间步长：∆t=min{𝑑𝑥𝑢𝑥,𝑑𝑧𝑢𝑧}；3.利用四阶龙哥库塔法求解微分方程，求得颗粒的速度𝑢𝑥,t和𝑢𝑧,t;4.计算颗粒相的位置：𝑥𝑡=𝑥0+0.5×(𝑢𝑥,0+𝑢𝑥,𝑡)𝑧𝑡=𝑧0+0.5×(𝑢𝑧,0+𝑢𝑧,𝑡)；5.以颗粒现所在的位置为原点重新建立网格,重复2-5步骤即可。计算结果1.颗粒相运动速度的变化1.固定网格法012340.00.10.20.30.40.50.60.70.80.91.01.11.21.31.41.51.61.71.81.9速度（m/s）时间（s）轴向速度径向速度-0.02-0.010.000.010.020.00.51.01.5上方曲线为轴向速度，下方为径向速度2.移动网格法-0.50.00.51.01.52.02.53.03.54.00.00.10.20.30.40.50.60.70.80.91.01.11.21.31.41.51.61.71.8速度（m/s）时间轴向速度径向速度-0.020.000.020.00.51.01.5上方轴向速度，下方径向速度2.运动轨迹1.固定网格法0.0000.0050.0100.0150.0200.0250.0300.0350.0400.0450.050012345轴向位移（m）径向位移（m）0.00800.00850.0090-0.10.00.12.移动网格法0.000.010.020.030.040.05012345轴向位移（m）径向位移（m）0.00800.00850.0090-0.10.00.1附录-程序代码1.0固定网格法%%%%%%%%初始数据%%%%%%%%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;%气体径向速度，为0trp=d^2*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;%网格起点%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%建立网格，储存网格内气体轴向速度%%%%%fori=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;endfork=2:(5/dz+1)y(k)=y(k-1)+dz;end%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%循环主体%%%%%j=84;l=2;m=1;whileand(and(rt=0,rt=0.05),and(zt=0,zt=5)),%判断是否出垂直管道ucz=u(j);%确定该网格内气体的轴向速度值whileand((y(l)-zt)=1e-7,x(j)-rt=1e-7),%满足条件则循环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*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;%输出数据excel%a2(m,1)=rt;a2(m,2)=zt;a2(m,3)=udrt;a2(m,4)=udzt;a2(m,5)=t;m=m+1;%更新初试速度%udr=0.5*（udrt+udr）;udz=0.5*(udzt+udr);fprintf('udz=%d\n',udz);fprintf('udr=%d\n',udr);fprintf('zt=%d\n',zt);fprintf('rt=%d\n',rt);fprintf('t=%d\n',t);end%判断下一网格位置%if(x(j)-rt=1e-7&&y(l)-zt1e-7)|rt=x(j)j=j+1;elseif(x(j)-rt=1e-7&&y(l)-zt1e-7)|zt=y(l)l=l+1;elsej=j+1;l=l+1;endend3.移动网格法%%%%%%%%初始数据%%%%%%%%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=d^2*pc/(18*uc);%颗粒的松弛时间udz=2*(1-(r0/0.05)^6);%颗粒初始轴向速度udr=udz/6;%颗粒初始径向速度rt=r0;%初始化径向位移zt=z0;%初始化轴向位移t=0;%初始化时间m=1%%%%%%%程序主体%%%%%%%whileand(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=(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+ud