基于蒙特卡羅算法的管內(nèi)智能封堵器減振結(jié)構(gòu)優(yōu)化設(shè)計

苗興園，趙 弘

(中國石油大學(xué)(北京) 機械與儲運工程學(xué)院，北京 102249)

管道運輸方式在石油、天然氣運輸中起著重要的作用。隨著使用時間的增長，管道性能逐漸下降，會出現(xiàn)腐蝕、泄漏等現(xiàn)象，對環(huán)境造成嚴重污染，甚至引發(fā)重大安全事故[1-3]。近年來，管內(nèi)智能封堵技術(shù)發(fā)展較快，許多企業(yè)和高校都致力于這方面的研究。張仕民等[4-6]完成了管內(nèi)智能封堵器的整體結(jié)構(gòu)設(shè)計，并對其連接機構(gòu)、密封特性和錨定性能進行了改進，使封堵器的總體性能更加完善；中油管道檢測技術(shù)有限責(zé)任公司的李振北等人研究了橡膠筒和卡瓦齒形角對封堵和錨定特性的影響[7-8]，同時設(shè)計出了單節(jié)雙封型封堵器[9]，提高了封堵器的密封特性；趙弘等[10-11]針對封堵器周圍流場進行了模擬分析，并在此基礎(chǔ)上對封堵器端面的形狀和結(jié)構(gòu)進行了優(yōu)化，同時研究了封堵器不同長度、直徑以及長徑比對管內(nèi)流體的影響[12-13]，根據(jù)模擬結(jié)果做出了相應(yīng)的改進。

在管內(nèi)智能封堵器工作過程中，周圍流場會產(chǎn)生較大的壓力脈動，尾部流體區(qū)域出現(xiàn)漩渦，造成大量的能量損失，引起封堵器的振動，甚至導(dǎo)致封堵失敗。通過對擾流板結(jié)構(gòu)的分析，發(fā)現(xiàn)擾流板的角度和面積對管內(nèi)流場存在較大影響。因此，采用響應(yīng)面法設(shè)計試驗方案，建立擾流板角度、面積與封堵前后壓力差的數(shù)學(xué)模型。利用強化學(xué)習(xí)中的蒙特卡羅算法，用封堵前后的壓力差對應(yīng)封堵器流場環(huán)境，將得到的數(shù)學(xué)模型作為流場環(huán)境與擾流板角度、面積之間的映射關(guān)系。根據(jù)該映射關(guān)系，在設(shè)定的參數(shù)范圍內(nèi)進行迭代搜索，獲得壓力差最小時對應(yīng)的擾流板角度和面積，并根據(jù)得到的參數(shù)建立優(yōu)化模型，進行數(shù)值模擬。

1 管內(nèi)智能封堵器減振結(jié)構(gòu)模型的建立

1.1 減振結(jié)構(gòu)及工作原理

管內(nèi)智能封堵器主要包括滑塊、擠壓碗、承壓頭、執(zhí)行器盤、密封圈和推筒等，如圖1所示。封堵開始時，封堵器接收到電磁脈沖信號，在管內(nèi)介質(zhì)的作用下沿管道向前行進，到達指定位置后由動力系統(tǒng)驅(qū)動推筒運動，使滑塊沿擠壓碗表面滑動，其表面螺紋嵌入管道內(nèi)壁，使整個機構(gòu)固定[14-15]。密封圈在受到擠壓力時徑向膨脹，與管道形成過盈配合，實現(xiàn)管道的密封。針對封堵器的尾部結(jié)構(gòu)進行了改進，在原有結(jié)構(gòu)的基礎(chǔ)上設(shè)計了3塊可折疊控制擾流板，α角為擾流板與底盤的夾角，工作時由動力系統(tǒng)帶動彈簧壓縮，支撐桿推動擾流板翻轉(zhuǎn)控制α角的大小。

1.2 控制方程的建立

為了方便分析計算，假設(shè)管內(nèi)介質(zhì)為單相不可壓縮流體，且工作過程中溫度不變。對于絕大多數(shù)流動情況，都需要質(zhì)量守恒方程和動量守恒方程來進行計算[16-17]。

質(zhì)量守恒方程：

(1)

式中：ρ為液體介質(zhì)的密度，kg/m3；ux、uy、uz為微元體在x、y、z方向上的速度分量，m/s。

動量守恒方程：

(2)

(3)

(3)

式中：X、Y、Z為作用在微元體的質(zhì)量力在x、y、z三個方向上的分量，N；v為運動黏性系數(shù)，m2/s。

1—壓縮彈簧；2—底盤；3—推筒；4—執(zhí)行器盤；5—滑塊；6—擠壓碗；7—密封圈；8—承壓頭；9—擾流板。

管內(nèi)流場雷諾數(shù)較高，將液體視為不可壓縮流體，采用壓力-速度耦合方法。為了使流場情況更加貼近實際，選擇標(biāo)準(zhǔn)k-ε湍流模型，利用Simple算法進行計算。該湍流方程包括湍流動能k方程和湍流耗散率ε方程：

(5)

(6)

式中：μ為流體動力粘性系數(shù)，Pa·s；ε為耗散率；κ為湍動能；μt為湍流粘度系數(shù)；其中σk、σε、c1、c2為常數(shù)，分別為1.0、1.3、1.44、1.92。

1.3 封堵器周圍流場模型

管內(nèi)智能封堵器周圍流場模型如圖2所示。

圖2 管內(nèi)智能封堵器周圍流場模型

為了方便建模，管道模型選擇直管形狀，長度為10 m，封堵器中心軸線與管道中心軸線重合，利用ANSYS中的SpaceClaim模塊建立三維簡化模型，參數(shù)如表1所示。數(shù)值模擬時管內(nèi)介質(zhì)為液態(tài)汽油，且為單相流體，其密度為830 kg/m3，黏度系數(shù)為0.332 mPa·s。入口邊界條件為速度入口，流速為3 m/s，水力直徑為0.6 m。出口邊界條件為壓力出口邊界，管壁和封堵器壁面為選擇無滑移壁面邊界條件。管內(nèi)操作壓力設(shè)置為5 MPa，考慮重力影響，方向設(shè)為z軸負向。

表1 管內(nèi)智能封堵器周圍流場模型參數(shù) mm

2 封堵器周圍流場數(shù)值模擬

本次數(shù)值模擬選擇ANSYS中的Meshing模塊進行網(wǎng)格劃分，此模塊操作簡便，精度較高，適用于復(fù)雜幾何模型。網(wǎng)格形狀為四面體結(jié)構(gòu)，將封堵器設(shè)置為固體域，管內(nèi)其它部分設(shè)置為流體域，網(wǎng)格類型為CFD網(wǎng)格，設(shè)定結(jié)點數(shù)量后生成體網(wǎng)格，如圖3所示。

圖3 管內(nèi)智能封堵器網(wǎng)格模型

為了更加清晰地觀察管內(nèi)流場的分布，選擇管道中心截面(z=0的平面)為觀測面。當(dāng)封堵完全時管道將被封死，無法模擬管內(nèi)流體的流動情況，因此選擇封堵狀態(tài)為99 %時進行研究。封堵過程中會發(fā)生壓力脈動，在封堵器尾部出現(xiàn)回流和漩渦現(xiàn)象，使封堵器能量損失嚴重，在封堵時發(fā)生振動現(xiàn)象，無法保證封堵的穩(wěn)定性。為了減少封堵器的能量損失，需要對封堵前后的壓力差進行控制，以降低封堵時的壓力脈動。由于封堵器尾部流場變化劇烈，因此選擇封堵器尾部截面壓力的絕對平均值作為觀測量。

管道中心截面的壓力分布如圖4所示，未開始封堵時管道上游壓力變化劇烈，呈現(xiàn)明顯的梯度分布。封堵器前后端面的壓差較大，為封堵器在管道內(nèi)前進提供了動力。封堵即將結(jié)束時管內(nèi)流場的壓力分布比較均勻，封堵器前端面所受壓力達到20 MPa，會對封堵器產(chǎn)生一定的沖擊，低壓區(qū)域主要集中在封堵器尾部，這與流場中形成的漩渦有關(guān)。封堵即將結(jié)束時的壓力遠遠大于未開始封堵時的壓力，且低壓區(qū)范圍逐漸增大，壓力脈動明顯，使得封堵器能量損失嚴重，封堵時產(chǎn)生振動。

圖4 管道中心截面的壓力分布

圖5為封堵器表面的壓力分布圖。隨著封堵比例的增大，封堵器受力逐漸增大，且表面各部分的壓力分界比較明顯。封堵器尾部與前端面的壓力差距逐漸增大，封堵即將結(jié)束時前端面的壓力要比尾部壓力高1個數(shù)量級，這是由于封堵進行過程中密封圈直徑增大，到達封堵器尾部的流量逐漸減小。封堵器在封堵前后的受力也存在較大差距，這會導(dǎo)致封堵器狀態(tài)不穩(wěn)定，封堵時出現(xiàn)振動現(xiàn)象，影響封堵的穩(wěn)定性。

圖5 封堵器表面的壓力分布

管道中心截面的速度分布如圖6所示。可以看出，未開始封堵時管內(nèi)流場的速度變化劇烈，在0～6 m/s波動。由于封堵器的阻隔，使封堵器前端面出現(xiàn)駐點，流體動能損失較大。流體在高壓作用下沿封堵器與管壁之間的狹長縫隙高速流過，到達下游后速度降低，且下游區(qū)域流體流動情況復(fù)雜，封堵器尾部附件有明顯的回流和漩渦現(xiàn)象，這會對封堵器造成一定的沖擊。封堵即將結(jié)束時管內(nèi)流場速度分布均勻，流速較低，此時管內(nèi)流體流動較為平穩(wěn)。

圖6 管道中心截面的速度分布

3 基于蒙特卡羅算法的管內(nèi)智能封堵器減振結(jié)構(gòu)優(yōu)化設(shè)計

通過分析管內(nèi)流場的分布規(guī)律，發(fā)現(xiàn)流體的流動情況與擾流板的面積及所處的角度有關(guān)，因此可以調(diào)整擾流板的角度和面積來使封堵器在管道內(nèi)部工作狀態(tài)達到最佳。

3.1 試驗設(shè)計與結(jié)果分析

針對管內(nèi)智能封堵器擾流板結(jié)構(gòu)的優(yōu)化問題，選擇擾流板的角度和面積為自變量，以封堵前后的壓力差為響應(yīng)量，利用CCD中心試驗法設(shè)計試驗方案。封堵前后的壓力差可以反映出封堵過程中的壓力脈動，壓差越大，說明管內(nèi)流場的壓力脈動就越大，這會導(dǎo)致封堵器尾部的漩渦和回流現(xiàn)象增多，封堵器受到的沖擊力變大，在管內(nèi)工作時的能量損失增加。因此選擇封堵前后的壓力差來衡量管內(nèi)流場的穩(wěn)定性，變量的范圍和響應(yīng)面試驗方案分別如表2～3所示。

表2 設(shè)計變量

表3 CCD中心試驗方案

根據(jù)設(shè)計的試驗方案建立模型，進行數(shù)值模擬試驗，邊界條件與上述相同。在響應(yīng)面試驗中會出現(xiàn)多組中心點，其目的是提高誤差估計的準(zhǔn)確性[18-19]。由于本次試驗沒有外界環(huán)境影響，因此中心點得出的結(jié)果相同。從試驗結(jié)果可以看出，各組試驗物理量的分布情況基本相同，加入擾流板后的壓力差相比于原始模型明顯降低，說明擾流板可以有效降低封堵時的壓力脈動，減少能量損失。但擾流板處于不同角度和面積時，封堵前后的壓力差也存在較大差異，這表明擾流板的角度和面積對封堵器的減振效果存在較大影響，通過優(yōu)化擾流板的參數(shù)對于封堵器的減振具有較高的可行性。

封堵前后的壓力差與擾流板角度、面積的關(guān)系如圖7所示。由圖7可以看出，關(guān)系曲面大致呈拋物面的形狀。從角度來看，當(dāng)擾流板面積較小時，隨著角度的增大，壓力差逐漸下降；當(dāng)擾流板面積較大時，隨著角度的增大，壓力差呈上升趨勢。從面積來看，當(dāng)擾流板角度較小時，隨著面積的增大，壓力差先減小后增大；當(dāng)擾流板角度較大時，隨著面積的增大，壓力差呈上升趨勢。因此，擾流板的角度和面積對封堵前后壓力差的影響較為顯著。

圖7 壓力差與擾流板角度和面積的關(guān)系

利用響應(yīng)面法對數(shù)值模擬試驗得出的數(shù)據(jù)進行擬合，得到Δp與α、S之間的數(shù)學(xué)模型。

Δp=1.329 63×105-227.424 63×α-
16.853 44×S+0.059 727×α×S+
0.178 89×α2+1.867 11×10-3×S2

(7)

式中：α為擾流板的翻轉(zhuǎn)角度；S為擾流板的面積；Δp為封堵前后封堵器尾部截面的壓力差。

通過對數(shù)學(xué)模型進行方差分析，發(fā)現(xiàn)其與試驗數(shù)據(jù)的擬合程度為97.92 %，說明該模型可以較好地預(yù)測試驗結(jié)果，具有較高的可靠性，可以作為目標(biāo)函數(shù)進行封堵器減振結(jié)構(gòu)的優(yōu)化。

3.2 管內(nèi)智能封堵器擾流板參數(shù)優(yōu)化

蒙特卡羅算法是一種統(tǒng)計模擬方法，這種方法可以在隨機采樣的基礎(chǔ)上得到近似的結(jié)果，并且采樣數(shù)越多，得到正確結(jié)果的概率就越大。它的基本思想是用某個事件出現(xiàn)的頻率來估計此事件的概率[20-21]，主要步驟為：首先構(gòu)造概率過程，對于隨機問題來說，需要準(zhǔn)確描述該過程，對于不是隨機性質(zhì)的問題，則需要構(gòu)造一個隨機的概率過程，將其轉(zhuǎn)化為隨機性質(zhì)的問題；其次就是要實現(xiàn)從概率分布中抽樣，已知的概率模型都是由各種概率分布組成的，而隨機變量就是蒙特卡羅方法的基本手段。對于均勻分布，隨機數(shù)就是隨機變量，從變量范圍內(nèi)產(chǎn)生隨機數(shù)就是從均勻分布抽樣的問題，隨機數(shù)是蒙特卡羅方法的基本工具；最后在完成抽樣的前提下，需要確定一個變量作為問題求解的對象，通過觀察結(jié)果，得到問題的解。蒙特卡羅算法適應(yīng)性強，可以求解模型復(fù)雜和計算困難的問題，方法簡單，便于編制程序，計算準(zhǔn)確性較高，且計算速度快。

管內(nèi)智能封堵器擾流板參數(shù)的優(yōu)化相當(dāng)于求解多變量、有約束的非線性最小化問題，其目標(biāo)函數(shù)為：

minΔp=f(α,S)

(8)

約束條件為30°≤α≤150°，1 500 mm2≤S≤6 000 mm2。采用蒙特卡羅算法進行最優(yōu)化參數(shù)求解，其優(yōu)化流程如圖8所示，利用Python語言進行算法程序的設(shè)計。將α和S看作指定范圍內(nèi)的均勻分布，生成隨機數(shù)作為程序的試驗點，實現(xiàn)概率分布中的抽樣。不斷代入隨機數(shù)求解目標(biāo)函數(shù)值，并將求得的值與上一次的值進行對比，保留最小函數(shù)值。為了保證計算的精度，迭代次數(shù)設(shè)置為108。迭代完成后，得到的最優(yōu)解參數(shù)為α=150 °，S=2 117 mm2。代入到式(7)后，得到此時封堵前后的壓力差為Δp=84 529.86 Pa。

4 優(yōu)化模型與原始模型對比

根據(jù)得到的優(yōu)化參數(shù)建立模型，將封堵器擾流板角度設(shè)置為150 °，面積為2 117 mm2，設(shè)置相同的邊界條件進行數(shù)值模擬，計算后的響應(yīng)量明顯低于各個試驗組，說明優(yōu)化后的模型可以很好地實現(xiàn)封堵器的減振。按照同樣的方法建立管內(nèi)智能封堵器的原始模型，將優(yōu)化模型與原始模型進行對比，結(jié)果如表4所示。可以看出，優(yōu)化模型的壓力差比原始模型小64.1%，最大壓力下降了47.8%，而最大湍動能比原始模型減少了將近1倍。這表明利用蒙特卡羅方法優(yōu)化的模型可以有效減少封堵過程中的壓力脈動，降低管內(nèi)流體對封堵器的沖擊，使流場變得更加穩(wěn)定，同時減少了漩渦和回流現(xiàn)象，降低封堵器的能量損失，達到減振的目的。

圖8 蒙特卡羅算法流程

參數(shù)封堵前后的壓力差/Pa最大壓力/Pa最大湍動能/(m2·s-2)原始模型130 31418 619 900420.56優(yōu)化模型79 42112 595 100251.81

圖9～10分別為原始模型與優(yōu)化模型封堵前后的壓力分布圖，從圖中看出優(yōu)化模型管內(nèi)流場的壓力比原始模型小，且壓力的分布情況大致相同，未封堵時優(yōu)化模型的壓力梯度范圍相對較小，因此對封堵器的影響相比于原始模型要小。封堵即將結(jié)束時管內(nèi)流場的壓力分布都比較均勻，但從數(shù)值上看，優(yōu)化模型封堵器前端面受力比原始模型小48 %左右，說明蒙特卡羅算法優(yōu)化后的模型對封堵器產(chǎn)生的壓力小，有利于保證封堵的穩(wěn)定性。

為了更加清晰地對比原始模型與優(yōu)化模型管內(nèi)智能封堵器在流場中的受力情況，繪制了不同封堵進度下封堵器尾部所受的壓力對比曲線，如圖11所示。由圖11可以看出，隨著封堵的進行，封堵器受力不斷增大。從初始狀態(tài)到封堵70 %過程中，封堵器受力增長的速度比較緩慢，且原始模型與優(yōu)化模型受力相差不多；封堵70 %以后，封堵器受力的增長速度變快，優(yōu)化模型與原始模型的差距逐漸增大，表明隨著封堵過程的進行，優(yōu)化模型表現(xiàn)出的作用越來越明顯，擾流板對于減少封堵器受力的作用也越來越大，優(yōu)化后的模型可以減少能量損失，實現(xiàn)封堵器的減振。

圖9 原始模型壓力分布云圖

圖10 優(yōu)化模型壓力分布云圖

圖11 原始模型與優(yōu)化模型封堵器受力對比曲線

5 結(jié)論

1) 為了解決管內(nèi)智能封堵器工作過程中的振動問題，在封堵器尾部結(jié)構(gòu)設(shè)計3塊可折疊擾流板。通過數(shù)值模擬發(fā)現(xiàn)該結(jié)構(gòu)對管內(nèi)流場有一定的改善。

2) 利用CCD中心試驗方法設(shè)計試驗方案，建立擾流板角度、面積與封堵前后壓力差的數(shù)學(xué)模型。采用蒙特卡羅算法將壓力差作為目標(biāo)函數(shù)，以角度和面積作為參數(shù)變量，從中進行抽樣、產(chǎn)生隨機數(shù)。通過不斷迭代求解出使目標(biāo)函數(shù)最小的參數(shù)，并根據(jù)得到的參數(shù)建立封堵器模型，進行數(shù)值模擬。

3) 優(yōu)化后的模型相比原始模型封堵前后的壓力差明顯降低，且最大壓力和最大湍動能也顯著減少，說明用蒙特卡羅算法優(yōu)化模型的管內(nèi)流場的壓力脈動降低，封堵器的能量損失減少，有效改善了管內(nèi)流場，使封堵過程更加穩(wěn)定，對于減少封堵器的振動具有重要意義。