“另類單元生死”在大壩等結構施工仿真中的應用

王 旭，賴國偉，王志鵬

（武漢大學水資源與水電工程科學國家重點實驗室，湖北 武漢 430072）

0 引言

采用隱式有限元方法求解復雜工程問題已經十分普及，但只適用于小位移、小轉動、小變形的界限內。隱式有限元方法會經常因局部材料破壞導致迭代不收斂，所以一旦需要模擬施工過程并且還要進行高度非線性分析時，隱式算法會產生嚴重的收斂困難問題。顯式有限元是求解高度非線性問題的一個非常有效的工具，在工程界和學術界已經取得高度的認可。Hooputra， H.等[1]利用顯式準靜態方法分析了薄壁鋁材逐步失效破壞，并與實驗數據進行了對比驗證；Linde P等[2]模擬纖維金屬成型得到了理想的結果；羅州[3]采用顯式有限元分析程序ABAQUS/Explicit提供的動力顯式算法，并加入與實際工程更相似的約束條件，從而得到了更為準確的結果。在混凝土工程方面，王素裹等[4]為了考察ABAQUS/Explicit顯式分析方法在鋼筋混凝土結構上應用的適用性和準確性，用顯式與隱式兩種分析方法對迭代算法和求解時間進行了對比，計算結果表明，顯式分析方法能夠在解決材料失效和破壞導致的收斂問題的基礎上，較為準確地模擬鋼筋混凝土結構響應。

ABAQUS/Explicit顯式有限元方法最顯著的特點是不需要在隱式方法中用到的整體剛度矩陣，顯式地前推模型的狀態，不需要迭代和收斂準則；而且在相同的模擬條件下，顯式方法比隱式方法需要更少的磁盤空間和內存。但是ABAQUS/Explicit顯式有限元沒有類似于隱式的 “單元生死”，所以無法模擬施工過程。如果不模擬施工過程，仿真模擬得到的結果會與現實的結果差別很大。

為了使用ABAQUS/Explicit顯式有限元強大的求解器來解決工程項目的諸多高度非線性問題并模擬施工過程，本文提出通過數據傳遞的方式來模擬單元生死，并將此方式定義為 “另類單元生死”法。基于該方法，模擬了小灣拱壩的施工過程，并將計算結果與Ansys隱式計算結果進行對比。

1 “另類單元生死”和顯式算法

1.1 “另類單元生死”

本文自定義的 “另類單元生死”原理是[5]，通過ABAQUS/Explicit模塊中有限元模型和應力變形場等數據傳遞的方法，亦即將結構前一荷載步的有限元模型和計算應力變形場等數據作為下一荷載步計算的初始輸入數據，來解決ABAQUS/Explicit顯式有限元計算的單元生死問題，達到模擬結構施工過程的目的。

“另類單元生死”法與隱式有限元的 “單元生死”法不同。隱式有限元的 “單元生死”是根據施工步的順序而改變相應單元剛度并且設定與單元相聯系的質量、阻尼、比熱及其他類似的效果的值。而 “另類單元生死”則是保留當前施工步僅有的有限元模型，刪除后面所有施工步出現的單元和結點，然后對當前施工步進行求解；之后，將當前施工步的有限元模型和計算結果傳遞到下一個施工步，作為下一個施工步定義初始狀態場，下一個施工步即在當前施工步的有限元模型和計算結果基礎上，增加相應的單元和結點，接著進行求解。

1.2 “另類單元生死”的實現方法

基于ABAQUS/Explicit顯式有限元，用 “另類單元生死”模擬施工步，把當前施工步的有限元模型和計算結果傳遞到下一個施工步的ABAQUS命令流文件，亦即inp文件格式如下所示：

*IMPORT，STEP=N，UPDATE=no：把當前施工步的計算結果傳遞到下一個施工步，其中，N表示當前施工步為第幾個施工步；

*IMPORT NSET：把當前施工步的結點組件傳遞到下一個施工步；

*IMPORT ELSET：把當前施工步的單元組件傳遞到下一個施工步。

采用 “另類單元生死”法需注意的是，ABAQUS計算所產生的后綴名為abq、stt、pac、prt和odb的文件必須保存在同一工作目錄下，才能保證計算順利進行。

1.3 顯式算法

ABAQUS/Explicit應用中心差分方法對運動方程進行顯式的時間積分[6]，由一個增量步的動力學條件計算下一個增量步的動力學條件。計算公式如下：

1)哪一種圖形密碼最安全；2)圖形密碼與傳統的六位數字密碼，哪一個更安全；3)改變構造圖案的規則，密碼的安全程度會有哪些改變。

該算法不需要迭代和收斂準則，為了得到精確的結果，時間增量步必須相當小，這樣在增量步中加速度幾乎為常數。由于時間增量步必須很小，一個典型的分析需要成千上萬個增量步，但在求解過程中，不必同時求解聯立方程組，故每一個增量步的計算成本并不高。

2 工程實例

2.1 小灣拱壩及計算模型

小灣大壩為混凝土雙曲拱壩，壩頂高程1 245 m，壩基面高程953 m，最大壩高292 m，壩頂寬12 m，最大壩底寬72.912 m，壩頂弧長901.771 m，大壩混凝土量842萬m3，弧高比3.056，弦高比2.74，厚高比0.25。拱壩承受的總水推力為1.66×108kN。有限元網格劃分見圖1。其中，壩體及巖體單元以六面體單元為主，四面體單元為輔。整個計算模型有48 609個單元，其中壩體有25 136個單元。模型整體坐標系X正向指向右岸，Y正向指向順水流向，Z正向為鉛直向上。

圖1 小灣拱壩計算模型

為檢驗 “另類單元生死”法的正確性，ABAQUS/Explicit顯式有限元和Ansys隱式有限元采用同一有限元模型。

2.2 基本資料

（1）水位。上游正常蓄水位1 240.0 m；下游水位1 004.0 m。

（2）泥沙。淤沙高程1 097 m，淤沙容重9.5 kN/m3， 內摩擦角 24°。

（3）依據混凝土設計參數為：容重γ=24 kN/m3，泊松比 γc=0．2， 線膨脹系數 α=8．26×10-6/℃， 抗剪強度 f′=1．4， c′=1．6 MPa。 不同混凝土分區 （A 區、 A0區， B 區， C 區） 180 d彈性模量分別為 23．1、 22．4、21．7 GPa，已按工程經驗考慮了混凝土的徐變影響。

（4）壩基采用均質模型，彈性模量21 GPa，泊松比 0．26， 基礎巖體抗剪強度 f′=1．4， c′=1．8 MPa。

2.3 材料本構關系

材料本構關系采用理想彈塑性模型[7，8]，屈服準則采用常用的Drucker-Prager準則，即

式中，I1為應力張量的第一不變量；J2為應力偏張量的第二不變量；α、k均為材料常數，與材料粘結力c、內摩擦角φ的關系取決于Drucker-Prager圓錐面與摩爾庫侖六棱錐面之間的相互關系。本文中α、k采用外接圓錐公式計算

2.4 荷載及工況

計算中考慮的荷載有拱壩自重、上下游壩面水壓力、上游泥沙壓力、運行期溫度荷載、水庫庫盆水壓力。自重和上游壩面水壓力按預備工況分期施加 （見表1）；下游壩面水壓力、上游泥沙壓力、運行期溫度荷載、水庫庫盆水壓力均一次性施加，考慮由面力代替模擬滲流場的體積力，乘以折減系數，折減系數取 0．8。

表1 施工過程模擬

3 計算結果及分析

考慮到Ansys隱式有限元計算[9]通過隱式的 “單元生死”法能夠模擬施工過程，為驗證 “另類單元生死”法的正確性，本文將分別采用ABAQUS/Explicit顯式有限元與Ansys隱式有限元進行對比計算分析，以拉應力為正。順河向位移等值線見圖2，第一、第三主應力等值線見圖3。壩體最大位移和最大應力及其位置見表2。

圖2 順河向位移等值線 （單位：m）

表2 壩體最大位移和最大應力及其位置

從圖2、3和表2可以看出：①順河向位移分布大致相同，而且都呈左右對稱分布，用兩種軟件計算的順河向最大位移最大值差別不大，并且出現的位置基本一致；②第一主應力在壩頂左右有少許不同，其中用ABAQUS/Explicit計算的數值稍大，但是整體的第一主應力分布基本相似，且都呈左右對稱分布，用兩種軟件計算的第一主應力最大值差別不大，并且出現的位置基本一致；③第三主應力在壩體1/4壩高處有少許差別，其中用ABAQUS/Explicit計算的數值稍偏大。雖然兩者局部有所不同，但是從整體上來看，第三主應力分布大致類似，且都呈左右對稱，除此之外，用兩種軟件計算的第三主應力最大值差別不大，并且出現的位置基本一致。

綜上所述，對比ABAQUS/Explicit顯式有限元與Ansys隱式有限元的位移、應力及其分布情況可以看出，基于 “另類單元生死”法的顯式有限元計算結果與傳統隱式有限元區別很小，表明ABAQUS/Explicit顯式有限元采用 “另類單元生死”法，通過數據傳遞命令能夠傳遞結構施工荷載步的位移、應力等信息，計算結果正確。

4 結論

圖3 第一、第三主應力等值線 （單位：MPa）

（1）“另類單元生死” 法的原理是， 通過ABAQUS/Explicit模塊中有限元模型和應力變形場等數據傳遞的方法，亦即將結構前一荷載步的有限元模型和計算應力變形場等數據作為下一荷載步計算的初始輸入數據，來解決ABAQUS/Explicit顯式有限元計算的單元生死問題，達到模擬結構施工過程的目的。

（2）在塑性破壞的情況下，ABAQUS/Explicit顯式有限元計算結果與Ansys隱式計算結果差別很小，表明ABAQUS/Explicit顯式有限元通過數據傳遞方式能夠模擬施工過程， “另類單元生死”法可以達到隱式的 “單元生死”要求。

（3）借助 “另類單元生死”法，廣大科研人員可以用ABAQUS/Explicit顯式有限元強大的非線性求解器對大壩等結構進行高度非線性仿真分析。

.

［1］ Hooputra H,Gese H,Dell H,etal.A Comprehensive Failure Model for Crashworthiness Simulation of Aluminium Extrusions［J］.International Journal of Crashworthiness,2004（9）:449-463.

［2］ Linde P,Pleitner J.Modelling and Simulation of Fiber Metal Laminates［A］.ABAQUS Users'Conference,2004.

［3］ 羅州.金屬環件軋制過程剛粘塑性動力顯式算法有限元模擬［D］.武漢：武漢理工大學，2004.

［4］ 王素裹，韓小雷，季靜.ABAQUS顯式分析方法在鋼筋混凝土結構中的應用［J］.科學技術與工程， 2009， 9（16）:4688-4692.

［5］ 曹金鳳，石亦平.ABAQUS有限元分析常見問題解答［M］.北京：機械工業出版社，2009.

［6］ 莊茁，由小川，等.基于ABAQUS的有限元分析和應用［M］.北京：清華大學出版社，2009.

［7］ 王金昌，陳葉開.ABAQUS在土木工程中的應用［M］.杭州：浙江大學出版社，2006.

［8］ 陳明祥.彈塑性力學［M］.北京：科學出版社，2007.

［9］ 王國強.實用工程數值模擬技術及其在ANSYS上的實踐［M］.西安：西北工業大學出版社，1999.