分层流管道持液率的计算

1、分层流管道持液率的计算专业：石油与天燃气工程:朱盼2012年1月2如图所示倾角为a的多相流管道截面管道直径为D,管道面积为A,图中界面以下为液体，以上为气体，液面髙度为 液体截面为A】，气体截面积为Ag,分层流的计算模型如下：式中：G 重力加速度；P.气相密度；C液相密度；J罠一一气相与管壁的剪切应力；5液相与管壁的剪切应力；C气液相间剪切应力。分层流模型中各参数有以下关系： 、液位高度同气液分解线与管截面圆心所圆心角“0 = 2 cos-1 1 -2-1 、截面含液率.与气液分解线与管截面圆心所圆心角的关系 、气液相所占管路横截面积A=(1_Hl)A A = HlA 、气相和液相的管壁湿周和

2、气液分界面的宽度0小(8、S=D兀一一，S. = Dsin 11 2丿1匕丿42 、气液相平均流速 、气液各相与管壁的单位长度剪切力及气液相间单位长度剪切力r =f空1r = f门帆岡Wg J Wg 2 wl J wl 2 、气液各相宁系数的求法当 ReW2000 时 /= ；当 血2000 时 / = 348_41og Re气液相间宁系数:当 wsg W wsg,t 时丄=1;股为光滑流转变为波状流的临界表观气速度，它的表达式为叫星=5101325现在假设对气液各相的Re均小于2000,就有人* =/WI=Y , Z =0.0142将 以上参数带入分层流计算公式中，结果为：16 $(伦 乎滋

3、门_壬一匹弘(4Ql 24J2*0142& 2 比口， H，nd2 Re 2 (1-HjnJ2 (1-/；)口， 16 4QL 24 sin 4d sin 沪丽話)+ S-2)gsin0其中： = 2cos-(l-2)若已知口，几，Re, d, Q则公式就是关于0的隐式函数，在0到d用二分法迭代求出厲。编程如下：#include#includemath. h#define A=2acos(l2h/d)define B=16*Ql*Ql/(h*3. 14*d*h*d*3. 14*h)define C=16*Pl*4*d*(3. 14-A/2)*B/(R*2*h*3. 14*d*d)define

4、D=16*P2*16*Q2*Q2/(lh)*3. 14*d*d*3. 14*(1-h)define F=4*d*sin(A/2)define equation C-D*2*d*A/(1H)*3. 14*d*d)+0. 0142/D*(F/(h*3. 14*d*d)+F/(l-h)/(3. 14*d*d) + (pl-p2) *9. 8*sina=0double pl,p2,Ql,Q2,d,h,a /*pl为液体密度，p2为气体密度，QI为液体流 量，Q2为气体流量，a为倾角，h为持液率，d为管道直径*/void find extent (double*hl,doubie*h2)double i

5、=0,j=0;*dl=0；if(equation(pl,p2,QI,Q2,d,h,a)0while (equation(pl,p2,Ql,Q2tdthFa,+i)0&equation(pl,p2,QI,Q2,d,h,a,-j)0)；*d2=equation(pl,p2,QI,Q2,d,h,a)=0?i：j；elsewhile(equation(p11 p2,Q1,Q2,d,h,a f +i)0&equation(pltp2tQl,Q2tht d,a,-j)=0?i:j；double get result (double hlt double h2,double accuracy)double

6、hl=(hl+h2)/2；while (fabs(equation(pl,p2,QI,Q2,d,a.h) accuracy)if(equation(pl,p2tQl,Q2,d,a,hl)*equation(pl,p2,Ql,Q2,d,a,h2)0)hl二h2；elseh2=hl；hl=(hl=h2)/2；)return hl；int main()double hl,h2；scanf (tt%l f%l f%l f%l f%l f%l f &p 1, &p2, &Q1, &Q2, &D, &a):fin extent(&hlt&h2):hl=get result(hl,h2,10e-6):pri

7、ntf (，r%lf，f,hl)；File fdit Run Compile Project Options Debug BEdk/v/dtchEdit Line 1 Col 17 Insert Indent Tab Fill Unindent * C:NONAME.Cttinclude ttinclude ttdefine A=2acos(l-2h/d) Hdefine B=16*Q1*Q1/(h*3.14*d*h*d*3.14*h) ttdefine C=16*Pl*4*d*(3.14-A/2)*B/(R*2*h*3.14*d*d) define D=16*P2*16*Q2*Q2/(1-h

8、)*3.14*d*d*3.14*(1-h) ttdefine F=4*d*sin(A/2)Hdefine equation C-D*2*d*A/(1-H)*3.14*d*d)+0.0142/D*(F/(h*3.14*d*d)+F/(l-h)/(3.14*d*d)*(pl-p2)*9.8*sina=0double pl,p乙Ql;Q2,d,h,avoid find extent(double*hl,doubie*h2) double i=0,j=0;*dl=0;if(equation(pl,p2,Ql.Q2td,h,a)0.汕ile(equat ion(pl.p2,Q2.d,h,d.+i)0溜F

9、 i leEdit Run ComDile Project Options DebugBeal0);*d2=equation(pl,p2.Ql,Q2.d,h,a =0?i:j; elsewhile(equation(pl,p2,QlfQ2,d,h,a,*i)0& equa t i on(pl,p2 ? QI,Q2,h.d,a,j j=0?i:j;) double get result(double hl.double h2,double accuracy) double hl=(hl*h2)/2;wh订e(fal)s(equation(Dl,p2,Ql,Q2,d,a,h)accuracy)if(equation(pl,p2,Ql,Q2,d,a,hl)*equation(pl,p2,Ql,Q2,d,a,h2)0) hl=h2;elseile Project Opt ions Debug Break/watchEdit Line 53 Col 1 Insert Indent Tab Fill Unindent M C:NONAME