中国石油大学(北京)瞬变流编程实践.docx

16 中国石油大学（北京）输油管道瞬变流 实验报告实验课程： 实验名称： 班级： 学号： 姓名： 实验台号： 成绩： 实验日期： 年 月 日 题目：条件：一管道输送柴油。管道由首站和末站组成。管径：4577.1mm站间距为130km。泵机组串联。密度842.1, e=0.01mm, 粘度7.0210-6 m2/s, k=1.39109pa，e=207109pa，kv=1.75, 3.83, 8.38, 18.32, 40.1, 87.62,191.6, 419.0, 916.2, 2003.4,2000000（）1#泵 转速 2970rpm2#泵 转速 2970rpm3#泵 转速 2970rpm流量m3/h扬程m功率kw流量m3/h扬程m功率kw流量m3/h扬程m功率kw496.66148.46241.17411.26280.39387.06350.4306.94385.32575.68146.94284.85634.45255.03458.6602.35295.16511.53855.39142.02325.77853.39232.04539.42855.11273.9644.49边界条件:首站给油泵定值30m; 末站进压0.4mpa; 事故工况，末站关阀; 阀门从全开到全关的时间为：10秒。假设：1、正常情况下首站给油泵定值30m，油品进入各主泵及其附件产生的局部磨阻为30m； 2、泵特性曲线中流量扬程曲线符合二次方程规律； 3、管道处于水平，即纵断面图为直线； 4、各阀之间用固定墩支撑，且管材泊松系数为0.3； 5、阀门开度与时间为一次函数关系，满足阀门打开t时间后，阀门开度为(1-0.01t); 6、下游边界为阀门和固定水位罐，其在固定水位罐管道末端阀门产生一恒定压力0.4mpa。简化流程图如下：利用matlab软件编程如下：ticclc;clear %求解各泵的特性曲线p1=496.66 148.46 575.68 146.94 855.39 142.02 ;p2=411.26 280.39 634.45 255.03 853.39 232.04 ;p3=350.4 306.94 602.35 295.16 855.11 273.9 ; %用二次方程求解各泵的a、b值r1=polyfit(p1(:,1)./3600,p1(:,2),2);r2=polyfit(p2(:,1)./3600,p2(:,2),2);r3=polyfit(p3(:,1)./3600,p3(:,2),2); %常数项g=9.81; %重力加速度rou=842.1; %密度od=0.457; %外径e=0.01*0.001; %绝对粗糙度sigma=7.1*0.001; %壁厚id=od-2*sigma; %内径mu=7.02*10(-6); %粘度mew=1.39*109; %水的弹性系数 modulus of elasticitymei=207*109; %钢管弹性模量dbs=130000; %站间距kesi=e/id; %相对粗糙度fi=1-0.32; %修正系数fsp=0.4*106; %末站进压hfsp=fsp/(g*rou); %末站压头t=100; %计算时间 %比例系数f=0.0246*mu0.25/id4.75; %列冰宗公式，单位长度摩阻系数a1=(r1(1)+r2(1)+r3(1); %二次方系数a2=(r1(2)+r2(2)+r3(2); %一次方系数a3=(r1(3)+r2(3)+r3(3); %常数项pf=8*1.75/(pi2*id4*g); %末站阀摩阻 %二分法求根pl=0; %左边界pr=1; %右边界delta=10(-9); %误差 maxl=1+round(log(pr-pl)-log(delta)/log(2); %最大计算步数for k=1:maxl mid=(pl+pr)/2; %求中间点 ypl=-f*dbs*pl1.75-hfsp-pf+a1*pl2+a2*pl+a3; %左边界值 ypr=-f*dbs*pr1.75-hfsp-pf+a1*pr2+a2*pr+a3; %右边界值 ymid=-f*dbs*mid1.75-hfsp-pf+a1*mid2+a2*mid+a3; %中间点边界值 if ymid=0 pl=mid; pr=mid; %已解出根 elseif ypl*ymid0 pl=mid; %左边界去除，求解区域减半 else ypr*ymid0; pr=mid; %右边界去除，求解区域减半 endif (pr-pl)delta, break; end %误差判断endq=(pl+pr)/2; %方程的根 %判断是否处于水力光滑区re=4*q/(pi*id*mu); %雷洛数lamda=0.3164/re0.25; %达西公式，摩阻系数re1=59.7/kesi(8/7); %雷洛准数一if rere1 disp(流动不处于水力光滑区);end %离散网格vow=(mew/(rou*(1+mew*id*fi/(mei*sigma)0.5; %velocity of wave 水击波传播速度 sd=100; %距离步长n=dbs/sd+1; %离散点数st=sd/vow; %特征线方程，时间步长m=floor(t/st); %时间节点数 %初始条件,即稳态时各点状态for i=1:n qt(i,1)=q; dop(i,1)=sd*(i-1); %各点与起点之间的距离 h(i,1)=-f*qt(i,1)1.75*dop(i,1)+a1*qt(i,1)2+a2*qt(i,1)+a3; %各点水头end %初始化阀kv=1.75, 3.83, 8.38, 18.32, 40.1, 87.62,191.6, 419.0, 916.2, 2003.4,2000000;kl=1./kv;tkv=0:10;t=0:st:m*st;ckv=spline(tkv,kl,t); %离散方程系数b=vow*4/(g*pi*id2);r=2*lamda*sd/(g*pi*id3); %瞬态计算for j=2:m for i=1:n if i=1 %左边界 cm=h(i+1,j-1)-b*qt(i+1,j-1)+r*qt(i+1,j-1)*abs(qt(i+1,j-1); if qt(i,j-1)=0 qt(i,j)=0; else qt(i,j)=1/(2*(a1-pf)*(b-a2)-(b-a2)2+4*(a1-pf)*(cm-a3)0.5); if qt(i,j)1&in %内点 cp=h(i-1,j-1)+b*qt(i-1,j-1)+r*qt(i-1,j-1)*abs(qt(i-1,j-1); cm=h(i+1,j-1)-b*qt(i+1,j-1)+r*qt(i+1,j-1)*abs(qt(i+1,j-1); h(i,j)=(cp+cm)/2; qt(i,j)=(cp-h(i,j)/b; else %右边界点 cp=h(i-1,j-1)+b*qt(i-1,j-1)+r*qt(i-1,j-1)*abs(qt(i-1,j-1); if (j-1)*st10 ckvv=1/ckv(j); qt(i,j)=g*pi*id2/(4*ckvv)*(-b+(b2-8*ckvv*(hfsp-cp)/(g*pi*id2)(0.5); else qt(i,j)=0; end h(i,j)=cp-b*qt(i,j); end endendtoc由程序可以得到各点相应的流量扬程随时间变化的数据，再根据如下画图程序可得个点的流量扬程图：ii=30;tt=(ii-1)*st;mmm=min(min(qt),max(max(qt);nnn=min(min(h),max(max(h);figure(1);ax,h1,h2=plotyy(1:n,qt(:,ii),1:n,h(:,ii); set(get(ax(1),ylabel),string,流量q(m3/s),color,r) %y1set(h1,linestyle,-);set(get(ax(2),ylabel),string,扬程h/m,color,k) %y2title(t=,num2str(tt),s 时管道内各点流量扬程图);legend(流量,扬程)box offset(ax(1),ytick,-0.3:0.1:0.7);set(ax(1),ylim,mmm); %设置y轴间隔set(ax(2),ytick,0:100:1500);set(ax(2),ylim,nnn)set(ax,xlim,1 n) % 设置x轴范围hh=figure(1);saveas(hh,t=,num2str(tt,4),s 时管道内各点流量扬程图,bmp);画出部分图像可得：模拟结果分析：1、稳态流量为0.2305m3/s。2、瞬态水击波传播速度为1093.24m/s。3、模拟情况与真实情况吻合度较高。由于管道的摩擦系数使用的是稳定时的量，产生的左边界波动无法避免。4、在模拟时间总长度为1000s, 对比节点间距的选取对流量最值、压力最值和程序计算时长的影响。下表为节点距离分别为0.1km、0.2 km、0.5 km、1 km所对应的流量最值、压力最值和程序计算时长。节点间距km最小流量m3/s最大流量m3/s最小压力pa最大压力pa计算持续时长s0.1-0.20070.6692400554.49278287.90.5471090.2-0.20050.6705400373.79271100.34.9025570.5-0.20090.6706400296.3927