可壓縮兩相流固耦合模型的間斷Galerkin 有限元方法1)

馬天然 沈偉軍 *?,2) 劉衛(wèi)群 Xu Hao

* (中國(guó)礦業(yè)大學(xué)力學(xué)與土木工程學(xué)院,江蘇徐州 221116)

? (中國(guó)礦業(yè)大學(xué)深部巖土力學(xué)與地下工程國(guó)家重點(diǎn)實(shí)驗(yàn)室,江蘇徐州 221116)

** (中國(guó)科學(xué)院力學(xué)研究所,北京 100190)

?? (中國(guó)科學(xué)院大學(xué)工程科學(xué)學(xué)院,北京 100049)

*** (勞倫斯伯克利國(guó)家實(shí)驗(yàn)室地球科學(xué)部,美國(guó)伯克利 94720)

引言

深入探索多孔介質(zhì)中多相流動(dòng)和固體變形特征對(duì)精準(zhǔn)認(rèn)識(shí)地下資源開發(fā)利用,例如石油勘探開發(fā)、煤層氣和頁(yè)巖氣開采、二氧化碳地質(zhì)封存等領(lǐng)域,具有重要的工程意義[1-4].在流體輸送過(guò)程中局部壓力和飽和度梯度下的滲流作用改變了孔隙結(jié)構(gòu)特征,打破了儲(chǔ)層內(nèi)部的水力平衡;固體空間結(jié)構(gòu)的變化影響了儲(chǔ)層內(nèi)部流體的運(yùn)移路徑和流通能力,具體表現(xiàn)為多孔介質(zhì)孔隙度和滲透率等物性參數(shù)的動(dòng)態(tài)響應(yīng).發(fā)展和建立兩相滲流與固體變形耦合模型及其相應(yīng)的數(shù)值方法是準(zhǔn)確描述和闡明儲(chǔ)層內(nèi)水力行為的有效手段[5-7].

針對(duì)兩相滲流和固體變形耦合模型的數(shù)值方法主要包括有限元法[8-10]、有限體積法[11-12]以及混合計(jì)算方法[13-14]等.對(duì)于兩相流問(wèn)題的離散和數(shù)值求解,傳統(tǒng)有限元能夠確保流體域內(nèi)整體的質(zhì)量平衡,但是無(wú)法實(shí)現(xiàn)流體在跨越單元間的局部質(zhì)量平衡,因此不適用于求解對(duì)流占主導(dǎo)的非線性不混溶問(wèn)題;有限體積法在確保局部流量平衡的同時(shí),具有較高的計(jì)算效率,因此在工業(yè)界和工程尺度模擬中得到了廣泛應(yīng)用.在通常情況下,有限體積方法的計(jì)算結(jié)果僅能滿足低階有效精度;混合計(jì)算方法是指對(duì)流體和力學(xué)控制方程分別利用不同的數(shù)值方法進(jìn)行求解,該方法集成了各種數(shù)值方法的優(yōu)勢(shì).Rutqivst[15]在耦聯(lián)流體軟件TOUGH2 和固體變形軟件FLAC3D的基礎(chǔ)上,開發(fā)了在油氣、地?zé)?、可燃冰開采以及二氧化碳地質(zhì)封存等領(lǐng)域廣泛使用的TOUGH-FLAC耦合模擬器.在此模擬器中,負(fù)責(zé)流體計(jì)算的TOUGH2 軟件采用的是有限體積法,FLAC3D 采用有限差分法計(jì)算固體變形,此方法在網(wǎng)格劃分和數(shù)據(jù)傳遞效率等方面具有一定的局限性.

如何在刑事立法政策的制定過(guò)程中辯證地認(rèn)識(shí)民意，正確引導(dǎo)民意、有效利用民意是一個(gè)需要亟待解決的問(wèn)題。然而民意的定義與范圍十分抽象、復(fù)雜，甚至真?zhèn)坞y辨，一旦研判失誤，極易誤導(dǎo)決策者，甚至可能造成社會(huì)混亂。因而對(duì)民意的考量必須加以合理的規(guī)制。

間斷伽遼金有限元(discontinuous Galerkin finite element method,DG-FEM)結(jié)合了有限體積與有限元兩種方法的優(yōu)點(diǎn),具有高階精度和局部單元物理守恒性的特點(diǎn),容易實(shí)現(xiàn)局部網(wǎng)格加密和各單元多項(xiàng)式的單獨(dú)選取[16-19].近些年來(lái),DG-FEM 在流體和彈性力學(xué)領(lǐng)域得到了廣泛的關(guān)注和持續(xù)的發(fā)展[20-22].Klieber 和Riviere[23]在解耦兩相滲流微分方程組的基礎(chǔ)上,提出了不可壓縮兩相流自適應(yīng)間斷有限元方法.Epshteyn 和Rivière[24]提出了多孔介質(zhì)不可壓縮兩相流完全隱式方程的不連續(xù)伽遼金方法,該方法在結(jié)構(gòu)化和非結(jié)構(gòu)化網(wǎng)格計(jì)算中都具有較好的魯棒性.李子然和吳長(zhǎng)春[25]、陳韻騁等[26]建立了基于線彈性力學(xué)問(wèn)題的局部間斷有限元方法.在流固耦合方面,Liu 等[27-28]提出了連續(xù)和間斷伽遼金耦合框架,利用自適應(yīng)加罰方法提高間斷伽遼金有限元方法在模擬大型多孔彈性問(wèn)題時(shí)的適用性和計(jì)算效率,并將此方法推廣至多孔介質(zhì)水力壓裂模擬中.目前關(guān)于間斷伽遼金有限元方法在計(jì)算可壓縮兩相流固模型方面的研究鮮有報(bào)道,還有待進(jìn)一步深入研究.

隨著社會(huì)經(jīng)濟(jì)的不斷發(fā)展以及人們生活水平的不斷提高，對(duì)道路橋梁工程的施工質(zhì)量的要求也越來(lái)越嚴(yán)格。為了滿足社會(huì)發(fā)展的需要，以及人們交通出行的需求，應(yīng)不斷提高道路橋梁工程的建設(shè)質(zhì)量。因此，保證道路橋梁工程原材料的質(zhì)量成為保障工程建設(shè)質(zhì)量的基礎(chǔ)與前提。

利用上述類似的推導(dǎo)過(guò)程和方式,可以得到式(2)的間斷伽遼金弱形式,如下所示

1 兩相流固耦合控制方程

1.1 兩相滲流方程

假設(shè)模型中濕潤(rùn)相和非濕潤(rùn)相為不混溶流體,兩者在多孔介質(zhì)中的流速滿足達(dá)西定律.在考慮毛細(xì)壓力、重力效應(yīng)和固體變形的基礎(chǔ)上,等溫可壓縮兩相滲流控制方程可表示為[2,29]

式中,下標(biāo)β 為w和 n w分別表示濕潤(rùn)相和非濕潤(rùn)相,Sβ為飽和度, ρβ和μβ分別為流體密度和黏度,Cβ為流體壓縮系數(shù),pβ為流體壓力,k為多孔介質(zhì)絕對(duì)滲透率張量,krβ為相對(duì)滲透率,Qβ為源匯項(xiàng), εv為體應(yīng)變,pc和 |pc|分別為毛細(xì)壓力和毛細(xì)壓力對(duì)濕相飽和度的導(dǎo)數(shù).

針對(duì)學(xué)校知識(shí)資產(chǎn)的存在形態(tài)與管理目標(biāo)，結(jié)合信息系統(tǒng)分析與設(shè)計(jì)、文檔一體化管理、關(guān)系數(shù)據(jù)庫(kù)系統(tǒng)等原理，高校機(jī)構(gòu)內(nèi)知識(shí)資源綜合管理平臺(tái)的結(jié)構(gòu)設(shè)計(jì)，整體上應(yīng)該包括系統(tǒng)用戶管理、各職能部門應(yīng)用子系統(tǒng)管理、知識(shí)檔案資源數(shù)據(jù)庫(kù)管理、知識(shí)庫(kù)增值服務(wù)管理等四大管理模塊。

為了求解方程,需要補(bǔ)充如下輔助公式

式中,pe為毛細(xì)管吸入壓力,Se為有效飽和度, λ 為孔徑分布指數(shù),Srw和Srnw分別為濕潤(rùn)相和非濕潤(rùn)相的殘余飽和度.

此外,選取Genuchten?Mualem 相對(duì)滲透率模型描述濕潤(rùn)相和非濕潤(rùn)相的流通能力,表達(dá)式如下所示[30]

模型的初始條件為

模型中Dirichlet 和Neumann 邊界條件如下所示

進(jìn)入防汛物資管理界面可對(duì)各倉(cāng)庫(kù)進(jìn)行防汛物資庫(kù)存情況查詢，網(wǎng)上申請(qǐng)、審批發(fā)放管理，對(duì)購(gòu)置入庫(kù)、維修保養(yǎng)、出入庫(kù)時(shí)間記錄列表上傳、下載進(jìn)行物資倉(cāng)庫(kù)管理。

1.2 兩相滲流方程的間斷伽遼金弱形式

圖1 為模擬幾何區(qū)域、邊界條件和網(wǎng)格示意圖.在計(jì)算域 Ω中,邊界 ? Ω由不與內(nèi)部單元接觸的外邊界 ? Ωo和內(nèi)部相鄰單元共享的內(nèi)邊界 ? Ωi兩部分組成.由圖1 可知,內(nèi)邊界 ? Ωi為相鄰單元E+和E?共同擁有;外邊界?Ωo由Dirichlet 邊界?ΩoD和Neumann 邊界 ? ΩoN兩者構(gòu)成,即 ? Ωo=?ΩoD∪?ΩoN.設(shè)是區(qū)域 Ω的正則剖分, ?E=?Ωo∪?Ωi表示為所有邊界的集合,則間斷有限元空間定義為

圖1 幾何區(qū)域、邊界條件和網(wǎng)格示意圖Fig.1 Domain with boundary conditions and mesh skeleton

式中,離散空間pr(E)為單元E上不超過(guò)r次的多項(xiàng)式的集合.

為了推導(dǎo)兩相滲流方程的間斷伽遼金弱形式,首先在式(1)等號(hào)的兩側(cè)乘以試函數(shù)并在單元計(jì)算域E內(nèi)進(jìn)行積分,可得到

利用分部積分和散度定理,式(15)中等號(hào)左側(cè)第二項(xiàng)可表示為

方程式(19),式(20)和式(27)是多孔介質(zhì)中可壓縮兩相滲流和固體變形的弱形式.本文利用通用有限元軟件COMSOL Multiphysics 內(nèi)置的Weak Form PDE 模塊求解可壓縮兩相流固耦合模型.基于二次不連續(xù)拉格朗日形函數(shù)構(gòu)建對(duì)解域的間斷伽遼金有限元逼近,選取向后差分公式(backward differentiation formula,BDF)隱式求解器.其中,BDF 求解器采用自適應(yīng)時(shí)間步長(zhǎng)算法和可變離散化階數(shù),離散化的BDF 精度從1～5 不等.基于LU 分解,采用MUMPS (multifrontal massively parallel sparse)直接求解器對(duì)線性系統(tǒng)進(jìn)行求解,并選取完全耦合的求解器(fully coupled solver)將線性求解器應(yīng)用于牛頓單步法的非線性問(wèn)題.默認(rèn)情況下,初始時(shí)間步長(zhǎng)設(shè)置為結(jié)束時(shí)間的0.1%,相對(duì)容差為0.01.

水還沒(méi)消下去，楊小水卻堅(jiān)持要回去。一路上看到的樹，樹梢上都掛滿了水草。第二天進(jìn)入遂平境內(nèi)，連樹都少見了。大的多伏在地上，小的，連根都拔走了。老鼠都圓滾滾的，像小孩子玩的皮球，也不怕人，在地上緩緩地滾動(dòng)。鐵路線這邊的路溝里，是她這輩子見過(guò)尸體最多的地方。層層疊疊地摞著，不計(jì)其數(shù)。附近的樹枝上落滿了蒼蠅，黑壓壓的，把樹都?jí)簭澚恕?/p>

將式(18)代入式(17)中,結(jié)合穩(wěn)定項(xiàng)、加罰函數(shù)項(xiàng)以及相應(yīng)邊界條件,可得

基于上述認(rèn)識(shí),本文首先建立可壓縮兩相流固耦合控制方程,其中包括考慮毛細(xì)壓力和重力作用的可壓縮兩相滲流方程以及線性多孔彈性方程.在此基礎(chǔ)上,推導(dǎo)出流固耦合方程的間斷伽遼金弱形式;然后,通過(guò)一維Terzaghi 流固問(wèn)題驗(yàn)證模型的準(zhǔn)確性;最后,分別開展二維平面應(yīng)變條件下和考慮重力效應(yīng)的三維流固數(shù)值算例,分析流體壓力、飽和度以及固體應(yīng)力、位移的分布特征,探討加罰因子對(duì)模擬結(jié)果的影響.

式(19)和式(20)為考慮固體變形作用下可壓縮兩相滲流方程的間斷伽遼金有限元方法的弱表達(dá)形式.

大學(xué)生作為獨(dú)特的社會(huì)群體，其思想活躍、充滿創(chuàng)意意識(shí)，并由此形成了別具特色的校園文化?；诖耍咝５慕逃虒W(xué)中，應(yīng)將社會(huì)主義核心價(jià)值觀融入到校園文化中來(lái)，利用校園文化活動(dòng)，為大學(xué)生的思想價(jià)值培養(yǎng)提供良好的外部環(huán)境。例如在校園文化建設(shè)中，學(xué)校可以組織社會(huì)主義核心價(jià)值觀的專題講座，做好理論宣傳；組織社區(qū)志愿活動(dòng)，讓學(xué)生在社區(qū)服務(wù)中強(qiáng)化社會(huì)責(zé)任感；組織專業(yè)實(shí)習(xí)，將學(xué)校教育與職業(yè)需求相結(jié)合，提高大學(xué)生的自我認(rèn)識(shí)；參與社會(huì)調(diào)查，從政治、文化、經(jīng)濟(jì)、生態(tài)等多角度了解社會(huì)，讓大學(xué)生從書本中走出來(lái)，提高社會(huì)實(shí)踐能力。

（4）止水帶壓板采用“凵”型鋼板，在螺帽和壓板之間設(shè)置彈性墊圈，調(diào)整定位后擰緊。螺帽以上緊為準(zhǔn)，原則上緊固力為3kg，一般不宜超過(guò)5kg，但不應(yīng)小于2.5kg。螺栓外露部分采用塑料套保護(hù)。此種鋼板剛度較“一”型成倍增加，擰緊螺栓時(shí)不會(huì)因壓板變形而產(chǎn)生漏水空洞。

1.3 固體變形方程及其弱形式

假設(shè)多孔介質(zhì)固體變形滿足彈性小變形條件,那么固體變形的平衡方程、幾何方程、本構(gòu)方程以及邊界條件如下所示

式中,總應(yīng)力 σ =σ′?αIp,則式(21)可改寫

2.3.1 產(chǎn)前遺傳咨詢及產(chǎn)前診斷情況 截止到2016年11月，通過(guò)電話隨訪及門診回訪232對(duì)夫妻進(jìn)行了進(jìn)一步的遺傳咨詢并對(duì)丈夫相關(guān)基因進(jìn)行了檢測(cè)。其中同基因攜帶為33例，其中2例暫未懷孕，18對(duì)夫妻經(jīng)知情同意選擇進(jìn)一步產(chǎn)前診斷，其余拒絕產(chǎn)前診斷。

式中, σ′為有效應(yīng)力張量, ε為應(yīng)變張量,u為位移張量,f為體積力張量, α 為Biot 常數(shù),D為4 階彈性張量,I為單位張量,平均孔隙壓力p=Snwpnw+Swpw.

本節(jié)利用間斷伽遼金方法開展在二維平面應(yīng)變條件下兩相滲流和固體變形耦合數(shù)值計(jì)算.在本算例中,分別選取水和氫氣作為濕潤(rùn)相和非濕潤(rùn)相.模型的幾何尺度、初始條件、邊界條件和監(jiān)測(cè)線如圖4 所示.模型的長(zhǎng)度和寬度均設(shè)置為10.00 m.多孔介質(zhì)的彈性模量和泊松比分別為10.00 GPa 和0.30;滲透率和孔隙度分別設(shè)定為 1 .00×10?14m2和0.20.模擬區(qū)域的初始?xì)怏w壓力為4.00 MPa,初始水飽和度為0.80.模擬左邊界和下邊界為位移約束,在上邊界和右邊界分別施加 ? 1.00 × 107MPa (負(fù)號(hào)表示壓應(yīng)力方向)的垂直和水平應(yīng)力.模型左下角為注入邊界,設(shè)置注入的氫氣壓力和水飽和度分別為5.00 MPa 和0.80;右上角為出口邊界,出口氣體壓力和水飽和度設(shè)置為4.00 MPa 和0.20.

1.4 數(shù)值方法的實(shí)現(xiàn)

將式(16)代入式(15),累加所有單元區(qū)域,可得

2 模型驗(yàn)證

圖2 描述了一維Terzaghi 固結(jié)問(wèn)題的幾何示意圖以及相應(yīng)的邊界條件.為了確?？紫秹毫ρ厣舷聝啥吮３譃榱?設(shè)置模型上下兩端為排水邊界.在模型上邊界瞬間施加均布載荷q=?5.00 MPa.本模型長(zhǎng)度 2h=10.00 m,其他參數(shù)的具體取值如下[2]:剪切模量為30.00 GPa,壓縮模量為50.00 GPa,多孔介質(zhì)滲透率為 1 .00×10?12m2,孔隙度為0.25,流體壓縮系數(shù)為 1 .00×10?12Pa?1, 流體黏度為 1 .00×10?3Pa·s,加罰因子 δf和 δs分別設(shè)置為1.00 和10.00[31].

干法電選是利用粉煤灰在高壓電場(chǎng)作用下，因灰與炭導(dǎo)電性能不同而進(jìn)行的分離。粉煤灰是非導(dǎo)體物料，炭粒是良好的導(dǎo)體物料，在圓形電暈電場(chǎng)中，當(dāng)粉煤灰獲得電荷后，炭粒因?qū)щ娦阅芰己茫芸斓貙⑺@電荷通過(guò)圓筒帶走，在重力慣性離心力作用下，脫離圓筒表面，被拋入導(dǎo)體產(chǎn)品槽；而非導(dǎo)體的粉煤灰所獲電荷在表面釋放速度較慢，故在電場(chǎng)力作用下，吸收在圓筒表面上，被旋轉(zhuǎn)圓筒帶到后部，由卸料毛刷排入非導(dǎo)體產(chǎn)品槽中，從而達(dá)到灰炭分離的效果。

隨著社會(huì)經(jīng)濟(jì)和醫(yī)藥產(chǎn)業(yè)的發(fā)展及公眾日益增長(zhǎng)的健康需求，醫(yī)藥新產(chǎn)品、新劑型層出不窮，例如中藥材的新型產(chǎn)地加工或炮制加工品種、中藥配方顆粒、中藥精制飲片超微細(xì)粉加工等，相應(yīng)的產(chǎn)品質(zhì)量標(biāo)準(zhǔn)研發(fā)需求旺盛。藥品檢驗(yàn)機(jī)構(gòu)作為藥品監(jiān)管和產(chǎn)業(yè)發(fā)展的重要技術(shù)支撐，承擔(dān)著大量的注冊(cè)檢驗(yàn)、監(jiān)督抽檢、仿制藥一致性評(píng)價(jià)復(fù)核檢驗(yàn)、質(zhì)量標(biāo)準(zhǔn)制修訂和補(bǔ)充檢驗(yàn)方法研制等檢驗(yàn)檢測(cè)和技術(shù)研究任務(wù)，是新時(shí)代市場(chǎng)監(jiān)管和醫(yī)藥產(chǎn)業(yè)發(fā)展不可或缺的技術(shù)力量。如何充分發(fā)揮藥檢機(jī)構(gòu)在質(zhì)量標(biāo)準(zhǔn)研究方面的技術(shù)優(yōu)勢(shì)，為市場(chǎng)監(jiān)管、產(chǎn)業(yè)發(fā)展和推動(dòng)“放管服”提供技術(shù)支撐，值得深入研究。

圖2 Terzaghi 固結(jié)問(wèn)題幾何示意圖Fig.2 Sketch of Terzaghi’s consolidation problem

圖3 為Terzaghi 固結(jié)問(wèn)題的理論解[32]和DG有限元數(shù)值解的比較.結(jié)果顯示,DG 有限元計(jì)算得到的數(shù)值解與理論解吻合較好,進(jìn)而驗(yàn)證了間斷伽遼金方法在求解流固耦合方程時(shí)的可行性和有效性.

在間斷伽遼金方法中,變量在單元內(nèi)邊界 ? Ωi處滿足不連續(xù)設(shè)定,因此式(17)中等號(hào)左側(cè)第三項(xiàng)可以改寫為

圖3 Terzaghi 固結(jié)問(wèn)題理論解和數(shù)值解的比較Fig.3 Comparison of the analytical and numerical results for Terzaghi’s consolidation problem

3 算例分析

3.1 二維平面兩相流固算例

利用類似式(19)和式(20)的推導(dǎo)過(guò)程,結(jié)合相應(yīng)的邊界條件,可以得到固體變形方程的弱形式,如下所示

圖4 模擬幾何尺度、初始和邊界條件Fig.4 Simulation geometrical configuration with boundary and initial conditions

數(shù)值模擬分兩個(gè)階段開展.首先,開展儲(chǔ)層力學(xué)和流體平衡的穩(wěn)態(tài)計(jì)算.然后,將第一步計(jì)算得到的應(yīng)力和孔壓等水力信息作為含水層氫氣注入模擬算例的初始條件.整個(gè)模擬時(shí)間為1500 s,其他變量的具體數(shù)值如表1 所示.

圖5 和圖6 分別給出了注氣第1500 s 氣體壓力pnw、水飽和度Sw、水 平位移u和水平有效應(yīng) 力的空間分布云圖以及變量沿監(jiān)測(cè)線的分布圖.由圖5和圖6 可知,隨著氫氣的持續(xù)注入,氣體壓力影響范圍擴(kuò)大至距注入邊界7 m 左右,并引起固體膨脹變形且變形影響范圍擴(kuò)散至邊界處.究其原因是因?yàn)樵诙嗫讖椥越橘|(zhì)中位移或者應(yīng)變可傳播到壓力前沿之前.受毛細(xì)壓力和氣水兩相物理屬性差異的共同影響,氫氣主要聚集于注入邊界附近.在多孔介質(zhì)內(nèi)孔壓在遠(yuǎn)離注入邊界的過(guò)程中呈現(xiàn)遞減趨勢(shì),而有效應(yīng)力沿監(jiān)測(cè)線增加.

圖6 第1500 s 時(shí)氣體壓力 pnw、水飽和度 Sw、水平位移 u 和水平有效應(yīng)力沿監(jiān)測(cè)線的分布Fig.6 The profiles of gas pressure pnw,water saturation S w ,horizontal displacement u and horizontal effective stress along the monitoring line at t = 1500 s

3.2 加罰因子的影響

圖7 給出了 δs=0.10,1.00 和10.00 時(shí)水平位移u和水平有效應(yīng)力沿監(jiān)測(cè)線的分布.圖8 給出了δs=0.10時(shí)水平位移u和水平有效應(yīng)力 σ′x的分布云圖.圖9 給出了 δf=0.01,0.10,1.00 和10.00 時(shí)氣體壓力pnw、水飽和度Sw沿監(jiān)測(cè)線的分布.結(jié)果表明,不同加罰因子條件下,變量u,,pnw和Sw沿著監(jiān)測(cè)線的變化趨勢(shì)大體保持一致.在計(jì)算固體變形時(shí),計(jì)算結(jié)果對(duì)加罰因子的選取更為敏感.當(dāng) δs= 0.10時(shí),應(yīng)力和位移均在局部出現(xiàn)明顯的波動(dòng);當(dāng) δf=0.01 時(shí),飽和度在小范圍內(nèi)發(fā)生輕微浮動(dòng).這主要由于加罰因子 δs和 δf的降低減少了自變量位移、壓力以及飽和度在單元處弱連續(xù)性的約束,造成了變量局部波動(dòng)的現(xiàn)象.通過(guò)適當(dāng)提高加罰因子可以加強(qiáng)變量在跨越單元交界面處的連續(xù)性,但是較高的加罰值會(huì)引發(fā)病態(tài)條件矩陣的可能性.

圖7 第1500 s 時(shí)不同加罰因子 δs條件下水平位移 u 和水平有效應(yīng)力 沿監(jiān)測(cè)線的分布Fig.7 The profiles of horizontal displacement u and horizontal effective stress along the monitoring line at t = 1500 s with different values of penalty parameter δs

圖8 第1500 s 時(shí) δs=0.01水平位移 u和水平有效應(yīng)力 分布云圖Fig.8 The distribution of horizontal displacement u and horizontal effective stress with δs=0.01 at t = 1500 s

圖9 第1500 s 時(shí)不同加罰因子 δf條件下氣體壓力pnw 、水飽和度Sw 沿監(jiān)測(cè)線(0,0.05)?(10,9.95)的分布Fig.9 The profiles of gas pressure pnw and water saturation S w along monitoring line at t = 1500 s with different values of δf

值得注意的是,受流體邊界條件、初始條件以及流體方程的非線性特征等因素的影響,傳統(tǒng)有限元方法在求解本節(jié)流體模型時(shí)無(wú)法收斂,此情況在一定程度上體現(xiàn)了間斷伽遼金有限元在求解可壓縮兩相流固方程上的優(yōu)越性.本節(jié)將傳統(tǒng)有限元的計(jì)算結(jié)果作為基準(zhǔn)方案,通過(guò)比較間斷伽遼金有限元和傳統(tǒng)有限元兩者計(jì)算得到的水平位移和水平有效應(yīng)力之間的相對(duì)誤差,分析加罰因子 δs對(duì)計(jì)算結(jié)果的影響.水平位移和水平有效應(yīng)力的計(jì)算誤差erru和errs的表達(dá)式為[33]

式中,l為幾何模型尺度分別為有限元方法計(jì)算得到的水平位移和水平有效應(yīng)力, ΔuCG和分別為間斷有限元方法計(jì)算得到的水平位移和水平有效應(yīng)力.

圖10 給出了不同加罰因子 δs條件下間斷伽遼金有限元和傳統(tǒng)有限元方法計(jì)算得到的水平位移和水平有效應(yīng)力的比較.結(jié)果表明,隨著 δs的增加,間斷有限元與傳統(tǒng)有限元計(jì)算得到的水平方向的變形和應(yīng)力的相對(duì)誤差隨之減少,這說(shuō)明了加罰項(xiàng)的增加可以有效抑制有限元函數(shù)在跨越單元時(shí)的不連續(xù)性.

圖10 不同加罰因子 δs 條件下間斷有限元和傳統(tǒng)有限元方法計(jì)算水平位移和水平有效應(yīng)力的比較Fig.10 Comparison of the horizonal displacement and horizonal effective stress between DG and the CG results with different values of penalty parameter δs

3.3 考慮重力效應(yīng)的兩相流固三維算例

在本節(jié)中,將開展含水層注氫氣的三維數(shù)值模擬.在此算例中,側(cè)重分析流體重力效應(yīng)對(duì)模擬結(jié)果的的影響.如圖11 所示,三維算例的幾何尺度為100.00 m×2.00 m×10.00 m.模擬上邊界和右邊界為應(yīng)力邊界且tx=tz= ?1.00 × 107MPa,其余邊界均設(shè)定為位移約束.氫氣從左側(cè)邊界注入,注入率gnw=1.00 × 10?4kg/(m2·s),右側(cè)為出口邊界,邊界上的氣體壓力和水飽和度分別為4.00 MPa 和0.20,其余邊界則為不流動(dòng)邊界.整個(gè)注氣過(guò)程持續(xù)1 d,其他變量的數(shù)值如表1 所示.

圖11 三維模型幾何尺度以及模擬初始和邊界條件Fig.11 Geometrical configuration with boundary and initial conditions of 3D case

圖12 給出了注氣1 d 后氣體壓力pnw、氣體飽和度Snw、垂直位移w和水平有效應(yīng)力的分布云圖.圖13 為在考慮和不考慮重力效應(yīng)兩種情況下點(diǎn)a,b和c處氣飽和度Snw隨時(shí)間變化圖.由圖12 可知,氫氣注入1 d 后,孔壓沿著水平x方向傳播至距左邊界約70.00 m 處,儲(chǔ)層邊界處氣體飽和度和孔隙壓力分別增加至約0.52 和6.50 MPa,孔壓的增加引起垂直位移的提高和水平有效應(yīng)力的降低.在注入過(guò)程中,由于氫氣和水兩者密度的差異,引起氣體上浮且聚集于模型頂部,造成頂部點(diǎn)c處氣體飽和度要明顯高于底部點(diǎn)a處氣體飽和度.

圖12 注氣1 d 后氣體壓力 pnw、氣體飽和度 Snw、垂直位移 w 和水平有效應(yīng)力的分布圖Fig.12 The distribution of gas pressurepnw,gas saturation Snw ,vertical displacement wand horizontal effective stress after 1 day of injection

圖13 考慮和不考慮重力效應(yīng)算例中點(diǎn)a,b 和c 處氣飽和度 S nw 隨時(shí)間變化圖Fig.13 The temporal evolution of gas saturation S nw at points a,b and c in the cases with and without gravity

4 結(jié)論

本文建立了考慮毛細(xì)壓力和重力效應(yīng)的可壓縮兩相滲流與孔隙介質(zhì)變形耦合作用的力學(xué)模型.在此基礎(chǔ)上,分別推導(dǎo)出變形孔隙介質(zhì)中兩相滲流方程的非對(duì)稱內(nèi)罰伽遼金弱形式以及固體變形方程的非完全內(nèi)罰伽遼金弱形式.主要結(jié)論如下:

(1)通過(guò)對(duì)比一維Terzaghi 固結(jié)問(wèn)題的解析解和數(shù)值解,驗(yàn)證了間斷伽遼金有限元方法在求解流固耦合問(wèn)題方面的可行性和準(zhǔn)確性;

(2)開展了二維平面應(yīng)變條件下和考慮重力效應(yīng)作用的三維兩相流固數(shù)值算例.結(jié)果顯示,隨著氣體的注入,氣體飽和度和孔壓隨之增加,導(dǎo)致多孔介質(zhì)膨脹變形和有效應(yīng)力降低.由于重力的影響,氣體會(huì)上浮聚集于頂部邊界;

本文在GP100-C3高頻感應(yīng)淬火設(shè)備上，利用設(shè)計(jì)制作的仿形感應(yīng)器，對(duì)制動(dòng)器壓盤錐窩部位進(jìn)行感應(yīng)淬火，得到如下結(jié)論：

(3)在二維平面算例中,討論了加罰因子對(duì)計(jì)算結(jié)果的影響.結(jié)果表明,在相同網(wǎng)格條件下,加罰因子 δf和 δs的降低會(huì)導(dǎo)致流體和固體信息在局部出現(xiàn)不同程度的波動(dòng).固體變形參量對(duì)加罰因子 δs更為敏感,提高加罰因子 δs可以有效抑制位移和應(yīng)力在跨越單元時(shí)的不連續(xù)性.