陳燈紅 劉云龍 林天成 潘子悅
(1.防災減災湖北省重點實驗室(三峽大學), 湖北 宜昌 443002;2.三峽大學 土木與建筑學院, 湖北 宜昌443002)
我國西南部地區(qū)是強震高發(fā)區(qū),同時我國大型300 m 級高混凝土壩主要集中在西南部地區(qū),如白鶴灘拱壩的設計地震加速度為0.325g,大崗山拱壩甚至達到了0.557 5g.一旦這些高壩大庫發(fā)生潰壩災變,后果不堪設想,因此,地震作用下的高混凝土拱壩抗震安全研究具有重大意義.
在實際工程中,如水工結構抗震領域,大量問題的求解往往伴隨著高度的非線性,這些問題被抽象為大規(guī)模偏微分方程組,而現有的數學理論往往無法求得解析解,數值算法便應運而生.經典有限元法作為一種成熟的偏微分方程的數值求解方法,已經廣泛應用于各類工程問題的求解,各種商用程序中也多采用經典的有限元方法,如商用程序ABAQUS,因其強大的非線性方程組求解能力而被學者和工程師廣泛使用.但隨著擬建的實際工程體量日益增加、工程服役的環(huán)境日益復雜,對工程分析的計算精度、計算效率要求亦相應提高,有必要對現有的工程結構動力分析方法進行進一步的探索.
比例邊界有限元法(the scaled boundary finite element,SBFEM)作為近年來的一種新型偏微分方程組求解算法,最早由Wolf、Song[1]提出,以其具有半解析解、降維求解等優(yōu)勢逐漸成為一種主流的數值計算方法.在無限域的理論研究中,比例邊界有限元法被廣泛用于無限域動力剛度計算[2],這是由于SBFEM 可以將研究域邊界向外擴張,進而準確、便捷地表征無限域.在無限域問題的工程應用中,土-結構相互作用[3]、壩體-庫水相互作用[4]領域有一大批成果,Qu等[4]提出了基于二維SBFEM 無限水域、半空間無限土體和重力壩結合的數值計算方法.許賀等[5]提出一種耦合SBFEM-FEM 方法用于研究壩體-庫水相互作用.Li等[6]將SBFEM 用于研究波在成層無限域的傳播問題,并成功應用于計算壩體地基的動力剛度[7].Xu等[8]研究了可壓縮庫水-壩體的非線性動力相互作用問題.在斷裂和損傷問題中,由于SBFEM 的半解析特性,應力強度因子可以直接獲得[9],Egger等[10]提出了一種提高應力強度因子計算精度的舒爾解耦和超級收斂恢復技術.Jiang等[11]將SBFEM 和四叉樹網格形式結合起來用于研究細觀尺度下的混凝土開裂過程.李建波等[12]應用基于水平集算法的擴展比例邊界有限元法研究了裂紋擴展問題.劉鈞玉等[13]研究了壩體-庫水-地基系統(tǒng)的動態(tài)斷裂問題.Zhang等[14]研究了SBFEM 非局部損傷模型中的應用,在消除網格敏感性的基礎上采用四叉樹網格加密了局部損傷區(qū)域.該方法在三維問題中同樣適用[15].杜成斌等[16]研究了SBFEM 在非局部宏觀與微觀尺度下的損傷模型并進行開裂模擬.Natarajan等[17]將相場模型與SBFEM 結合用于研究脆性材料在動力作用下的斷裂.在接觸非線性問題中,Zhang和Song[18]提出一種將非匹配網格轉化為匹配網格且在局部進行網格細化的技術,增加了將多個獨立計算模型結合用于數值仿真的靈活性.Xing等[19]將該方法拓展到三維問題中,采用點對點接觸模型對匹配網格的接觸面進行離散,降低了求解難度.Zhang等[20]將SBFEM 應用于開裂面的接觸分析,并對靜動力作用下的裂紋開展問題進行了研究[21].Chen等[22]對結構與半空間外域接觸問題進行了研究.
近年來,將SBFEM 與商業(yè)軟件結合逐漸稱為一種主流觀點,Ya 等[23]率先對基于ABAQUS 的SBFEM 進行了研究.孔憲京等[24]自主研發(fā)了大型巖土工程分析軟件,并對高土石壩工程進行了精細化數值分析.但目前與商用程序相結合的研究還有待進一步加強,因此考慮建立一種基于ABAQUS的比例邊界有限元工程結構動力分析方法有一定意義.本文在Ya等[23]的研究基礎上發(fā)展了一種不引入額外自由度的蒙層單元技術,將SBFEM 與ABAQUS中的接觸非線性模型結合,建立了考慮橫縫接觸非線性的壩體-庫水-地基相互作用模型.通過與ABAQUS 計算結果對比,說明該方法為工程結構動力分析提供了一種選擇.
為增強自身的兼容性,商用程序ABAQUS提供了一系列的二次開發(fā)接口,其中子程序接口UEL(user defined element)為用戶實現自己的理論提供了便捷的方法,文章即是采用該接口將比例邊界有限元法嵌入ABAQUS中.
比例邊界有限元法控制方程為:
系數矩陣為:
式中:[D]為彈性矩陣;ρ為密度;[Jb]為Jacobian矩陣行列式,矩陣[Bi]為形函數矩陣,意義同經典有限元理論.
圖1展示了有兩個單元多面體單元組成的模型圖,圖2為圖1的補充單元拓撲信息.該簡單模型的坐標原點為節(jié)點3,兩個單元均為七面體單元,如果僅僅依據普通的單元節(jié)點信息,是無法抽取多面體的所有面節(jié)點信息的,補充了額外的拓撲信息以后,通過索引的方式,所有的面節(jié)點信息可以很方便地得到.注意單元信息中的符號,該下部單元的第1、5和6個單元的符號為負,以第一個面為例,相對于比例中心坐標(0.5,0.5,0.5),節(jié)點1、2、3與比例中心滿足右手法則,因此符號為負號,對于其他符號為正的面,它們的面節(jié)點編號順序相對于比例中心不滿足右手法則.
圖1 多面體單元示例
HHT(Hilber-Hughes-Taylor)時程積分法是ABAQUS內部采用的一種隱式積分方法,該方法是Newmark積分法的拓展之一.在HHT 時程積分法中多質點體系的運動方程表示為:
ABAQUS與UEL子程序之間傳遞數據的過程中需要遵循嚴格的規(guī)則.ABAQUS 傳進子程序的LFLAGS數組包含了關鍵的分析類型信息,表明了當前增量步ABAQUS 主程序的分析狀態(tài),UEL 計算的關鍵量即是根據該數組的不同取值而輸出到主程序的.UEL 子程序返回給主程序的關鍵量中最重要的為AMARTIX 矩陣和RHS矩陣,這二者的計算結果根據LFLAGS的值不同而不同.本文主要涉及到LFLAGS(3)和LFLAGS(1)的取值,AMARTIX矩陣和RHS矩陣相應的具體取值如圖3所示.
圖3 不同分析狀態(tài)對應的關鍵矩陣
其中,關鍵矩陣的不同取值取決于上下標LFLAGS數組的取值,具體的表達式如圖4所示.
圖4 不同分析狀態(tài)的關鍵矩陣取值
作為成熟的商用程序,ABAQUS有著嚴格的程序流程設計,采用了 UEL 子程序接口之后,ABAQUS內部的部分計算流程會被替換,自編程序會與ABAQUS主程序有著大量數據交互.圖5為基于SBFEM 開發(fā)的UEL用戶子程序計算流程圖,圖6為SBFEM_UEL程序中函數的調用關系.
圖5 基于UEL的ABAQUS計算流程圖
圖6 基于SBFEM 的UEL計算流程圖
文章提出了一種不引入額外自由度的蒙皮單元技術實現了耦合SBFEM-FEM.蒙皮技術是用戶自定義單元與ABAQUS豐富的內置功能相結合的一種思路,這種思路的根本目的是利用ABAQUS自帶單元庫的一些功能以簡化用戶子程序接口的一些繁瑣操作,規(guī)避了直接在用戶子程序接口中進行大量繁瑣的編程,通過直接定義蒙皮單元便可以完成復雜荷載的施加.以下介紹蒙皮技術的兩種具體應用方法,分別是全覆蓋蒙皮單元技術和通過直接耦合完成用戶自定義單元與ABAQUS內置單元庫的結合計算.
圖7為通過蒙皮單元給用戶自定義單元施加重力的示意圖,圖中等式左端第一項意為用戶自定義單元,藍色為單元節(jié)點,為示意繪制了單元的線,實際上ABAQUS是不識別的,左端第二項為ABAQUS 內置單元庫中相同幾何拓撲信息的單元,等式右端即為蒙皮單元與用戶自定義單元結合后的示意圖,結合后的單元力學量來自于內置單元和用戶自定義單元,顏色為藍-紅相間.
通過蒙皮技術給用戶自定義單元施加重力的數學原理見式(4).
式中:下標S表示該力學量來自用戶自定義單元,下標F表示該力學量來自經典的有限元法,即ABAQUS的內置單元,下標G 表示重力項,下標ext表示節(jié)點的其余迭代不平衡量.
由式(4)可知,方程右端量已經具有了重力項,這是內置單元的形函數插值的結果.這里需注意的是內置單元與用戶自定義單元的屬性定義的細節(jié),內置單元的屬性需要包括密度,這是為了重力的計算與插值,彈性模量與泊松比的數值則為極小值,這就保證內置單元對方程左端的剛度矩陣無貢獻,即KF為極小量.而用戶自定義單元的屬性則只包括彈性模量與泊松比,密度數值則定義為極小量,因為方程左端的質量矩陣項已由MF提供,因此MS為極小值,這里用戶自定義單元只提供了剛度矩陣項KS,完成了剛度矩陣項與質量矩陣項的定義后,對于阻尼項,若采用瑞麗阻尼,則由剛度矩陣項與質量矩陣項線性組合而成.
用戶自定義單元與ABAQUS 內置單元庫結合示意如圖8所示.
圖8 用戶自定單元結合ABAQUS內置單元
用戶自定義單元與內置單元的結合在幾何上可以簡單理解為兩種幾何體相加的布爾運算.需要注意的是這兩種幾何體結合的界面不可以存在懸掛節(jié)點,否則會導致力學量插值出錯.背后的數學原理見式(5)、(6),所有的下標含義同式(4).
上式兩種原理表達式是可以互相轉化的,不同于給用戶自定義單元施加重力,內置單元與外置單元結合的時候二者的計算是互不關聯(lián)的.同時,二者的屬性定義也不同于給用戶自定義單元施加重力,這兩種單元所有的力學參數都需要正常定義,這種結合方式在豐富了ABAQUS內置單元庫的同時即保留了內置單元庫參與用戶自定義單元計算的可能,又增加了用戶自定義單元在ABAQUS中的適用性.
蒙皮技術只是將用戶自定義單元與ABAQUS內置功能結合的一種思路,用戶自定義單元與內置單元結合參與計算只是其中一個分支,以上只是介紹了其中兩種具體實施方法.另外,蒙皮技術的采用也解決了用戶自定義單元在ABAQUS中無法可視化的部分問題.用戶自定義單元無法在ABAQUS可視化在一定程度上影響了后處理工作,通過解析ODB 文件,將數據格式轉化為其他可視化軟件可識別的文件是一種思路,但是往往需要大量的編程,而且數據量較大的ODB文件解析非常耗時,因此探索一種簡便的用戶自定義單元可視化思路是必要的.蒙皮技術通過使用戶自定義單元與內置單元共享節(jié)點數據在一定程度上達到了位移可視化的目的,在后處理時只需要隱藏用戶自定義單元即可.
考慮精細梯形鍵槽模型,研究鍵槽受動荷載嚙合接觸分析,單個鍵槽平面如圖9所示,考慮帶有兩組鍵槽的塊體幾何模型,詳細尺寸如圖10所示.
圖9 單個鍵槽平面圖(單位:cm)
圖10 幾何模型尺寸圖(單位:cm)
如圖11所示,邊界面S1為加載面,S2為固定約束面,S3為兩塊體接觸面.采用六面體單元對兩個塊體進行離散,計算模型如圖12~13所示,其中圖12為協(xié)調網格,圖13為非協(xié)調網格.
圖11 加載、邊界面示意圖
圖12 協(xié)調網格離散示意圖
圖13 非協(xié)調網格離散示意圖
兩塊體的力學參數見表1.
表1 模型力學參數
在模型邊界上施加三角位移荷載,該位移荷載幅值如圖14所示,考慮模型尺寸,最終加載面上的位移值為時程圖中的2倍,施加的位移量為幅值曲線的2倍,即峰值位移荷載為2 cm,荷載作用時長為10 s,摩擦因數μ=0.6,采用動力隱式算法,增量步時長自適應.
圖14 荷載時程(時域)
對該模型進行加載,將位移荷載沿垂直鍵槽方向施加,即垂直于加載面S1,對加載情況下的模型力學響應進行分析.設置4種工況分別為:
工況1:采用ABAQUS自帶的單元,使用協(xié)調網格計算;
工況2:采用基于SBFEM 的用戶自定義單元,使用協(xié)調網格計算,通過蒙皮技術實現接觸屬性定義;
工況3:采用ABAQUS自帶的單元,使用非協(xié)調網格計算;
工況4:采用基于SBFEM 的用戶自定義單元,使用非協(xié)調網格計算,通過蒙皮技術實現接觸屬性定義.
在接觸面上設置3處參考點,如圖15所示.
圖15 觀測點位置示意圖
提取參考點的位移時程曲線和接觸壓力時程曲線如圖16所示.發(fā)現計算結果與ABAQUS 自帶單元吻合較好,采用蒙皮技術后,用戶自定義單元的可視化效果與ABAQUS自帶單元保持一致.通過該算例驗證了耦合SBFEM-FEM 的非線性分析方法可以用于結構動力計算.
圖16 觀測點動力響應
圖17 1 s時刻位移云圖
以NG5拱壩為研究對象,NG5水電站對南俄河的控制流域達483 km3,占總河流域的2.9%,最大壩高99.00 m,壩頂寬度6.0 m,拱冠梁底寬42.00 m,心線弧長234.838 m,弧高比2.42,拱頂中心角度92.794°,該水電站總庫容為3.41 億m3,死庫容為0.63 m3,水庫正常蓄水時候的庫容量為3.14億m3,庫容系數達34.9%.依據《水工建筑物抗震設計標準》[25],該水利工程按照500年一遇洪水(P=0.2%)涉及,按照10 000年一遇洪水(P=0.01%)校核,水庫正常蓄水位高度1 100.00 m.該拱壩平面布置如圖18所示.
圖18 NG5拱壩平面布置圖
建立了考慮橫縫接觸非線性的混凝土壩體-庫水-地基動力相互作用模型.其中壩體部分設置一條中縫,幾何模型示意如圖19所示.采用六面體單元離散,壩體與壩基采用實體C3D8單元模擬,庫水采用聲學單元AC3D8模擬,計算網格如圖20所示,計算模型網格信息見表2.
表2 計算模型網格單元信息
圖19 壩體-庫水-地基動力相互作用幾何模型
圖20 NG5拱壩計算模型(單位:m)
模型材料力學參數見表3.采用黏彈性人工邊界來模擬遠域地基的輻射阻尼效應,以此來考慮壩體-地基相互作用,黏彈性人工邊界的物理原件示意如圖21所示.對于橫縫非線性,采用ABAQUS中的指數接觸模型,其模型表達如圖22 所示.接觸界面采用ABAQUS內部的面-面接觸模型定義界面接觸屬性,接觸截面摩擦因數μ=0.6,根據大量試算[26],界面參數c0應取適當微量,p0取值應足夠大,那么在動力計算中可保證界面不發(fā)生相互貫入現象,文中取c0=5μm,p0=50 GPa.
表3 材料力學參數
圖21 黏彈性邊界示意圖
圖22 指數接觸模型示意圖
動荷載考慮地震作用,采用人工擬合地震動,峰值加速度為0.93 m/s2,輸入地震動的加速度時程曲線如圖23所示.
為兩種計算模型設計3種工況,見表4,其中工況3為無質量地基對照模型.
表4 計算工況信息
在地震作用下,對壩體位移、橫縫開度以及動水壓力等結果進行了對比分析.在壩體關鍵位置設置若干參考點,參考點的位置如圖24所示.
圖24 參考點位置
提取了壩體頂部的橫縫開度時程如圖25所示,3種工況下,橫縫開度峰值分別為2.21、2.70、5.27 mm.對比工況1和工況2發(fā)現,兩組數據吻合良好.并且采用了黏彈性人工邊界后,壩體的橫縫開度時程響應較無質量地基模型有較大降幅,降幅為58.06%.
圖25 點A-B橫縫開度時程曲線
繪制壩體壩頂與壩踵部位順河向相對位移曲線如圖26所示,提取順河向相對位移與相對速度峰值見表5,并給出工況2、3較工況1的增長幅度.
表5 壩頂-壩踵順河向峰值數據
圖26 順河向相對位移時程曲線
其中,工況3相對位移及相對速度峰值較工況1分別增長了115.44%、141.10%,說明了無質量地基模型相較有質量地基模型會放大結構地震響應.且工況1和工況2的相對位移峰值以及相對速度峰值的相對誤差分別為2.40%和2.44%,說明SBFEM 與FEM 兩種方法計算結果一致.
考慮壩面左岸和右岸中部位置的兩個關鍵點,提取了關鍵點部位的動水壓力時程曲線,如圖27~28所示,并匯總動水壓力峰值于表6.
表6 壩面中部關鍵點動水壓力峰值 (單位:Pa)
圖27 C點動水壓力時程曲線
圖28 D 點動水壓力時程曲線
對于工況1和工況2,點C、D 二者的動水壓力峰值相對誤差分別為0.54%、0.42%,結合時程曲線,發(fā)現兩種工況吻合良好,說明兩種算法的計算效果一致.且由表6可知無質量地基模型動水壓力峰值較有質量地基模型,具有放大效應.
本文結合比例邊界有限元法理論與商用程序ABAQUS的優(yōu)點,通過用戶子程序接口,將SBFEM嵌入到商用程序ABAQUS中.開發(fā)了一種不引入額外自由度的蒙皮單元技術.在此基礎上以NG5拱壩為研究對象,采用黏彈性動力人工邊界條件,考慮接觸非線性,建立了耦合比例邊界有限元法-經典有限元法的壩體-庫水-地基動力相互作用模型,進行非線性地震響應分析,可以得出以下結論:
1)本文提出蒙皮單元技術可以實現比例邊界有限元與ABAQUS自帶單元、模型等功能相結合,在保留了該理論計算性能的基礎上仍能發(fā)揮ABAQUS強大的非線性計算性能.
2)采用蒙皮技術后,用戶自定義單元的可視化效果與ABAQUS 自帶單元保持一致,驗證了耦合SBFEM-FEM 的非線性分析方法可以用于結構動力計算.
3)耦合SBFEM-FEM 方法,可以充分考慮大壩-庫水-地基的動力相互作用模型、橫縫接觸非線性以及地震動輸入.計算結果與ABAQUS自帶單元的結果吻合良好,計算效率高,為復雜工程結構的動力分析提供一種新的方法.