王 旭，賴國偉，王志鵬

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

0 引言

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

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

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

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

1.1 “另類單元生死”

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

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

1.2 “另類單元生死”的實現(xiàn)方法

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

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

*IMPORT NSET：把當(dāng)前施工步的結(jié)點組件傳遞到下一個施工步；

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

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

1.3 顯式算法

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

1)哪一種圖形密碼最安全；2)圖形密碼與傳統(tǒng)的六位數(shù)字密碼，哪一個更安全；3)改變構(gòu)造圖案的規(guī)則，密碼的安全程度會有哪些改變。

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

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。有限元網(wǎng)格劃分見圖1。其中，壩體及巖體單元以六面體單元為主，四面體單元為輔。整個計算模型有48 609個單元，其中壩體有25 136個單元。模型整體坐標(biāo)系X正向指向右岸，Y正向指向順?biāo)飨?，Z正向為鉛直向上。

圖1 小灣拱壩計算模型

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

2.2 基本資料

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

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

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

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

2.3 材料本構(gòu)關(guān)系

材料本構(gòu)關(guān)系采用理想彈塑性模型[7，8]，屈服準(zhǔn)則采用常用的Drucker-Prager準(zhǔn)則，即

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

2.4 荷載及工況

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

表1 施工過程模擬

3 計算結(jié)果及分析

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

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

表2 壩體最大位移和最大應(yīng)力及其位置

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

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

4 結(jié)論

圖3 第一、第三主應(yīng)力等值線 （單位：MPa）

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

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

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

.

［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］ 羅州.金屬環(huán)件軋制過程剛粘塑性動力顯式算法有限元模擬［D］.武漢：武漢理工大學(xué)，2004.

［4］ 王素裹，韓小雷，季靜.ABAQUS顯式分析方法在鋼筋混凝土結(jié)構(gòu)中的應(yīng)用［J］.科學(xué)技術(shù)與工程， 2009， 9（16）:4688-4692.

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

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

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

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

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