r0=0.05/6;%颗粒初始径向位置,单位m
z0=0;%颗粒初始轴向位置
dr=0.0001;%网格径向长度，单位m
dz=20*dr;%网格轴向长度，单位m
d=6.6e-5;%颗粒的直径，单位m
pc=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
%%%%%%%程序主体%%%%%%%
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;
pc=66;%气体密度
pd=2066;%颗粒密度
uc=17.9e-6;%气体的动力粘度
ucr=0;%气体径向速度，为0
trp=d^2*pc/(18*uc);%颗粒的松弛时间
udz=2*(1-(r0/0.05)^6);%颗粒初始轴向速度
udr=udz/6;%颗粒初始径向速度
rt=r0;%初始化径向位移
4.计算颗粒相的位置：
；
5.判断颗粒的位置xt和zt是否落在网格边界上（当颗粒到网格边界的距离小于10-7m时即认为已经到达网格边界上），如果落到网格边界上或已出网格，此步计算结束，进入下一网格进行计算；如果落到该计算网格内部，则从新选择时间步长：
之后重复上述2-5步骤即可。
2.
基本思路是不划定位置确定的网格，颗粒每经过一个时间步长 后，再以颗粒现所在的位置为原点重新建立网格（网格大小始终保持一样），确定下一步的时间步长，进而计算颗粒的速度和位移。具体步骤如下：
2.移动网格法
上方轴向速度，下方径向速度
2.
1.固定网格法
2.移动网格法
附录
1.0固定网格法
%%%%%%%%初始数据%%%%%%%%
r0=0.05/6;%颗粒初始径向位置,单位m
z0=0;%颗粒初始轴向位置
dr=0.0001;%网格径向长度，单位m
dz=20*dr;%网格轴向长度，单位m
d=6.6e-5;%颗粒的直径，单位m
要求：1.给出该颗粒运动速度的变化。
2.给出该颗粒的运动轨迹（颗粒到达壁面或者出口视为颗粒运动结束）。
3.提供计算的编程。
4.提供纸质版。
初始数据和条件
本人学号14S******，n=6,所以流体的轴向速度分布为：
，其中 ，D=0.05m；
径向速度：
；
颗粒的初始位置：
；
入口处轴向速度：
，其中 ；
入口径向速度：
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;
1.划分网格，因为本题目径向与轴向尺寸差异较大，且径向与轴向速度也相差较大，为了保证计算精度和计算速度，采用径向宽度和轴向高度不相等的长方形网格，径向宽度dx=10-4m，轴向高度dz=2X10-3m；
2.初选时间步长：
和 分别为颗粒到网格边界的距离;
3.利用四阶龙哥库塔法求解微分方程，求得颗粒的速度 和 ;
end
if udz<0
dz=-dz;
end
end
end
%判断下一网格位置%
if (x(j)-rt<=1e-7 && y(l)-zt>1e-7)|rt>=x(j)
j=j+1;
elseif (x(j)-rt>=1e-7 && y(l)-zt<1e-7)|zt>=y(l)
l=l+1;
else
j=j+1;l=l+1;
end
end
3.移动网格法
%%%%%%%%初始数据%%%%%%%%
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;
x(i)=x(i-1)+dr;
end
for 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)),%判断是否出垂直管道
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时刻后的颗粒径向速度
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时刻后的颗粒径向速度
；
物性参数为：
，
，
；
气体粘度取常温下空气的粘度：
。
解题思路和步骤
直角坐标系下，颗粒相速度满足如下偏微分方程：
所以当给定初始速度、位移和合适的时间步长后，可对其后的速度和位移进行求解。
1
将整个计算区域划分成均匀的计算网格，以单一网格作为基本计算区域，确定时间步长，计算颗粒的运动速度，判断颗粒的位置。不断缩小时间步长，直至颗粒落到网格的边界上，进入下一网格计算。具体步骤如下：
计算条件：计算管长为5.0m，管直径为50mm。
颗粒直径取为学号数的最后2位数，m（例如********，即颗粒直径为55m）。颗粒密度为学号数的最后4位数，kg/m3。气体密度为学号数的最后3位数，kg/m3。
进口气体速度为：
,
其中，ulz为入口轴向速度分量，ulzo为管中心轴向速度，取为5.0m/s。n取学号数的最后1位数（当0时，取学号的最后第2位数）。ulr为入口径向速度分量。
ucz=u(j);%确定该网格内气体的轴向速度值
while and((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;
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;
2014年秋季学期研究生课程考核
（课程考核报告）
考核科目
:多相流
学生所在院（系）
:能源科学与工程学院
学生所在学科
:动力工程及工程热物理
学生姓名
:xxx
学号
: 14S0020xx
学生类别
:学术型
考核结果
阅卷人
研究生课程：多相流
考试和作业
1．流体携带颗粒的流动过程。假设流体和颗粒具有相同的温度，两者之间无质量交换，颗粒在流体携带下通过一个垂直管道，见图所示。
1.选定网格大小，和固定网格法一样，径向宽度dx=10-4m，轴向高度dz=2X10-3m；
2.选定时间步长：
；
3.利用四阶龙哥库塔法求解微分方程，求得颗粒的速度 和 ;
4.计算颗粒相的位置：
；
5.以颗粒现所在的位置为原点重新建立网格,重复2-5步骤即可。
计算结果
1.
1.固定网格法
上方曲线为轴向速度，下方为径向速度
kz2=(ucz-(udz+kz1*dt/2))/trp-(pd-d)/pd*9.8;