蘇 巖，劉 彪，王正中，董 鵬，劉銓鴻，楊 誠

(1.陜西省引漢濟(jì)渭工程建設(shè)有限公司， 陜西 西安 710011;2.西北農(nóng)林科技大學(xué) 水利與建筑工程學(xué)院，陜西 楊凌 712100)

陜西省以秦嶺為界分為長江流域和黃河流域，處于黃河流域的關(guān)中地區(qū)和陜北地區(qū)存在季度缺水問題，且長江流域水資源豐富且不能被充分利用，目前的水資源配置無法滿足陜西省經(jīng)濟(jì)發(fā)展的需求[1]。為緩解關(guān)中渭河沿線城市和工業(yè)用水短缺的問題，引漢濟(jì)渭工程被提出，該工程建成后可為西安、咸陽、寶雞、渭南4個城市、沿渭河兩岸11個縣及6個工業(yè)園提供各類用水[2]。引漢濟(jì)渭二期輸配水工程是引漢濟(jì)渭工程的重要組成部分，由黃池溝配水樞紐、南干線、北干線及相應(yīng)的輸水支線組成[3-4]。二期工程從關(guān)中配水節(jié)點(diǎn)黃池溝起，輸水干線東到渭南市華州區(qū)，西到楊凌示范區(qū)，南到鄠邑區(qū)，北到富平縣，輸配水區(qū)域東西長約163 km，南北寬約84 km，總面積約1.4萬km2。圖1為引漢濟(jì)渭二期工程總平面布置示意圖。

圖1 引漢濟(jì)渭二期工程總平面布置示意圖

南干線黃池溝至灞河段全長102.07 km，由隧洞、倒虹、渡槽、管橋、箱涵及分退水設(shè)施組成。南干線工程主要穿越秦嶺嶺北中低山區(qū)、山前洪積扇區(qū)、河谷階地、黃土臺塬四大地貌單元，出露地層主要為下元古界寬坪群變質(zhì)巖、燕山期花崗巖、新近系沉積巖及四系松散堆積層。地下水以基巖裂隙水、第四系松散巖類孔隙水、新近系碎屑巖類裂隙孔隙水為主。黃午隧洞、神禾塬隧洞、少陵塬隧洞和白鹿塬隧洞這四段隧洞總長90.05 km，占南干線近90%。

神禾塬隧洞、少陵塬隧洞為黃土隧洞，白鹿塬隧洞為泥巖、砂礫巖隧洞。鑒于黃土隧洞及泥巖隧洞成洞條件均十分困難，且隧洞埋深不大，地表有村落房屋等建筑物，在高地下水位的飽和性黃土洞段及泥砂巖洞段，常規(guī)鉆爆法施工風(fēng)險極高，因此，擬采用盾構(gòu)法施工。引漢濟(jì)渭二期工程中神禾塬隧洞、少陵塬隧洞、白鹿塬隧洞處于高地下水且地質(zhì)條件復(fù)雜的黃土、泥巖環(huán)境中，在國內(nèi)外非常罕見，缺乏工程經(jīng)驗(yàn)、相關(guān)理論和設(shè)計方法，采用盾構(gòu)法對這類高地下水復(fù)雜地層隧洞施工的安全性有待進(jìn)一步研究。

為了明確引漢濟(jì)渭二期工程輸水隧洞的安全性，必須在施工之前對引漢濟(jì)渭二期工程輸水隧洞的安全性進(jìn)行深入研究，而數(shù)值模擬方法是進(jìn)行輸水隧洞安全性預(yù)測的重要手段。隨著軟件、硬件及算法的發(fā)展，數(shù)值模擬方法在水利工程領(lǐng)域得到了廣泛認(rèn)可[5-9]。輸水隧洞的開挖過程數(shù)值模擬較為復(fù)雜，采用一般的方法進(jìn)行數(shù)值模擬存在工作效率低、出錯率高的缺點(diǎn)。為解決這一問題，本文將編程語言應(yīng)用于輸水隧洞開挖過程的數(shù)值模擬中，以期提出一種高效、可靠、通用的數(shù)值模擬方法，從而提高輸水隧洞數(shù)值模擬的效率。

1 Python在Abaqus二次開發(fā)中的應(yīng)用

Abaqus具有強(qiáng)大的非線性數(shù)值分析功能，含有8類共562種單元，足以對水利工程中的問題進(jìn)行求解[10]。對于輸水隧洞而言，需要能夠進(jìn)行流固耦合的數(shù)值模擬，而其中提供了能夠考慮滲流與應(yīng)力耦合的單元C3D8RP。因此，學(xué)者們[11-14]利用Abaqus對隧洞開挖進(jìn)行了相關(guān)研究，證明了其在隧洞工程中的可靠性。

Python語言目前應(yīng)用于各個領(lǐng)域，且Abaqus的多個模塊基于Python語言編制而成，因此，Python語言可以很好地應(yīng)用于Abaqus的各個過程。學(xué)者們已經(jīng)運(yùn)用Python語言對Abaqus進(jìn)行了部分研究，如郭玲等[15]將Python語言應(yīng)用于薄壁管彎曲起皺的研究中；夏敏等[16]利用Python語言對Abaqus計算平臺進(jìn)行二次開發(fā)，建立了混凝土框架火災(zāi)受損分析系統(tǒng)。這些研究都表明了將Python語言與工程結(jié)構(gòu)領(lǐng)域的數(shù)值模擬軟件相結(jié)合能夠極大提高數(shù)值模擬的效率。

鑒于輸水隧洞開挖過程的復(fù)雜性，并基于Abaqus計算平臺解決非線性問題的能力及Python語言在Abaqus中應(yīng)用的優(yōu)勢。本文將基于Abaqus計算平臺，利用Python語言提出輸水隧洞開挖過程的參數(shù)化分析方法，以提高輸水隧洞開挖過程數(shù)值模擬的效率，為解決引漢濟(jì)渭二期工程輸水隧洞安全性面臨的不確定性問題做出貢獻(xiàn)。

2 Abaqus中隧洞開挖建模過程

隧洞開挖的數(shù)值模擬基本過程包括前處理、計算及后處理，其中操作復(fù)雜且占用時間較多的在于前處理過程。前處理過程基本包括以下步驟：

(1) 建立零件。首先根據(jù)隧洞尺寸繪制橫截面的草圖并對其進(jìn)行擠出，即可建立隧洞開挖模型的整體幾何模型。由于隧洞的開挖幾何模型、開挖過程中應(yīng)變、應(yīng)力及變形等反應(yīng)具有對稱性，因此，本文采用1/2的模型進(jìn)行開挖過程的有限元分析，采用簡化模型可以提高分析效率。本文在繪制草圖的同時對橫截面進(jìn)行了切割，分別按照管片外徑、灌漿圈外徑、距洞較近的圓、圓心與方形角點(diǎn)連線對橫截面草圖進(jìn)行切割，如圖2所示。切割后可保證在內(nèi)部存在圓形隧洞的長方體模型在網(wǎng)格劃分時能夠得到質(zhì)量較好的網(wǎng)格。

圖2 隧道草圖及切割情況

(2) 定義材料并指定截面。為幾何模型指定材料和截面屬性。定義材料時需要指定巖土體的密度、彈性模量、孔隙度、本構(gòu)模型等；同時需要定義襯砌的彈性模量、泊松比、密度等參數(shù)、盾構(gòu)的彈性模量、泊松比、密度等參數(shù)、注漿料彈性模量、泊松比、密度等參數(shù)。在定義好材料屬性之后需定義截面屬性，并將截面屬性指定給各部件。

(3) 裝配零件。將零件組裝成分析用的模型。

(4) 設(shè)置分析步。根據(jù)分析目標(biāo)設(shè)置不同的分析步驟和分析方法，首先，需要定義第一個分析步，即地應(yīng)力平衡分析步；其次，需要定義開挖過程的分析步，對于開挖過程需要利用生死單元技術(shù)設(shè)置用于單元移除與單元重生的幾十個甚至上百個分析步，每一個分析步向前掘進(jìn)1管片的距離為縱向長度。在每一步選擇單元的過程中，不僅要求精準(zhǔn)選擇每一步要移除單元，而且要求同時精準(zhǔn)選擇要激活的單元，此過程中若出現(xiàn)一個單元的選擇錯誤將導(dǎo)致整個分析過程的失敗。

(5) 設(shè)置接觸。單元移除和單元激活采用接觸行為中的model change屬性定義，每一個分析步向前掘進(jìn)一個管片的長度，其過程需要設(shè)置兩個接觸，一個接觸行為是移除，另一個接觸行為是激活。第一步的地應(yīng)力平衡分析步中需要設(shè)置移除隧洞巖土體以外的其他所有單元，即需要利用GUI在關(guān)鍵字中添加*model change, remove字符[17]，并列出需要移除的單元集合，Abaqus在第一步計算地應(yīng)力平衡式就會不計入這些被移除的單元；其次，對于后續(xù)開挖分析步的每一步都需要設(shè)置兩個接觸，即移除該步所需要開挖的土體單元和管片單元，隨后添加該步所對應(yīng)的襯砌單元與注漿料單元。對于某一步來說，移除與激活接觸行為的設(shè)置看似簡單，然而當(dāng)開挖過程較多時，這兩個移除與激活接觸行為的設(shè)置過程就變得非常復(fù)雜。其原因是在移除與激活時需要編寫所要處理的單元集合名稱，這就需要在設(shè)置移除與激活接觸行為之前先定義需要被處理的集合，這些集合的數(shù)量可能上百個。集合的定義過程極易出錯，移除與激活接觸行為字符的編寫過程也極易出錯，且要在不同的分析步中移除不同的接觸，同時激活不同的接觸。可見，這一步的過程非常復(fù)雜，且容易出錯。

(6) 施加荷載及邊界條件。首先，需要在初始條件中設(shè)置地應(yīng)力及孔隙水壓力；其次，在第一步中需要施加重力荷載和上部覆土壓力，在隨后的分析步中保持不變；最后，需要根據(jù)實(shí)際情況設(shè)置模型的邊界條件。對于法線方向?yàn)閄方向的兩個面限制X方向的位移，對于法線方向?yàn)閅方向的兩個面限制Y方向的位移，對于法線方向?yàn)閆方向的兩個面限制Z方向的位移。

(7) 設(shè)置并劃分網(wǎng)格。根據(jù)分析目標(biāo)設(shè)置網(wǎng)格劃分規(guī)則，網(wǎng)格大小、單元類型等，并對網(wǎng)格進(jìn)行劃分。對開挖區(qū)域周圍較近的土體采用較細(xì)的網(wǎng)格，在操作時需切割出細(xì)化網(wǎng)格的區(qū)域，這一步可以在(1)中建立部件的草圖中完成，也可在裝配以后完成。隨后需對不同區(qū)域的網(wǎng)格設(shè)置尺寸，本文中對開挖方向設(shè)置的網(wǎng)格大小均為管片的寬度，目的是與實(shí)際施工過程相對應(yīng)，即開挖一步正好完成該步管片的施工。在模型的橫向，距洞較近的區(qū)域設(shè)置較密的網(wǎng)格，在距凍較遠(yuǎn)的區(qū)域設(shè)置相對較大的網(wǎng)格，這樣即可保證計算精度，亦可保證計算速度。圖3為劃分網(wǎng)格后的盾構(gòu)隧洞有限元模型。

圖3 盾構(gòu)隧洞有限元模型

(8) 檢查模型并提交計算。利用數(shù)據(jù)檢查功能進(jìn)行模型檢查，檢查無誤后提交計算。

3 基于Python語言的隧洞開挖過程參數(shù)化建模方法

3.1 手動建模存在的問題

對于隧洞開挖過程的數(shù)值模擬，如果模型幾何尺寸發(fā)生變化，則上述過程中(1)需要重新修改尺寸，且將導(dǎo)致(5)中定義的數(shù)百個集合都需要重新被選擇和定義，同時也導(dǎo)致(7)中的網(wǎng)格發(fā)生變化和重新劃分。其中，單獨(dú)就(5)而言，重新定義上百個集合、對不同集合在不同的分析步設(shè)置為移除或激活的工作量是非常大的，再加上其他步驟的工作量，這是非常占用時間的。且由于接觸的定義極為復(fù)雜，所以出錯率非常高。以上過程還只是說要修改幾何尺寸的情況，并沒有考慮其他參數(shù)的修改。

針對上述復(fù)雜問題，非常有必要提出一種可靠且通用的數(shù)值模型，從而解決出錯率高、分析效率低的問題。為此，本文將Python語言應(yīng)用于Abaqus計算平臺中，以提出一種高效、可靠的隧洞開挖過程數(shù)值模擬方法。

3.2 基于Python語言的隧洞開挖建模過程

第2節(jié)所述GUI界面操作的各個步驟，均可通過Python語言代替，Python語言建模的所有操作可在幾秒內(nèi)完成且錯誤率為0，可解決界面操作的過程復(fù)雜且出錯率高。下文將詳細(xì)解釋第2節(jié)每一步通過Python編程實(shí)現(xiàn)的方法，并給出部分關(guān)鍵的代碼。

(1) 建立零件

(2) 定義材料并指定截面

定義巖土體密度：Sd_MX_Tu.Density(table=((1590.0, ), ))

彈模和泊松比：Sd_MX_Tu.Elastic(table=((1.5e8, 0.125), ))

摩爾-庫侖模型：Sd_MX_Tu.MohrCoulombPlasticity(table=((18.0, 0.0), ))

Sd_MX_Tu.mohrCoulombPlasticity.MohrCoulombHardening(table=((300000.0,0.0),))

Sd_MX_Tu.mohrCoulombPlasticity.TensionCutOff(dependencies=0,table=((0.0, 0.0), ), temperatureDependency=OFF)

設(shè)置滲透系數(shù)：Sd_MX_Tu.Permeability(inertialDragCoefficient=0.142887, specificWeight=10000.0, table=((1.1e-12, 0.67), ))

體積模量：Sd_MX_Tu.PorousBulkModuli(table=((0.0, 2e9), ))

襯砌、漿料及盾構(gòu)機(jī)均采用彈性材料，設(shè)置方法與巖土體相同。

(3) 裝配零件

(4) 設(shè)置Step

(5) 設(shè)置接觸

首先，定義集合： 編制單元選擇代碼和循環(huán)定義集合代碼將每步需要操作的單元定義為名稱規(guī)則的集合，Sd_MX.rootAssembly.Set(elements=所選擇的單元, name=集合名稱)；

其次，通過修改關(guān)鍵詞定義每一步要移除的單元：Model change, remove；

同時，通過修改關(guān)鍵詞定義每一步要激活的單元集合: Model change, add；

(6) 施加荷載及邊界條件

(7) 設(shè)置并劃分網(wǎng)格

定義邊種子：Sd_MX.rootAssembly.seedEdgeBySize(constraint=FINER, deviationFactor=0.1, edges=邊, minSizeFactor=0.1, size=網(wǎng)格尺寸)

劃分網(wǎng)格：Sd_MX.rootAssembly.generateMesh(regions=(需要劃分網(wǎng)格的部件名稱, ))

(8) 檢查模型并提交計算

將上述過程編制成Python文件，對于不同的模型只需要修改不同參數(shù)的數(shù)值大小，然后在Abaqus平臺中運(yùn)行該文件就可以在數(shù)秒內(nèi)完成復(fù)雜的隧洞開挖支護(hù)的建模工作。更重要的是采用Python進(jìn)行建模具有零錯誤率的特點(diǎn)，且應(yīng)用時只需要修改某些參數(shù)的數(shù)值就可以迅速完成數(shù)值模型的建立及計算過程。

通過對比可知，手動建模操作過程繁瑣、工作效率低、準(zhǔn)確率差，而基于Python語言的建模方法具有可重復(fù)性高、工作效率高、零出錯率的特點(diǎn)。因此，基于Python語言的建模方法在復(fù)雜的隧洞開挖問題中的應(yīng)用前景廣闊。

4 參數(shù)化模擬方法的應(yīng)用算例

本文以白鹿原隧洞為例，采用上述方法分析洞周圍橫向網(wǎng)格尺寸變化對盾構(gòu)隧洞開挖與支護(hù)過程的影響，以驗(yàn)證本文所述方法的可靠性。首先，利用參數(shù)化模擬方法建立不同橫向網(wǎng)格尺寸的隧洞開挖有限元模型；其次，利用各模型進(jìn)行開挖和支護(hù)的全過程模擬，對比結(jié)果的差異。

4.1 不同橫向網(wǎng)格尺寸的有限元模型

本文所建立的不同橫向網(wǎng)格尺寸的有限元模型如圖1所示，所有模型中開挖方向網(wǎng)格尺寸與管片寬度相同，洞徑方向網(wǎng)格分為兩個區(qū)域，即距洞較遠(yuǎn)處的粗網(wǎng)格區(qū)域和距洞較近的細(xì)化網(wǎng)格區(qū)域。圖4(a)中的黃色區(qū)域即為細(xì)化網(wǎng)格區(qū)域，圖4(b)—圖4(f)模型中細(xì)化區(qū)域的網(wǎng)格尺寸分別為20 cm、50 cm、80 cm、110 cm、140 cm，所分析的五個模型只有細(xì)化區(qū)域網(wǎng)格尺寸變化，其他參數(shù)均相同。

圖4 細(xì)化區(qū)域不同網(wǎng)格尺寸的有限元模型

4.2 有限元分析結(jié)果

盾構(gòu)隧洞開挖與支護(hù)有限元分析完成后取出位于距離模型前端30 m處的管片，對比不同細(xì)化網(wǎng)格尺寸下該處襯砌結(jié)構(gòu)的應(yīng)力分布情況，結(jié)果如圖5所示。圖5(a)為所分析的管片位置的示意圖，圖5(b)—圖5(f)為細(xì)化區(qū)域不同網(wǎng)格尺寸時管片應(yīng)力的分布情況。

圖5 不同橫向網(wǎng)格尺寸時開挖并支護(hù)后距端面30 m處管片應(yīng)力

從以上分析結(jié)果可以看出，當(dāng)洞周圍細(xì)化區(qū)域網(wǎng)格尺寸分別為20 cm、50 cm、80 cm、110 cm、140 cm時，開挖支護(hù)后的管片結(jié)構(gòu)應(yīng)力分布沒有變化。然而，圍巖壓應(yīng)力隨著網(wǎng)格尺寸變小呈減小趨勢，且變化幅度不大。當(dāng)細(xì)化區(qū)域網(wǎng)格尺寸為20 cm、50 cm、80 cm、110 cm、140 cm時，圍巖最大壓應(yīng)力分別為2.563 MPa、2.580 MPa、2.591 MPa、2.606 MPa、2.617 MPa。可見，洞周圍細(xì)化區(qū)域網(wǎng)格尺寸對開挖支護(hù)后管片與圍巖受力幾乎沒有影響。上述算例也證明了本文所提參數(shù)化分析方法在進(jìn)行參數(shù)敏感性分析方面的優(yōu)勢，該方法大幅減少了復(fù)雜的重復(fù)性工作，提高了工作效率和準(zhǔn)確性。

5 結(jié) 論

為解決引漢濟(jì)渭二期工程輸水隧洞開挖過程、支護(hù)過程的隧洞圍巖及支護(hù)結(jié)構(gòu)受力性能、安全性分析中所面臨的復(fù)雜問題，本文將Python編程語言應(yīng)用于隧洞開挖過程的數(shù)值模擬中。本文基于Python語言編制了隧洞開挖過程的代碼，可代替界面中手動地操作以完成復(fù)雜和重復(fù)性的建模過程。算例表明本文方法能夠極大提高隧洞安全性評價的工作效率，可為加快引漢濟(jì)渭二期工程的進(jìn)度做出貢獻(xiàn)，能夠用于類似隧洞工程的數(shù)值模擬，具有實(shí)用價值和參考意義。