哈工大兩相流作業(yè)

1、2014 年 秋 季學期研究生課程考核（課程考核報告）考核科目: 多相流 學生所在院（系）: 能源科學與工程學院學生所在學科: 動力工程及工程熱物理學 生 姓 名: xxx學 號: 14S0020xx學 生 類 別: 學術型考核結(jié)果閱卷人研究生課程：多相流考試和作業(yè)1流體攜帶顆粒的流動過程。假設流體和顆粒具有相同的溫度，兩者之間無質(zhì)量交換，顆粒在流體攜帶下通過一個垂直管道，見圖所示。計算條件：計算管長為5.0m，管直徑為50mm。顆粒直徑取為學號數(shù)的最后2位數(shù)，mm（例如，即顆粒直徑為55mm）。顆粒密度為學號數(shù)的最后4位數(shù)，kg/m3。氣體密度為學號數(shù)的最后3位數(shù)，kg/m3。進口氣體速度為

2、：, 其中，ulz為入口軸向速度分量，ulzo為管中心軸向速度，取為5.0m/s。n取學號數(shù)的最后1位數(shù)（當0時，取學號的最后第2位數(shù)）。ulr為入口徑向速度分量。并且假設管內(nèi)氣體速度分布與入口具有相同的速度分布。入口處顆粒的初始位置：n取學號的最后第2位數(shù)，若大于50，取為學號的最后第1位數(shù)。入口處顆粒軸向和徑向速度分量為，其中，n取學號數(shù)的最后1位數(shù)（由于r=0和1.0為壁面，學號尾數(shù)為0和50時，取學號的最后第2位數(shù)）。Uszo取為2.0m/s。計算中所需要的其他參數(shù)自行確定。要求：1. 給出該顆粒運動速度的變化。 2. 給出該顆粒的運動軌跡（顆粒到達壁面或者出口視為顆粒運動結(jié)束）。 3

3、. 提供計算的編程。 4. 提供紙質(zhì)版。初始數(shù)據(jù)和條件本人學號14S，n=6,所以流體的軸向速度分布為： ，其中，D=0.05m；徑向速度：； 顆粒的初始位置：；入口處軸向速度：，其中；入口徑向速度：； 物性參數(shù)為：，；氣體粘度取常溫下空氣的粘度：。解題思路和步驟直角坐標系下，顆粒相速度滿足如下偏微分方程：dud,xdt=(uc,x-ud,x)rpdud,ydt=(uc,y-ud,y)rpdud,zdt=(uc,z-ud,z)rp+g所以當給定初始速度、位移和合適的時間步長后，可對其后的速度和位移進行求解。1固定網(wǎng)格法將整個計算區(qū)域劃分成均勻的計算網(wǎng)格，以單一網(wǎng)格作為基本計算區(qū)域，確定時間步長

4、，計算顆粒的運動速度，判斷顆粒的位置。不斷縮小時間步長，直至顆粒落到網(wǎng)格的邊界上，進入下一網(wǎng)格計算。具體步驟如下：1. 劃分網(wǎng)格，因為本題目徑向與軸向尺寸差異較大，且徑向與軸向速度也相差較大，為了保證計算精度和計算速度，采用徑向?qū)挾群洼S向高度不相等的長方形網(wǎng)格，徑向?qū)挾萪x=10-4m，軸向高度dz=2X10-3m；2. 初選時間步長： t1=minx1ux,z1uzx和z分別為顆粒到網(wǎng)格邊界的距離;3. 利用四階龍哥庫塔法求解微分方程，求得顆粒的速度ux,t和uz,t;4. 計算顆粒相的位置：xt=x0+0.5(ux,0+ux,t) zt=z0+0.5(uz,0+uz,t)；5. 判斷顆粒的

5、位置xt和zt是否落在網(wǎng)格邊界上（當顆粒到網(wǎng)格邊界的距離小于10-7m時即認為已經(jīng)到達網(wǎng)格邊界上），如果落到網(wǎng)格邊界上或已出網(wǎng)格，此步計算結(jié)束，進入下一網(wǎng)格進行計算；如果落到該計算網(wǎng)格內(nèi)部，則從新選擇時間步長： t2=minx2(ux+ux，t)/2,z2(uz+uz,t)/2之后重復上述2-5步驟即可。2.移動網(wǎng)格法基本思路是不劃定位置確定的網(wǎng)格，顆粒每經(jīng)過一個時間步長t后，再以顆粒現(xiàn)所在的位置為原點重新建立網(wǎng)格（網(wǎng)格大小始終保持一樣），確定下一步的時間步長，進而計算顆粒的速度和位移。具體步驟如下：1. 選定網(wǎng)格大小，和固定網(wǎng)格法一樣，徑向?qū)挾萪x=10-4m，軸向高度dz=2X10-3m；

6、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. 以顆?，F(xiàn)所在的位置為原點重新建立網(wǎng)格,重復2-5步驟即可。計算結(jié)果1. 顆粒相運動速度的變化1. 固定網(wǎng)格法上方曲線為軸向速度，下方為徑向速度2.移動網(wǎng)格法上方軸向速度，下方徑向速度2.運動軌跡1.固定網(wǎng)格法2. 移動網(wǎng)格法附錄-程序代碼1.0固定網(wǎng)格法%初始數(shù)據(jù)%r0=0.05/6;%顆粒初始徑向位置,單位mz0=0;%顆粒初始軸向位置dr=0.0001;%網(wǎng)

7、格徑向長度，單位mdz=20*dr;%網(wǎng)格軸向長度，單位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;%網(wǎng)格起點y(1)=0;%網(wǎng)格起點%建立網(wǎng)格，儲存網(wǎng)格內(nèi)氣體軸向速度%for i=2:(0.05/dr+1) u(i)=0.5*(5*(1-(x(

8、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%循環(huán)主體%j=84;l=2;m=1;while and(and(rt=0,rt=0,zt=1e-7,x(j)-rt=1e-7),%滿足條件則循環(huán) 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)/tr

9、p-(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時刻

10、后的顆粒徑向速度 zt=zt+0.5*(udz+udzt)*dt;%計算dt時刻后的顆粒軸向位置 rt=rt+0.5*(udr+udrt)*dt;%計算dt時刻后的顆粒徑向位置 t=t+dt; %輸出數(shù)據(jù)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=%dn,udz); fprintf(udr=%dn,udr); fprintf(zt=%dn,zt); fprintf

11、(rt=%dn,rt); fprintf(t=%dn,t); end %判斷下一網(wǎng)格位置% if (x(j)-rt1e-7)|rt=x(j) j=j+1; elseif (x(j)-rt=1e-7 & y(l)-zt=y(l) l=l+1; else j=j+1;l=l+1; end end3. 移動網(wǎng)格法%初始數(shù)據(jù)%r0=0.05/6;%顆粒初始徑向位置,單位mz0=0;%顆粒初始軸向位置dr=0.0001;%網(wǎng)格徑向長度，單位mdz=20*dr;%網(wǎng)格軸向長度，單位md=6.6e-5;%顆粒的直徑，單位mpc=66;%氣體密度pd=2066;%顆粒密度uc=17.9e-6;%氣體的動力粘度

12、ucr=0;trp=d2*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,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);%確定該網(wǎng)格內(nèi)氣體的軸向速度值 %四階Rongue-Kutta% kz1=(ucz-udz)/

13、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; %輸出表格% a1(m,1)=rt; a1(m,2)=z