陳 汭,楊海霞,賈金生,鄭璀瑩
(1.河海大學(xué) 力學(xué)與材料學(xué)院,江蘇 南京 211100; 2.中國水利水電科學(xué)研究院,北京 100038)
膠結(jié)砂礫石壩是區(qū)別混凝土壩和土石壩的新壩型[1],有就地取材方便、安全經(jīng)濟(jì)、漫頂不潰的優(yōu)點,同時能較好地適應(yīng)軟弱地基。 砂卵石土在我國西南、西北地區(qū)有較廣分布[2],探究砂卵石覆蓋層上膠結(jié)砂礫石壩的應(yīng)力變形特性、地基承載狀況,對該壩型的推廣有著積極作用[3]。 國內(nèi)關(guān)于膠結(jié)砂礫石壩的研究主要集中于材料特性和壩體穩(wěn)定性研究層面[4-7]。 賈金生等[6]進(jìn)行了膠結(jié)砂礫石壩斷面形式和材料配比對膠結(jié)砂礫石料力學(xué)性能影響的研究,楊世鋒等[8]利用有限元軟件對膠結(jié)砂礫石壩進(jìn)行了穩(wěn)定分析研究。但考慮覆蓋層上的膠結(jié)砂礫石壩應(yīng)力變形特征的研究較少,關(guān)于膠結(jié)砂礫石壩優(yōu)化設(shè)計的研究也較少。
筆者以一座設(shè)計階段的膠結(jié)砂礫石壩為例,其壩址處地基覆蓋層為砂卵石土層,壩體為膠結(jié)砂礫石料,在該工程設(shè)計初期采用材料力學(xué)法,進(jìn)行了結(jié)構(gòu)穩(wěn)定性校核,發(fā)現(xiàn)壩基應(yīng)力超過地基承載力。 由于膠結(jié)砂礫石料本構(gòu)關(guān)系呈非線性,材料力學(xué)法計算結(jié)果不夠準(zhǔn)確,因此需選擇鄧肯-張模型(經(jīng)典模型與修正模型)作為材料本構(gòu)模型,以了解結(jié)構(gòu)真實的應(yīng)力分布規(guī)律[8]。 筆者首先運用ANSYS 參數(shù)化設(shè)計命令流(APDL)定義非線性本構(gòu)模型,基于荷載分步施加進(jìn)行有限元計算,獲得并分析應(yīng)力變形規(guī)律,然后依據(jù)材料力學(xué)法與有限元分析所得結(jié)果,選取優(yōu)化設(shè)計目標(biāo),通過ANSYS 與MATLAB 的數(shù)據(jù)交互,優(yōu)化壩體斷面尺寸,達(dá)到應(yīng)力狀態(tài)的改良。
運用APDL 建模,順河向為X軸,垂直河向為Y軸,垂直方向為Z軸。 模型基于設(shè)計需求簡化后,分為壩體、防滲墻、覆蓋層土層、基巖四部分,其中壩體劃分為兩階,具體尺寸如圖1 所示。 模型在Z方向的厚度取單位厚度1 m,可視為假三維的二維模型。 模型下界面添加全約束,壩基在X方向的兩個界面添加X方向約束,模型在Z方向的兩個剖面添加Z方向約束。 模型采用映射網(wǎng)格劃分,共劃分4 910 個網(wǎng)格單元,如圖2 所示。
圖1 壩體橫斷面與分區(qū)(單位:m)
圖2 網(wǎng)格單元劃分
根據(jù)廣義胡克定律建立彈性常數(shù)矩陣D,其包含的彈性模量E和泊松比μ都不是常量,而是隨應(yīng)力狀態(tài)改變的量[9]。 覆蓋層為砂卵石土層,采用鄧肯-張經(jīng)典模型[2]。 該工程壩體材料為膠結(jié)砂礫石料,應(yīng)力應(yīng)變曲線具有雙曲線特征,但與傳統(tǒng)土石壩不同,需采用經(jīng)虛加剛性彈簧法修正的鄧肯-張模型(后文簡稱鄧肯-張修正模型),該模型解決了應(yīng)力應(yīng)變曲線的軟化問題,且與其他膠結(jié)砂礫石材料本構(gòu)模型相比,在有限元計算中具有優(yōu)勢[10-14]。
1.2.1 鄧肯-張經(jīng)典模型
判斷單元狀態(tài):(σ1- σ3)<(σ1- σ3)0且S <S0時,單元處于卸載狀態(tài),此時彈性模量用Eur表示;反之,單元處于加載狀態(tài),此時彈性模量用Et表示。 其中:σ1和σ3分別為單元的大主應(yīng)力、小主應(yīng)力;(σ1-σ3)0為歷史上單元達(dá)到的最大偏應(yīng)力狀態(tài);S為單元的應(yīng)力水平;S0為歷史上單元達(dá)到的最大應(yīng)力水平。
對于加載狀態(tài),切線彈性模量Et表達(dá)式為
式中:c為材料凝聚力;φ為材料內(nèi)摩擦角;Pa為標(biāo)準(zhǔn)大氣壓力;Rf為破壞比;K為彈性模量基數(shù);n為彈性模量指數(shù)。
對于卸載狀態(tài),切線彈性模量Eur表達(dá)式為
式中:Kur為卸載和再加載時的彈性模量基數(shù);nur為卸載和再加載時的彈性模量指數(shù)。
模型的泊松比μ表達(dá)式為
其中
式中:G為切線泊松比基數(shù);F為反映初始切線泊松比隨小主應(yīng)力減小的參數(shù);D為側(cè)向應(yīng)變與軸向應(yīng)變關(guān)系曲線上小漸進(jìn)值的倒數(shù)。
1.2.2 基于虛加剛性彈簧法的鄧肯-張修正模型
圖3(a)中①為實際應(yīng)力應(yīng)變曲線,②為彈簧應(yīng)力應(yīng)變線性關(guān)系,在適當(dāng)位置加虛擬彈簧(如圖3(b)所示),①和②疊加即得曲線③,為雙曲線,沒有應(yīng)力軟化階段,可用鄧肯-張模型表示。 則膠結(jié)砂礫石真實切線彈性模量為
圖3 3 種應(yīng)力應(yīng)變曲線和虛加彈簧示意
式中:E1為實際膠結(jié)砂礫石的切線彈性模量;E2為虛加剛性彈簧的彈性模量,取實際應(yīng)力應(yīng)變曲線最大負(fù)彈性模量的絕對值;E3為虛擬雙曲線的切線彈性模量[11]。
模型參數(shù)經(jīng)過虛加剛性彈簧法處理后便可適用于膠結(jié)砂礫石料。 修正模型與經(jīng)典模型在表達(dá)式上一致,但不區(qū)分加載、卸載的參數(shù)差別,即切線彈性模量Et表達(dá)式為式(1)、泊松比μ表達(dá)式為式(3)、A的表達(dá)式為式(4)。
可以通過三軸固結(jié)排水試驗得到2 種材料本構(gòu)模型參數(shù)(見表1)。
表1 本構(gòu)模型參數(shù)
1.3.1 建立鄧肯-張宏命令
創(chuàng)建本構(gòu)關(guān)系的宏命令文件后,可使各單元在每個荷載步后自動根據(jù)應(yīng)力狀態(tài)更新彈性常數(shù),簡化計算流程[15]。
1.3.2 荷載逐級加載的模擬實現(xiàn)[9,16-17]
逐級加載與瞬間加載效果不同,非線性計算更需分步加載。 計算荷載施加分為9 步,模擬壩體填筑和蓄水過程,填筑過程考慮重力荷載,蓄水過程還考慮靜水壓力與揚壓力,通過ANSYS 的生死單元命令實現(xiàn)該功能。 計算流程如圖4 所示。
圖4 ANSYS 計算流程
1.3.3 重啟動技術(shù)
ANSYS 非線性計算在荷載步求解完畢后,都需將文件存儲目錄里的“rst”“rxxx”“l(fā)dhi”這3 個文件刪除,通過單點重啟動保證后續(xù)計算正常運行[18]。
1.3.4 約束荷載平衡法
在荷載步計算前讀入初始應(yīng)力文件,產(chǎn)生與原荷載大小相同而方向相反的附加位移,使結(jié)構(gòu)初始位移為零,消除基巖區(qū)因自重沉降而產(chǎn)生的影響[19-20]。
在工程設(shè)計初期,通過材料力學(xué)法進(jìn)行了壩體穩(wěn)定性校核,結(jié)果見表2,可知壩體抗滑穩(wěn)定是安全的,但壩踵處應(yīng)力超過了砂卵石地基的允許承載力,所以設(shè)計單位希望通過有限元分析進(jìn)行復(fù)核。 一方面可以通過有限元分析實現(xiàn)非線性本構(gòu)的定義,計算所得應(yīng)力結(jié)果更符合結(jié)構(gòu)真實應(yīng)力狀態(tài),結(jié)構(gòu)各特征值極值也可為安全性判斷提供參考;另一方面可依據(jù)材料力學(xué)法與有限元分析結(jié)果,選取優(yōu)化設(shè)計目標(biāo),通過調(diào)整斷面尺寸,改善結(jié)構(gòu)應(yīng)力狀態(tài)。
表2 材料力學(xué)法校核結(jié)果
由于單元材料間的彈性常數(shù)存在差異,因此ANSYS 呈現(xiàn)結(jié)果云圖時,非線性材料區(qū)呈網(wǎng)格狀,線彈性材料區(qū)則呈片狀。 防滲墻會對揚壓力和單元彈性常數(shù)的連續(xù)性產(chǎn)生影響,有限元建模時對防滲墻的存在進(jìn)行了保留,結(jié)果分析仍只討論壩體和覆蓋層結(jié)構(gòu)。
選取竣工期與蓄水期結(jié)果進(jìn)行對比,關(guān)鍵特征值的極值見表3。 對于X向位移,結(jié)構(gòu)極值都不超過1 cm,本文不再分析該特征值。 對于Y方向沉降,壩體極值因揚壓力作用更顯著而變小,覆蓋層極值因靜水壓力的傳遞效果而變大,但數(shù)值都不超過6.5 cm,不影響結(jié)構(gòu)安全。 對于最大拉應(yīng)力,壩體擠壓覆蓋層,在接觸面處產(chǎn)生了極值,該數(shù)值不影響結(jié)構(gòu)安全。 對于最大壓應(yīng)力,壩體極值為竣工期的2.78 MPa,小于壩體材料的抗壓強(qiáng)度(6 MPa),位置變化原因是上游水壓力改變了壩體一、二階結(jié)構(gòu)接觸點的擠壓效果,使上游接觸點擠壓減小而下游接觸點擠壓增大;覆蓋層極值為蓄水期時的1.09 MPa,極限點位置與材料力學(xué)法結(jié)果一致,都為壩踵處。 考慮到剪切破壞是土體主要破壞形式,提取蓄水期覆蓋層各節(jié)點抗剪強(qiáng)度與最大剪應(yīng)力差值相對于抗剪強(qiáng)度的比值(稱為抗剪強(qiáng)度相對值),該值最大為-0.038,可知覆蓋層存在剪切破壞點。
表3 特征值極值及位置
靜水壓力和揚壓力對結(jié)構(gòu)特征值極值影響不大:一是壩體本身斷面較厚,且與地基接觸面寬,壩型穩(wěn)定;二是蓄水深度不大,且因防滲墻的存在而對揚壓力進(jìn)行了折減,兩種荷載作用有限。
優(yōu)化設(shè)計涉及目標(biāo)函數(shù)、設(shè)計變量、約束條件3 個要素,需通過數(shù)學(xué)建模,將優(yōu)化問題轉(zhuǎn)換為求函數(shù)極值問題。
3.1.1 目標(biāo)函數(shù)
由材料力學(xué)法結(jié)果可知,優(yōu)化目標(biāo)應(yīng)包含地基壓應(yīng)力。 根據(jù)有限元分析結(jié)果可知,結(jié)構(gòu)變形、結(jié)構(gòu)受拉、壩體受壓都在安全范圍內(nèi),但覆蓋層存在剪切破壞點,且其受壓極值在有限元法中并無具體控制值。 考慮蓄水期靜水壓力會加大地基受力,本優(yōu)化設(shè)計目標(biāo)定為蓄水期覆蓋層最大壓應(yīng)力絕對值、覆蓋層抗剪強(qiáng)度相對值,通過提取各節(jié)點的目標(biāo)數(shù)值進(jìn)行大小排序,對最劣點進(jìn)行優(yōu)化,實現(xiàn)地基應(yīng)力狀態(tài)的改良。
3.1.2 設(shè)計變量
設(shè)計變量為壩體一、二階的上下游坡比,即X=[X1,X2,X3,X4]T,其中X1和X2為壩體二階結(jié)構(gòu)的下游和上游坡比,X3和X4為壩體一階結(jié)構(gòu)的下游和上游坡比。
3.1.3 約束條件
優(yōu)化設(shè)計的約束條件參考《膠結(jié)顆粒料筑壩導(dǎo)則》《混凝土重力壩設(shè)計規(guī)范》《水閘設(shè)計規(guī)范》[21-23],結(jié)合材料本構(gòu)特性與設(shè)計需求,設(shè)定如下。
(1)幾何約束條件。 ①壩體一階的上下游坡比范圍取1 ∶1.5~1 ∶0.6,壩體二階的上下游坡比范圍分別取1 ∶0.4~1 ∶0.7 和1 ∶0.5 ~1 ∶0.8;②優(yōu)化壩型斷面面積S優(yōu)化不大于初始壩型斷面面積S初始。
(2)應(yīng)力約束條件。 ①壩踵垂直應(yīng)力不出現(xiàn)拉應(yīng)力;②壩體上游面不出現(xiàn)拉應(yīng)力;③壩體最大主壓應(yīng)力不大于混凝土的最大允許壓應(yīng)力;④壩踵與壩趾的應(yīng)力比值不大于允許值,該允許值取1.5。
(3)穩(wěn)定約束條件。 各種荷載組合下壩基面的最小抗滑安全系數(shù)Ks不小于允許值[Ks], [Ks]取1.05。
APDL 本身具備優(yōu)化功能,但在進(jìn)程控制和邏輯框架編制上,操作難度較大[24],且優(yōu)化算法只有零階和一階兩種,在優(yōu)化算法種類、問題解決適用性上,不如具備豐富優(yōu)化工具箱的MATLAB。 因此,本優(yōu)化設(shè)計運用了MATLAB 和ANSYS 的混編,即優(yōu)化主模塊在MATLAB 中進(jìn)行,將ANSYS 有限元分析視為計算子程序進(jìn)行調(diào)用。 軟件混編的實現(xiàn),結(jié)合了MATLAB強(qiáng)大的優(yōu)化功能和ANSYS 良好的有限元分析性能,提供了一種有限元優(yōu)化設(shè)計的借鑒思路。 數(shù)據(jù)交互流程如圖5 所示。 MATLAB 優(yōu)化計算使用MATLAB 優(yōu)化工具箱中全局優(yōu)化工具Patten Searsh。
圖5 MATLAB 與ANSYS 數(shù)據(jù)交互流程
表4 列出了初始壩型和優(yōu)化壩型的特征值,優(yōu)化壩型一和優(yōu)化壩型二分別為覆蓋層抗壓、抗剪優(yōu)化壩型。 優(yōu)化壩型一和優(yōu)化壩型二斷面對比見圖6、主應(yīng)力云圖對比見圖7。
圖6 優(yōu)化壩型一和優(yōu)化壩型二斷面對比
圖7 優(yōu)化壩型一和優(yōu)化壩型二主應(yīng)力云圖對比(單位:Pa)
表4 各特征值優(yōu)化效果
優(yōu)化壩型一的4 個坡比都顯著增大,優(yōu)化原理是通過減小壩體體量,降低覆蓋層受壓效果。 優(yōu)化壩型二的二階坡比增大而一階坡比變化很小,優(yōu)化原理是通過地基接觸面寬度與斷面尺寸的平衡,實現(xiàn)了覆蓋層抗剪狀態(tài)的改善。
解決材料力學(xué)法結(jié)果體現(xiàn)的地基壓應(yīng)力過大問題,是優(yōu)化設(shè)計的目的之一。 在覆蓋層壓應(yīng)力優(yōu)化上,優(yōu)化壩型一效果更好,且優(yōu)化壩型一、壩型二的極值點位置都在壩踵處,契合材料力學(xué)法所得結(jié)果,具有設(shè)計指導(dǎo)意義。 剪切破壞雖是土體主要破壞形式,但優(yōu)化壩型一抗剪狀態(tài)特征值為正值,可知實現(xiàn)了土體剪切破壞點數(shù)量的清零。 在其余特征值的對比上,優(yōu)化壩型一都比優(yōu)化壩型二更具優(yōu)勢。 綜上,優(yōu)化壩型一實現(xiàn)了斷面面積減少11.3%、覆蓋層最大壓應(yīng)力降低9.2%和抗剪破壞點清零,達(dá)到了材料用量減少和地基應(yīng)力狀態(tài)改良的目的,更符合設(shè)計需求。
首先,本文從膠結(jié)砂礫石料和砂卵石土的材料本構(gòu)出發(fā),介紹了鄧肯-張經(jīng)典與修正模型通過ANSYS實現(xiàn)的思路,并簡要介紹了ANSYS 二次開發(fā)的細(xì)節(jié),基于此,對初設(shè)尺寸的壩型進(jìn)行了非線性有限元計算,獲得了竣工期與蓄水期的特征值,分析其應(yīng)力變形規(guī)律;其次,依據(jù)材料力學(xué)法和有限元分析結(jié)果,以覆蓋層抗壓、抗剪性能作為優(yōu)化設(shè)計目標(biāo),運用數(shù)據(jù)交互方式,實現(xiàn)ANSYS 與MATLAB 的混編,得到兩種優(yōu)化壩型;最后,將優(yōu)化壩型與初始壩型的各項特征值進(jìn)行對比,評判優(yōu)化效果。
優(yōu)化目標(biāo)的設(shè)定應(yīng)符合工程設(shè)計需求,設(shè)計變量也不只局限于壩型尺寸,可擴(kuò)展到覆蓋層尺寸與結(jié)構(gòu)材料參數(shù)。 本文探究了ANSYS 有限元分析與MATLAB優(yōu)化計算結(jié)合的思路,所得優(yōu)化壩型與特征值可供相關(guān)設(shè)計單位借鑒,后期會針對不同的設(shè)計變量對各目標(biāo)函數(shù)的敏感度做進(jìn)一步的優(yōu)化設(shè)計。