C題數(shù)學(xué)建模2013.doc

題目C:有桿抽油系統(tǒng)的數(shù)學(xué)建模及診斷摘要:2關(guān)鍵詞：3引言31.問(wèn)題一的分析和求解41.1問(wèn)題一陳述41.2問(wèn)題分析與解決62.問(wèn)題二：泵功圖計(jì)算93問(wèn)題三第一問(wèn)陳述173.1模型一173.2 模型二194.問(wèn)題四第二問(wèn)分析20參考文獻(xiàn)21摘要: 有桿抽油系統(tǒng)是一種廣泛應(yīng)用于開(kāi)采深井原油的重要方法。電機(jī)旋轉(zhuǎn)運(yùn)動(dòng)通過(guò)四連桿機(jī)構(gòu)轉(zhuǎn)變?yōu)槌橛蜅U的垂直運(yùn)動(dòng).假設(shè)驢頭的外輪廓線是部分圓弧,電機(jī)勻速運(yùn)動(dòng),懸點(diǎn)E下只掛光桿（光桿下不接其他桿，不抽油，通常用來(lái)調(diào)試設(shè)備），即在附錄4的圖4中以O(shè)為原點(diǎn)建立直角坐標(biāo)系，（參考1）根據(jù)各點(diǎn)坐標(biāo)和E點(diǎn)與A點(diǎn)、B點(diǎn)的運(yùn)動(dòng)規(guī)律求出懸點(diǎn)E的一個(gè)沖程的位移函數(shù)，速度函數(shù)和加速度函數(shù)。并用MATLAB軟件畫(huà)出位移函數(shù)的圖像，與由附件1給出的數(shù)據(jù)擬合的圖像比較，我們的實(shí)際值與理論值很接近.第二問(wèn)中，Gibbs模型是一種有阻尼的波動(dòng)方程（參考6），根據(jù)變量分離法可以求出只有一級(jí)桿和三級(jí)桿的懸點(diǎn)示功圖轉(zhuǎn)化為泵功圖的詳細(xì)計(jì)算過(guò)程。并利用MATLAB軟件模擬出泵功圖（參考4），泵功圖是一個(gè)波動(dòng)圖，圖1，圖2是在不同位移時(shí)的泵功圖。并利用附件1和附件2給出的一級(jí)桿和三級(jí)桿的油井參數(shù)模擬出了這兩口井的泵功圖。問(wèn)題三中第一問(wèn)，模型一參考3（略）。泵的工作的過(guò)程是由三個(gè)基本環(huán)節(jié)組成的，即活塞在泵內(nèi)讓出容積，原油進(jìn)泵和從泵內(nèi)排出原油。在理想的情況下，活塞上，下一次進(jìn)入和排出液體體積都等于活塞體積V，據(jù)此，我們建立模型二，并解出每日的排油量（參考1）。問(wèn)題四第二問(wèn)，通過(guò)解問(wèn)題二中的波動(dòng)方程得到的解，運(yùn)用數(shù)學(xué)方法來(lái)處理逼近真實(shí)的C.關(guān)鍵詞：懸點(diǎn)，Gibbs模型，懸點(diǎn)示功圖，泵功圖引言有桿抽油系統(tǒng)是開(kāi)采原油廣泛使用的方法，電機(jī)旋轉(zhuǎn)運(yùn)動(dòng)轉(zhuǎn)化為抽油桿上下往返周期運(yùn)動(dòng)，帶動(dòng)設(shè)置在桿下端的泵的兩個(gè)閥相繼開(kāi)閉，從而將地下上千米深處蘊(yùn)藏的石油抽到地面上來(lái)。鋼制抽油桿有很多節(jié)連接而成，具有相同直徑的歸為一級(jí)，級(jí)數(shù)從上到下進(jìn)行編號(hào)，可多達(dá)5級(jí)，從上端點(diǎn)到下端點(diǎn)可能長(zhǎng)達(dá)上千米。描述抽油桿中的任意水平截面處的基本信息的通用方法是示功圖。抽油桿上端點(diǎn)稱為懸點(diǎn)，在一個(gè)沖程期間，儀器以一定的時(shí)間僅間隔測(cè)得懸點(diǎn)處的一系列位移數(shù)據(jù)和荷載數(shù)據(jù)，據(jù)此建立懸點(diǎn)示功圖。泵是由柱塞、游動(dòng)閥、固定閥、部分油管等幾個(gè)部件構(gòu)成的抽象概念，泵中柱塞處的示功圖稱為泵功圖，因?yàn)槭艿街T多因素的影響，在同一時(shí)刻，懸點(diǎn)處受的力與柱塞處受的力是不同的，同樣的，在同一時(shí)刻，懸點(diǎn)處的相對(duì)位移與柱塞處的相對(duì)位移也不一樣概念，因此，懸點(diǎn)示功圖與泵功圖是不同的示功圖包含了很多信息，其中就有有效沖程，泵的有效沖程是指泵中柱塞在一個(gè)運(yùn)動(dòng)周期內(nèi)真正實(shí)現(xiàn)從出油口排油的那段沖程。工程上一般根據(jù)示功圖的形狀與理論示功圖進(jìn)行對(duì)比來(lái)判斷抽油機(jī)的工作狀態(tài)。通過(guò)懸點(diǎn)示功圖可以初步診斷該井的工作狀況，如產(chǎn)量、氣體影響、閥門(mén)漏夜、堵沙等等。要精確診斷油井的工作狀況，最好采用泵功圖。然而，泵在地下深處，使用儀器測(cè)試其示功數(shù)據(jù)實(shí)現(xiàn)困難大、成本高。因此，通過(guò)數(shù)學(xué)建模，把懸點(diǎn)示功圖轉(zhuǎn)化為桿上任意點(diǎn)的示功圖并最終確定泵功圖，以準(zhǔn)確診斷該井的工作狀況，是一個(gè)很有價(jià)值的實(shí)際問(wèn)題。近年來(lái)，功圖量油技術(shù)在我國(guó)油田得到了快速發(fā)展，特別是在西部油田得到了大規(guī)模的推廣應(yīng)用，并取得了較好的應(yīng)用效果。由于示功圖煉油技術(shù)可以實(shí)現(xiàn)24小時(shí)的連續(xù)監(jiān)測(cè)計(jì)產(chǎn)，因此對(duì)于產(chǎn)量波動(dòng)較大的油井，共圖計(jì)產(chǎn)更具有明顯的優(yōu)越性。為此，本文以懸點(diǎn)示功圖為基礎(chǔ)，建立了基于井下多相流動(dòng)的油井產(chǎn)量和有效率的泵功圖模型。1.問(wèn)題一的分析和求解1.1問(wèn)題一陳述 光桿懸點(diǎn)運(yùn)動(dòng)規(guī)律： 電機(jī)旋轉(zhuǎn)運(yùn)動(dòng)通過(guò)四連桿機(jī)構(gòu)變?yōu)槌橛蜅U的垂直運(yùn)動(dòng)，假設(shè)驢頭外輪廓線為部分圓弧、電機(jī)勻速運(yùn)動(dòng)，懸點(diǎn)E下只掛光桿（光桿下不接其它桿，不抽油，通常用來(lái)調(diào)試設(shè)備）。請(qǐng)按附錄4給出四連桿各段尺寸，利用附件1的參數(shù)，求出懸點(diǎn)E點(diǎn)的一個(gè)沖程的運(yùn)動(dòng)規(guī)律：位移函數(shù)、速度函數(shù)、加速度函數(shù)。并與有荷載的附件1的懸點(diǎn)位移數(shù)據(jù)進(jìn)行比較。 BOADEOO 懸點(diǎn)E的運(yùn)動(dòng)過(guò)程：t=0時(shí)刻，曲柄滑塊D位于上頂點(diǎn)（）,AB平行于水平面,E對(duì)應(yīng)坐標(biāo)原點(diǎn)（稱為E的下死點(diǎn)），E的位移為0；D運(yùn)動(dòng)到下頂點(diǎn)（）時(shí)，E的的位移達(dá)到最大（稱為E的上死點(diǎn)）；D接著運(yùn)動(dòng)到上頂點(diǎn)（）時(shí)，E又回到0的位置，完成一個(gè)周期（即一個(gè)沖程）。前臂AO=4315mm,后臂BO=2495mm,連桿BD=3675mm,曲柄半徑OD=R=950mm1.2問(wèn)題分析與解決A，B點(diǎn)做的都是部分圓弧運(yùn)動(dòng)，A點(diǎn)從水平位置向上運(yùn)動(dòng)的軌跡是部分圓弧，B點(diǎn)從水平位置向下運(yùn)動(dòng)的軌跡也是一圓弧，且A，B兩點(diǎn)繞同一點(diǎn)o+運(yùn)動(dòng)，所成圓心角相同，以O(shè)點(diǎn)為原點(diǎn)建立直角坐標(biāo)系，如圖所示： OB點(diǎn)坐標(biāo)為(bcos-1,r+l-bsin),D點(diǎn)坐標(biāo)為(rsin,rcos),r=OD表示曲柄的半徑，l=BD表示連桿的長(zhǎng)度，a=AO表示前臂長(zhǎng)，b=BO表示后臂長(zhǎng)。經(jīng)分析有如下公式(bcos-1-rsin)2+(r+l-bsin-rcos)2=l2， （1.1）化簡(jiǎn)得b2+brsincos+r+lb-brcossin=b2+r2+brsin+rl-r+lrcos，sin+(b2+brsin)2+(r+b-brcos）2=b2+r2+brsin+rl-r+lrcos，其中 tan=b2+brsinr+lb-brcos.+=arcsinb2+r2+brsin+rl-r+lrcos(b2+brsin)2+(r+b-brcos）2，=arcsinb2+r2+brsin+rl-r+lrcos(b2+brsin)2+(r+b-brcos）2-，=arcsinb2+r2+brsin+rl-r+lrcos(b2+brsin)2+(r+b-brcos）2-arctan(b2+brsinr+lb-brcos) （1.2）其中=t;E點(diǎn)的位移：s=a圖0是用MATLAB軟件計(jì)算的E點(diǎn)的實(shí)際位移曲線和附件1的理論位移曲線的比較： 圖0E點(diǎn)的速度函數(shù)： v=dsdt （結(jié)果見(jiàn)附表1的程序結(jié)果x4）E點(diǎn)的加速度函數(shù)：a=dvdt (結(jié)果見(jiàn)附表1的程序結(jié)果x5)補(bǔ)充：E點(diǎn)的位移函數(shù)可以近似為： Y=2.4sin(606.27.6(t-607.64)2.問(wèn)題二：泵功圖計(jì)算 井下泵示功圖轉(zhuǎn)化是功圖計(jì)產(chǎn)的基礎(chǔ)。井下泵示功圖轉(zhuǎn)化技術(shù)是以實(shí)測(cè)的懸點(diǎn)示功圖為基礎(chǔ)，即根據(jù)測(cè)得的懸點(diǎn)處運(yùn)動(dòng)的規(guī)律去計(jì)算每級(jí)抽油桿頂部及抽油泵柱塞上的載荷和位移，并繪出相應(yīng)的井下示功圖，然后根據(jù)地面及井下示功圖分析抽油系統(tǒng)的工況。1966年由Gibbs給出了早期的懸點(diǎn)示功圖轉(zhuǎn)化為井下示功圖的模型?，F(xiàn)在我們根據(jù)Gibbs模型，用帶粘滯阻尼的波動(dòng)方程作為描述油桿動(dòng)態(tài)的基本的微分方程(2.1)來(lái)研究井下泵功圖的轉(zhuǎn)化。 下面是完成整個(gè)沖程時(shí)的模型： (2.1) 這里是粘滯阻尼系數(shù)， 為抽油桿柱斷面不同時(shí)間的位移 應(yīng)力傳播速度 時(shí)間 深度 已知由變量分離的方法可得波動(dòng)方程的特解: (2.2)這里分別是關(guān)于和函數(shù)。對(duì)(2.2)微分，并代入(2.1)得： 由于上式每一邊僅包含獨(dú)立的變量，因此 (2.3) 這里是一個(gè)常數(shù)。 分離上述偏微分方程產(chǎn)生兩個(gè)常微分方程 (2.4) (2.5) 特征方程(2.3)解滿足方程(2.1)得： 為任意常數(shù)。將（2.5）代入（2.3）得： Tt+cTt+k2l22a2TkT=0, k=1,2. (2.6) 解此常微分方程得： TKt=Aker1t+Bker2t, k=1,2 (2.7) 其中：r1=-c+d2, r2=-c-d2, d=c2-4k2d2d2l2 則有： u(x,t)=k=1Xk(x)Tk(t) u(x,t)=k=1(Aker1t+Bker2t)sinklx （2.71） u(x,0)=k=1(Ak+Bk) sinklx=(x) ut(x,0)= k=1(Akr1+Bkr2)sinklx=(x) 在Fourier級(jí)數(shù)討論中我們知道，如果f(x)是一個(gè)確定在區(qū)間（-l,l）上的連續(xù)可微函數(shù)，則可以展開(kāi)為Fourier級(jí)數(shù)。 fx=a02+k=1(akcosklx+bksinklx) (2.8) 其中： ak=1l-11f(x)cosklxdx, k=0,1. (2.9) bk=1l-11f(x)sinkl xdx, k=0,1. (2.10)如果f(x)是奇函數(shù)則由于(2.8)中的被積函數(shù)是奇函數(shù)，所以 ak=0（k=0,1,2.）f(x)可以按正弦函數(shù)展開(kāi)為： f(x)=k=1bksinklx (2.11) bk=1l-11f(x)sinkl xdx=2l01f(x)sinklxdx, k=1,2. (2.12)如果f(x)是偶函數(shù)，則由于(2.10)中的被積函數(shù)是奇函數(shù)，所以： bk=0（k=0,1,2) 可以按余弦函數(shù)展開(kāi)為： f(x)=a02+k=1akcosklx (2.13) ak=1l-11f(x)cosklxdx=2l01f(x)sinklxdx,k=0,1,2. (2.14)所以(x)和x若可以作奇偶延拓到區(qū)間（-l,0）上，只要它們可微就可以按公式(2.11)， (2.12)展成正弦級(jí)數(shù)。因此取： Ak+Bk=2l02(x)sinkldAkr1+Bkr2=2ka01()sinkld解得： Ak=1r2-r1(2r2l0lsinkld-2ka01()sinkld Bk= 1r2-r1(2ka01sinkld-2r1l0lsinkld 由胡克定律得： F(x,t)=E fru(x,t)x, E 表示抽油桿的彈性模量，fr表示抽油桿的橫截面積。參考4的算法，用MATLAB 畫(huà)出該模型的泵功圖。 圖1是利用MATLAB模擬出的不同位移的泵功圖：圖1圖2是不同于圖1位移的泵功圖：圖2注：由于該模型是一動(dòng)態(tài)方程，所以MATLAB擬合出來(lái)的圖像也是變化的，故提交的數(shù)據(jù)與該圖不一致。（程序見(jiàn)附表2）圖3是附件1給出的數(shù)據(jù)模擬出的圖形：圖3圖4是附件2給出的數(shù)據(jù)模擬的圖形：圖4第二種情況，未完成一個(gè)沖程的情況：2u（x,t)t2=a22u(x,t)x2-cux,ttu0,t=0, ux0,t=u0ux,0=u1, utx,0=u2由變量分離的方法可得波動(dòng)方程的特解這里分別是關(guān)于和函數(shù)。對(duì)(2.2)微分，并代入(2.1)得：由于上式每一邊僅包含獨(dú)立的變量，因此這里是一個(gè)常數(shù)。 解下面特征值問(wèn)題：由常微分方程可知，有以下三種情況，（下面給出的，均為任意常數(shù)） 時(shí)， 時(shí)，由 推出： 由推出：設(shè)：可解得： (舍去)解得： 時(shí)，由 推出：由推出：解得：當(dāng)時(shí)， 當(dāng)時(shí)， 其中 ，當(dāng)時(shí)， 3問(wèn)題三第一問(wèn)陳述建立兩個(gè)不同的由泵功圖估計(jì)油井產(chǎn)量的模型，其中至少一個(gè)要利用“有效沖程”；并利用附件1，附加2的數(shù)據(jù)分別估算兩口井一天的產(chǎn)液量。3.1模型一泵示功圖的載荷項(xiàng)包括了液體載荷和泵的摩擦載荷，因此，由泵功圖所計(jì)算得到的功率為泵的輸入功率，即泵的指示功率。在指示功率中減去泵的摩擦功率，即有效功率。游動(dòng)閥隙的水力阻力Fvd可按如下公式計(jì)算：Fvd=pd(Ap-Avd+Adr) (3.1)其中：Ap表示抽油泵柱塞橫截面積 Avd表示游動(dòng)閥座孔過(guò)流面積 Ar表示最下部抽油桿柱的橫截面積 pd表示游動(dòng)閥的水力損失，pa,p由下式計(jì)算： pd=nv112ApAv2VP22g (3.2)其中：vp表示柱塞瞬時(shí)運(yùn)動(dòng)速度 1表示流量系數(shù) nv表示游動(dòng)閥個(gè)數(shù)固定閥隙的水力阻力Fvs計(jì)算如下： Fvs=psAp-Avs+Ar (3.3) 其中： Ap表示抽油泵柱塞橫截面積 Avs表示游動(dòng)閥座孔過(guò)流面積 Ar表示最下部抽油桿柱的橫截面積 ps表示固定閥的水力損失，pa,p由下式計(jì)算： ps=112ApAv2VP22g (3.4)根據(jù)液壓流體力學(xué)中縫隙流動(dòng)理論，所建立的柱塞和泵筒之間的摩擦力的計(jì)算公式： Ff=LpDP2Lp+Vp11-2 (3.5) 其中 ：Lp表示柱塞長(zhǎng)度 ，p表示柱塞壓差 表示動(dòng)力粘度 ， 表示偏心比 表示柱塞和泵筒間的偏心矩 vp表示柱塞運(yùn)動(dòng)速度考慮泵閥水力阻力、柱塞與泵筒之間摩擦力對(duì)控油桿柱軸向力的影響，泵輸入節(jié)點(diǎn)抽油桿柱軸向力由下式計(jì)算：柱塞上沖程：ppump=Pp(t)+Ff+Fvs (3.6)柱塞下沖程：ppump=Pp(t)+Ff-Fvd (3.7)抽油泵輸入節(jié)點(diǎn)平均功率,即有效功率:Npump1T=0Tppump1000Vpumpdt (3.8)其中 Vpump表示柱塞順時(shí)運(yùn)動(dòng)速度將式子（3.8）代入仿真模型： Ne=Q1+RHmsg86400 中，可得油井實(shí)際產(chǎn)量模型： Q=86400Npump1+RHmsg (3.9)3.2 模型二 泵的理論排量： 泵的工作的過(guò)程是由三個(gè)幾基本的環(huán)節(jié)所組成，即活塞在泵內(nèi)讓出容積，原油進(jìn)泵和從泵內(nèi)排出原油。在理想的情況下，活塞上，下一次進(jìn)入和排出液體體積都等于活塞體積V: V=式中 活塞截面積 光桿沖程每分鐘的排量為： 式中：沖次，次分每日排量 式中 泵的理論排量，米日 沖數(shù)，次分引用附表1中的數(shù)據(jù)，可以得到每日排量為： =1440噸4. 問(wèn)題四第二問(wèn)分析 Gibbs模型在數(shù)學(xué)上簡(jiǎn)化為“波動(dòng)方程”：2u2t=a22u2x-cut其中a為已知常數(shù)，c為阻尼系數(shù)。大多數(shù)的阻尼系數(shù)并不能完整地反應(yīng)實(shí)際情況。如果能從方程本身和某些數(shù)據(jù)出發(fā)用數(shù)學(xué)方法估計(jì)參數(shù)c，貢獻(xiàn)是很大的，下面詳細(xì)的計(jì)算給出了c的理論推導(dǎo)過(guò)程： 因?yàn)閹в凶枘嵯禂?shù)的波動(dòng)方程為： 而且 其中 從而可以推出阻尼系數(shù)為：附表1：clc;clear all;omiga=7.6*2*pi/60;syms t ;x=omiga*t;a=4.315;b=2.495;l=3.675;r=0.95;x1=-atan(b2+b*r*sin(x)./(b*l+b*r-b*r*cos(x);x2=x1+asin(r2+b2+r*l+b*r*sin(x)-r*(r+l)*cos(x)./sqrt(b2+b*r*sin(x).2+(r+l)*b-b*r*cos(x).2);x3=a*x2;x4=diff(x3,t)x5=diff(x4,t)附表2：clc;clear all;N=800; dx=0.01;dt=0.0005;c=400*dt*dt/dx/dx;g=5;u(1:100,1)=0;x=linspace(0,1,100);u(2:99,1)=1/(2:99)+(1/(2:99).2)*dt;u(2:99,2)=2*u(2:99,1)+c*(u(3:100,1)-2*u(2:99,1)+u(1:98,1);u(2:99,2)=1