摘要： FEM是目前許多行業數值計算的主要方法。本文結合二灘地下廠房洞室群的圍巖穩定分析，主要討論了進行大型有限元建模及計算時的相關問題，并對洞室群的圍巖穩定進行了力學分析，給出了圍巖整體是穩定的結論。
關鍵詞： FEM 圍巖穩定 相關問題
許多工程分析問題，如固體力學中位移場和應力場分析、熱學中的溫度場分析、流體力學中的流場分析等，都可借助于計算機數值模擬技術來求出它們的數值解，其中有限單元法（Finite Elements Method，以下簡稱 FEM）是工程中最常用的數值模擬技術之一。本文結合二灘水電站地下廠房洞室群的圍巖穩定分析，對使用有限元軟件進行大規模數值模擬計算時的相關問題進行探討，給出了建議。
1 建模與有限元網格劃分 1.1 工程概況二灘地下廠房洞室群設置在左岸下游的山體內，上覆巖體厚度為250~350m，水平埋深300m，洞軸線方向NW6°。洞室群所在區域主要為正長巖和蝕變玄武巖，圍巖巖體有80％質量達到A,B級。地下廠房洞室群由三大洞室和母線洞以及尾水洞組成。其中主廠房：長×寬×高＝280.29m×25.5m×65.78m；主變室：214.9m×18.3m×25.0m；尾調室：203.0m×19.8m×69.8m。主廠房主要采用噴錨支護，在主副廠房和尾調室高邊墻采用了一定數量的預應力錨索。地下洞室群水平區域布置如圖1。
1.2 建模 一般的有限元計算，70%~80%的時間花在建立模型上面，因此選擇合適的建立模型的方法是很重要的，否則將有可能事半功倍。
二灘洞室群位于地下，結構復雜，操作起來比較困難。為此把大模型分解為許多很小的模型，劃分單元以后再把它們拼裝成一個大的模型；并且不能把它們合成一個整體，否則將無法劃分單元。具體的作法是：先建立洞室群的模型，然后在長、寬、高方向上分別向外延伸，形成一個長×寬×高=404.2m×283.7m×185.3m的巨大巖體，用該巖體挖去洞室群，就是所需要的模型。由于有限元軟件的自適應網格劃分功能只能劃分規整的模型，如六面體或規則的棱柱體（如三角棱柱體），而地下廠房洞室群則是縱橫交錯，極其不規則的，建模難度可想而知。經過不斷的修改和試驗，最終建立的模型共有3647個體組成。
盡管如此，本文模型與洞室的布置相比，還是進行了一些簡化，主要有：（1）出于單元形狀不致過分奇異以及單元數量方面的考慮，忽略了尾調室和主變室之間的夾角，將它們的軸線設為平行。（2）由于僅考慮三大洞室及其之間的圍巖變形與破壞，忽略了進水洞和兩個尾水洞。（3）忽略了洞徑較小并且僅有局部影響的施工支洞、交通洞、電纜斜井和排水廊道等。
建模過程中，需要注意的有以下幾點：
1)對模型內部操作（如空洞）時，可以先生成空洞部分，再和外部模型體組合成所需要的模型。
2)對于任何體的某一個面，都應當是凸多邊形（即圖形必須在任何一邊或其延長線的一側），且其邊數不應超過4，否則將會出現拓撲錯誤，無法劃分單元。
3)合并重復的實體圖元，如相同位置的兩個關鍵點或線，以節省存儲空間。
4)不能出現扭曲的面，但允許是弧面。
1.3 網格劃分
此工程模型中間是地下廠房的洞室群，為此，先給洞室群部分劃分單元，然后給廠房鄰近的模型體劃分單元，這里要注意的一點是要保證它們的接觸部分的節點能夠吻合。最初采用的是四面體單元，之所以使用該種單元，主要是它的自適應劃分能力較強，比較適用于這種復雜的模型。但是當劃分網格完成后，單元數量過多（達10萬多個），并且出現一些非常小的尖角（小于15度），對應的矩陣難以收斂，故更改為6面體8結點單元。又做了一些簡化，并調整單元長度，最終劃分單元34224個，結點數量為37632個，計算規模仍然是十分巨大的（劃分單元后整體模型如圖2，洞室群部分模型如圖3）。在劃分單元時應注意的問題：
1)必須要保證所劃分的單元和鄰近體的單元能夠吻合，即使不能完全吻合，也要使誤差控制在一定的范圍內。
2)劃分單元完成后，應合并整個模型的所有節點。否則，由于模型是各個獨立的體，加荷載計算后，在體與體的結合處將會出現開裂。
3)壓縮關鍵點、線、面、單元、節點等編號，以提高運算效率，節省機時。
4)在應力集中的地方單元網格密度適當加密；反之，則可以加大劃分單元的邊長，以達到不致太多的影響計算精度，又提高了計算速度。
地下洞室的開挖過程是地應力逐步釋放的過程。在本工程中，主要采用應力函數法和有限元分析法相結合[1]。具體做法如下：
第1步：在計算模型邊界的某些點處試探施加給定位移得到參考計算結果；然后采用應力函數法在有限元模型的范圍內進行回歸。回歸計算中，顧及上述參考結果以及實測測點應力的最大、最小和平均值，在計算模型的邊界上設定若干點的應力分量的值，并使這些邊界條件同時參與回歸計算；進而得到在建模范圍內地應力的第1步近似結果。
第2步：利用第1步的結果計算出在建模邊界上所有單元結點的應力，采用類似邊界荷載調整法的做法把它們作為遠場外荷載施加到有限元模型的邊界上。此時由于邊界荷載采用的是滿足彈性力學中變形協調方程和平衡方程的應力場所求得的，邊界上的正應力和剪應力可以共同保證構成平衡力系，并且解決了求剪應力場所遇到的困難。利用上述遠場邊界荷載，借助有限元程序計算建模范圍內的地應力，然后將此結果在實測地應力的測點上與實測地應力進行比較，并且考慮在洞室群范圍內地應力的分布情況結合工程概念考察回歸計算結果的合理性。根據所得到的趨勢，調整在第1步計算中施加在邊界上的有限各點的力學邊界條件，再次回歸計算內部各點的地應力，并由此獲得全部邊界上由遠場位移所引起的所有應力分量。
重復以上兩步的分析計算，直到得到滿意的應力場為止。
在本次計算中，采用雙參數準則中的Drucker-Prager準則來判斷復雜應力狀態下巖體的破壞。該準則是Mises準則的推廣，因此又被稱為廣義Mises準則[2]。該準則認為，處于彈性區的應力滿足式(1)時，則認為巖體處于破壞的臨界狀態。
（1）
其中：
（2）
（3）
C,φ分別為凝聚力和摩擦角；Rc為單軸抗壓強度。
當巖體材料處于彈性階段時，其本構關系為：
（4）
其中：為彈性矩陣，分別為應力矩陣和應變矩陣。
當材料進入塑性以后采用增量理論的本構模型。即認為總應變增量可以表示為彈性應變增量與塑性應變增量之和： 。
在施工期間，二灘水電站地下廠房的1號尾水洞和2號尾水洞之間曾發生過巖爆及塌方，故取2號機組所在的橫斷面作為破損區的示意圖[3]。由圖4可以看出：
1）橫斷面上破壞主要集中在主廠房下游邊墻1號和2號兩個母線洞之間，對應的兩個尾水洞之間，以及主變室上游邊墻對應的兩個母線洞之間的部位，這與該處的地應力較大及結構復雜有關。
2）母線洞底板破壞的位置集中在母線洞中部到主廠房一側。
圍巖的垂直位移分布規律十分明顯：頂拱向下，底板向上。這是由于開挖引起地應力釋放所造成的。圍巖總位移的分布規律則基本上是按1號尾水洞到6號尾水洞的方向依次遞減（見圖5）。各主要洞室的位移分布如下：
1）主廠房上游邊墻位移的最大值出現在邊墻的中部，下游邊墻水平位移的最大值出現邊墻上部和頂拱位置。恰好在此位置上，巖爆所造成的破壞最嚴重。
2）主變室高度較低，約為主廠房高度的2/5，但是其最大位移值卻大于主廠房。
3）尾調室的最大位移發生在1號和2號尾水洞之間。
與破損區圖相對照，可以看出是否破壞以及破壞程度與絕對位移的大小不一定有確定的對應關系。最明顯的是在主廠房的上游拱肩的位移比下游拱肩的位移為小，但是下游拱肩并沒有破壞，而上游拱肩則發生了破壞，這主要是由于下游拱肩的剛體位移較大，形變位移卻不大的原因。
由以上可看出，地下廠房洞室群的圍巖整體是穩定的，破損區是局部的，并且由于該巖體的完整性較好，故可通過適當噴錨處理。圍巖的破壞從1號尾水洞到6號尾水洞依次遞減，故應加強1號和2號尾水洞之間巖體的穩定性觀測，遇到問題及時處理。
圖1二灘地下廠房區域水平布置示意圖
|
|
圖4 2號機組所在面破損區等值線示意圖 圖5主廠房下游邊墻總位移示意圖
3.3 計算時應注意的問題
在計算前，可以編制計算文件，計算時直接調入該文件執行即可。例如要模擬開挖的過程就可以編制相應的步驟，在計算不同的開挖步驟時，只需調入不同的文件執行即可。這樣就減少了人工干預，提高了計算效率。在計算后，若出現模型開裂的情況，則說明在不同的小模型之間尚未完全合并，需調整合并圖元的容差，直至所有小模型合并成一個整體為止。
4 結束語 利用有限元軟件進行數值模擬計算已成為當前許多行業數值計算的主流，本文結合具體工程，討論了在運用此類軟件進行大規模數值計算時的相關問題，主要有以下幾點：
1）建立模型時，整體考慮選擇合適的模型和建模方法，盡量符合工程原型且要便于計算分析。
2）把大模型分解為小模型，做到化繁為簡。
3）計算時，可編制計算步驟文件，在不同的載荷步調用不同的文件即可。
參考文獻：
[1] 張有天,胡惠昌.地應力場的趨勢分析.水利學報.1984,(4).
[2] 徐秉業,劉信聲.應用彈塑性力學.北京:清華大學出版社,1995.513~514.
[3] 二灘水電站地下廠房洞室群模擬開挖過程的三維非線性有限元分析.昆明理工大學電力工程學院,2000.
