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

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

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

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

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

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

6、度dx=10-4m，軸向高度dz=2X10-3m；2. 選定時(shí)間步長(zhǎng)：t=mindxux,dzuz；3. 利用四階龍哥庫(kù)塔法求解微分方程，求得顆粒的速度ux,t和uz,t;4. 計(jì)算顆粒相的位置：xt=x0+0.5×(ux,0+ux,t) zt=z0+0.5×(uz,0+uz,t)；5. 以顆?，F(xiàn)所在的位置為原點(diǎn)重新建立網(wǎng)格,重復(fù)2-5步驟即可。計(jì)算結(jié)果1. 顆粒相運(yùn)動(dòng)速度的變化1. 固定網(wǎng)格法上方曲線為軸向速度，下方為徑向速度2.移動(dòng)網(wǎng)格法上方軸向速度，下方徑向速度2.運(yùn)動(dòng)軌跡1.固定網(wǎng)格法2. 移動(dòng)網(wǎng)格法附錄-程序代碼1.0固定網(wǎng)格法%初始數(shù)據(jù)%r0=0.05/6;%顆

7、粒初始徑向位置,單位mz0=0;%顆粒初始軸向位置dr=0.0001;%網(wǎng)格徑向長(zhǎng)度，單位mdz=20*dr;%網(wǎng)格軸向長(zhǎng)度，單位md=6.6e-5;%顆粒的直徑，單位mpc=66;%氣體密度pd=2066;%顆粒密度uc=17.9e-6;%氣體的動(dòng)力粘度ucr=0;%氣體徑向速度，為0trp=d2*pc/(18*uc);%顆粒的松弛時(shí)間udz=2*(1-(r0/0.05)6);%顆粒初始軸向速度udr=udz/6;%顆粒初始徑向速度rt=r0;%初始化徑向位移zt=z0;%初始化軸向位移t=0;%初始化時(shí)間x(1)=0;%網(wǎng)格起點(diǎn)y(1)=0;%網(wǎng)格起點(diǎn)%建立網(wǎng)格，儲(chǔ)存網(wǎng)格內(nèi)氣體軸向速度%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%循環(huán)主體%j=84;l=2;m=1;while and(and(rt>=0,rt<=0.05),and(zt>=0,zt<=5),%判斷是否出垂直管道 ucz=u(j);%確定該網(wǎng)格內(nèi)氣體的軸向速度值 while and(y(l)-zt)>=1e-7,x(j)-rt>=1e-7),%滿足條件

9、則循環(huán) dt=min(abs(x(j)-rt)/udr),abs(y(l)-zt)/udz);%選取最小時(shí)間步 %四階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);%計(jì)算dt時(shí)刻后的顆粒軸向速度 udrt=udr+dt/6*(kr1+2*kr2+2*kr3+kr4);%計(jì)算dt時(shí)刻后的顆粒徑向速度 zt=zt+0.5*(udz+udzt)*dt;%計(jì)算dt時(shí)刻后的顆粒軸向位置 rt=rt+0.5*(udr+udrt)*dt;%計(jì)算dt時(shí)刻后的顆粒徑向位置 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)

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 %判斷下一網(wǎng)格位置% 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. 移動(dòng)網(wǎng)格法%初始數(shù)據(jù)%r0=0.05/6;%顆粒初始徑向位置,單位mz0=0;%顆粒初始軸向位置dr=0.0001;%網(wǎng)格徑向長(zhǎng)度，單位mdz=20*dr;%網(wǎng)格軸向長(zhǎng)度，單位md=6.6e-5;%顆粒的直徑，單位mpc=66;%氣體密度pd=2066;%顆粒密度uc=17.9e-6;%氣體的動(dòng)力粘度ucr=0;trp=d2*pc/(18*uc);%顆粒的松弛時(shí)間udz=2*(1-(r0/0.05)6);

13、%顆粒初始軸向速度udr=udz/6;%顆粒初始徑向速度rt=r0;%初始化徑向位移zt=z0;%初始化軸向位移t=0;%初始化時(shí)間m=1%程序主體%while and(and(rt>=0,rt<=0.05),and(zt>=0,zt<=5),%判斷是否出垂直管道 dt=min(dr/udr,dz/udz);%選取最小時(shí)間步 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)/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);%計(jì)算dt時(shí)刻后的顆粒軸向速度 udrt=udr+dt/6*(kr1+2*kr2+2*kr3+kr4);%計(jì)算dt時(shí)刻后的顆粒徑向速度 zt=zt+0.5*(udz+udzt)*dt;%計(jì)算dt時(shí)刻后的顆粒軸向位置 rt=rt+0.5*(udr+udrt)*dt;%計(jì)算dt時(shí)刻后的顆粒徑向位置 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;