1、2020/8/25,中國(guó)石油大學(xué),1,管道瞬變流動(dòng)分析,C語(yǔ)言編程大作業(yè),2020/8/25,中國(guó)石油大學(xué)儲(chǔ)建研07-2班,2,問(wèn)題的提出,一條2195mm的水庫(kù)引水管道，L= 5000m， 水庫(kù)液位表面相對(duì)于末端閥門的高度為100m，穩(wěn)定運(yùn)行時(shí)，管內(nèi)流速為0.5m/s，水力坡降i=0.002，a=1100m/s。如果末端閥門瞬時(shí)關(guān)閉，試求管內(nèi)流動(dòng)瞬變狀態(tài)，即流量和流速,2020/8/25,中國(guó)石油大學(xué)儲(chǔ)建研07-2班,3,背景知識(shí),當(dāng)閥門瞬間關(guān)閉時(shí)，由于水流的慣性作用，在閥門的上游部分產(chǎn)生壓力升高，下游部分產(chǎn)生壓力降低。同時(shí)這種壓力波動(dòng)會(huì)沿著管道向上下游傳播，引起管道內(nèi)的瞬變流動(dòng)，這種水流的

2、不穩(wěn)定現(xiàn)象也稱為水擊現(xiàn)象。 求解水擊中壓力流量變化的方法有多種，其中特征線法比較常用。,2020/8/25,中國(guó)石油大學(xué)儲(chǔ)建研07-2班,4,特征線法Characteristic,假設(shè)管道中某一點(diǎn)產(chǎn)生一個(gè)擾動(dòng)，擾動(dòng)會(huì)以波的形式，以dx/dt=a的速度向下游方向傳播，同時(shí)又會(huì)以-a的速度向上游傳播。,0,x,t,P,A,B,C+,C-,如圖所示，C+、C-兩條特征線可以理解為向上游、下游傳播的運(yùn)動(dòng)曲線，斜率為1/a。,2020/8/25,中國(guó)石油大學(xué)儲(chǔ)建研07-2班,5,簡(jiǎn)單管路特征結(jié)點(diǎn)圖,求解特征線的壓頭，有下列公式： C+:Hpi=CP-BQpi (1) C-: Hpi=CM+BQpi CP

3、 CM為常數(shù) (2) CP=Hi-1+BQi-1-RQi-1|Qi-1| CM=Hi-1-BQi-1+RQi+1|Qi+1| (1)和(2)可得 Hpi=(CP+CM)/2 代入任何一式即可得到Qpi,x,1,2,3,N,N+1,t=0,t,2t,3t,P1,P2,P3,PN,PN+1,C+,C-,2020/8/25,中國(guó)石油大學(xué)儲(chǔ)建研07-2班,6,邊界條件,x,1,2,3,N,N+1,t=0,t,2t,3t,P1,P2,P3,PN,PN+1,C+,C-,圖中邊界處的結(jié)點(diǎn)的參數(shù)，無(wú)法根據(jù)上一層結(jié)點(diǎn)的條件求得 可以利用兩個(gè)特征方程，上游視為恒液位，QP1=(HP1-CM)/B 下游端點(diǎn)流量為0

4、，即QPN+1=0.,2020/8/25,中國(guó)石油大學(xué)儲(chǔ)建研07-2班,7,# include # include # define D 0.209 # define g 9.81 # define pai 3.14 # define n 10 main() int H0=100; int m,k,L=5000,a=1100,IPR=8; float dx,dt,Ad,B,R,f,T,Cp,Cm,Tmax=200; float Q0,v=0.5,i=0.002; float Hn+2,Qn+2,Hpn+2,Qpn+2; FILE *fp; fp=fopen(運(yùn)行結(jié)果顯示（03122608）.d

5、at,w); if(fp=fopen(運(yùn)行結(jié)果顯示.txt,w)=NULL) printf(Errorn); exit(0); dx=L/n; dt=dx/a; Ad=pai*D*D/4; Q0=v*Ad; B=a/g/Ad; f=2*g*i*D/v/v; R=f*dx/2/g/D/Ad/Ad;,2020/8/25,中國(guó)石油大學(xué)儲(chǔ)建研07-2班,8,for(m=1;mn+2;m+) Hm=H0-(m-1)*i*dx; Qm=Q0; fprintf(fp,n穩(wěn)態(tài)變量H n); for(m=1;mn+2;m+) fprintf(fp,%f ,Hm); fprintf(fp,n穩(wěn)態(tài)變量Q n); f

6、or(m=1;mn+2;m+) fprintf(fp,%f ,Qm);,求解穩(wěn)態(tài)時(shí)的參數(shù):,2020/8/25,中國(guó)石油大學(xué)儲(chǔ)建研07-2班,9,T=0.0; k=0; while(TTmax) fprintf(fp,nn輸出到此結(jié)束); else for(m=2;mn+1;m+) /求內(nèi)部節(jié)點(diǎn)的H,Q值 Cp=Hm-1+B*Qm-1-R*Qm-1*fabs(Qm-1); Cm=Hm+1-B*Qm+1+R*Qm+1*fabs(Qm+1); Hpm=(Cp+Cm)/2.0; Qpm=(Hpm-Cm)/B; ,2020/8/25,中國(guó)石油大學(xué)儲(chǔ)建研07-2班,10,Cm=H2-B*Q2+R*Q2*

7、fabs(Q2); Hp1=H1; /求前端點(diǎn)的H,Q值 Qp1=(Hp1-Cm)/B; Cp=Hn+B*Qn-R*Qn*fabs(Qn); Hpn+1=Cp; /求后端點(diǎn)的H,Q值 Qpn+1=0.0; for(m=1;mn+2;m+) Hm=Hpm; Qm=Qpm; ,2020/8/25,中國(guó)石油大學(xué)儲(chǔ)建研07-2班,11,if(k/IPR*IPR=k) /控制輸出次數(shù) fprintf(fp,nnn時(shí)間T：%fn,T); fprintf(fp,瞬態(tài)變量Hpn); for(m=1;mn+2;m+) fprintf(fp,%f ,Hm); fprintf(fp,n瞬態(tài)變量Qpn); for(m

8、=1;mn+2;m+) fprintf(fp,%f ,Qm); fclose(fp); ,2020/8/25,中國(guó)石油大學(xué)儲(chǔ)建研07-2班,12,穩(wěn)態(tài)變量H 100.000000 99.000000 98.000000 97.000000 96.000000 95.000000 94.000000 93.000000 92.000000 91.000000 90.000000 穩(wěn)態(tài)變量Q 0.017145 0.017145 0.017145 0.017145 0.017145 0.017145 0.017145 0.017145 0.017145 0.017145 0.017145 時(shí)間T：3.636364 瞬態(tài)變量Hp 100.000000 99.000000 98.000000 149.568848 149.067429 149.564972 149.064774 149.563599 149.064056 149.563507 149.064133 瞬態(tài)變量Qp 0.017145 0.017145 0.017145 0.001069 0