基于ANSYS 和MATLAB 的膠結(jié)砂礫石壩斷面優(yōu)化

陳 汭，楊海霞，賈金生，鄭璀瑩

（1.河海大學(xué) 力學(xué)與材料學(xué)院，江蘇 南京 211100； 2.中國(guó)水利水電科學(xué)研究院，北京 100038）

膠結(jié)砂礫石壩是區(qū)別混凝土壩和土石壩的新壩型［1］，有就地取材方便、安全經(jīng)濟(jì)、漫頂不潰的優(yōu)點(diǎn)，同時(shí)能較好地適應(yīng)軟弱地基。 砂卵石土在我國(guó)西南、西北地區(qū)有較廣分布［2］，探究砂卵石覆蓋層上膠結(jié)砂礫石壩的應(yīng)力變形特性、地基承載狀況，對(duì)該壩型的推廣有著積極作用［3］。 國(guó)內(nèi)關(guān)于膠結(jié)砂礫石壩的研究主要集中于材料特性和壩體穩(wěn)定性研究層面［4－7］。 賈金生等［6］進(jìn)行了膠結(jié)砂礫石壩斷面形式和材料配比對(duì)膠結(jié)砂礫石料力學(xué)性能影響的研究，楊世鋒等［8］利用有限元軟件對(duì)膠結(jié)砂礫石壩進(jìn)行了穩(wěn)定分析研究。但考慮覆蓋層上的膠結(jié)砂礫石壩應(yīng)力變形特征的研究較少，關(guān)于膠結(jié)砂礫石壩優(yōu)化設(shè)計(jì)的研究也較少。

筆者以一座設(shè)計(jì)階段的膠結(jié)砂礫石壩為例，其壩址處地基覆蓋層為砂卵石土層，壩體為膠結(jié)砂礫石料，在該工程設(shè)計(jì)初期采用材料力學(xué)法，進(jìn)行了結(jié)構(gòu)穩(wěn)定性校核，發(fā)現(xiàn)壩基應(yīng)力超過(guò)地基承載力。 由于膠結(jié)砂礫石料本構(gòu)關(guān)系呈非線性，材料力學(xué)法計(jì)算結(jié)果不夠準(zhǔn)確，因此需選擇鄧肯－張模型（經(jīng)典模型與修正模型）作為材料本構(gòu)模型，以了解結(jié)構(gòu)真實(shí)的應(yīng)力分布規(guī)律［8］。 筆者首先運(yùn)用ANSYS 參數(shù)化設(shè)計(jì)命令流（APDL）定義非線性本構(gòu)模型，基于荷載分步施加進(jìn)行有限元計(jì)算，獲得并分析應(yīng)力變形規(guī)律，然后依據(jù)材料力學(xué)法與有限元分析所得結(jié)果，選取優(yōu)化設(shè)計(jì)目標(biāo)，通過(guò)ANSYS 與MATLAB 的數(shù)據(jù)交互，優(yōu)化壩體斷面尺寸，達(dá)到應(yīng)力狀態(tài)的改良。

1 有限元模型與材料本構(gòu)關(guān)系

1.1 二維有限元模型的建立

運(yùn)用APDL 建模，順河向?yàn)閄軸，垂直河向?yàn)閅軸，垂直方向?yàn)閆軸。 模型基于設(shè)計(jì)需求簡(jiǎn)化后，分為壩體、防滲墻、覆蓋層土層、基巖四部分，其中壩體劃分為兩階，具體尺寸如圖1 所示。 模型在Z方向的厚度取單位厚度1 m，可視為假三維的二維模型。 模型下界面添加全約束，壩基在X方向的兩個(gè)界面添加X(jué)方向約束，模型在Z方向的兩個(gè)剖面添加Z方向約束。 模型采用映射網(wǎng)格劃分，共劃分4 910 個(gè)網(wǎng)格單元，如圖2 所示。

圖1 壩體橫斷面與分區(qū)（單位：m）

圖2 網(wǎng)格單元?jiǎng)澐?/p>

1.2 鄧肯－張本構(gòu)模型

根據(jù)廣義胡克定律建立彈性常數(shù)矩陣D，其包含的彈性模量E和泊松比μ都不是常量，而是隨應(yīng)力狀態(tài)改變的量［9］。 覆蓋層為砂卵石土層，采用鄧肯－張經(jīng)典模型［2］。 該工程壩體材料為膠結(jié)砂礫石料，應(yīng)力應(yīng)變曲線具有雙曲線特征，但與傳統(tǒng)土石壩不同，需采用經(jīng)虛加剛性彈簧法修正的鄧肯－張模型（后文簡(jiǎn)稱鄧肯－張修正模型），該模型解決了應(yīng)力應(yīng)變曲線的軟化問(wèn)題，且與其他膠結(jié)砂礫石材料本構(gòu)模型相比，在有限元計(jì)算中具有優(yōu)勢(shì)［10－14］。

1.2.1 鄧肯－張經(jīng)典模型

判斷單元狀態(tài)：（σ1－ σ3）＜（σ1－ σ3）0且S ＜S0時(shí)，單元處于卸載狀態(tài)，此時(shí)彈性模量用Eur表示；反之，單元處于加載狀態(tài)，此時(shí)彈性模量用Et表示。 其中：σ1和σ3分別為單元的大主應(yīng)力、小主應(yīng)力；（σ1－σ3）0為歷史上單元達(dá)到的最大偏應(yīng)力狀態(tài)；S為單元的應(yīng)力水平；S0為歷史上單元達(dá)到的最大應(yīng)力水平。

對(duì)于加載狀態(tài)，切線彈性模量Et表達(dá)式為

式中：c為材料凝聚力；φ為材料內(nèi)摩擦角；Pa為標(biāo)準(zhǔn)大氣壓力；Rf為破壞比；K為彈性模量基數(shù)；n為彈性模量指數(shù)。

對(duì)于卸載狀態(tài)，切線彈性模量Eur表達(dá)式為

式中：Kur為卸載和再加載時(shí)的彈性模量基數(shù)；nur為卸載和再加載時(shí)的彈性模量指數(shù)。

模型的泊松比μ表達(dá)式為

其中

式中：G為切線泊松比基數(shù)；F為反映初始切線泊松比隨小主應(yīng)力減小的參數(shù)；D為側(cè)向應(yīng)變與軸向應(yīng)變關(guān)系曲線上小漸進(jìn)值的倒數(shù)。

1.2.2 基于虛加剛性彈簧法的鄧肯－張修正模型

圖3（a）中①為實(shí)際應(yīng)力應(yīng)變曲線，②為彈簧應(yīng)力應(yīng)變線性關(guān)系，在適當(dāng)位置加虛擬彈簧（如圖3（b）所示），①和②疊加即得曲線③，為雙曲線，沒(méi)有應(yīng)力軟化階段，可用鄧肯－張模型表示。 則膠結(jié)砂礫石真實(shí)切線彈性模量為

圖3 3 種應(yīng)力應(yīng)變曲線和虛加彈簧示意

式中：E1為實(shí)際膠結(jié)砂礫石的切線彈性模量；E2為虛加剛性彈簧的彈性模量，取實(shí)際應(yīng)力應(yīng)變曲線最大負(fù)彈性模量的絕對(duì)值；E3為虛擬雙曲線的切線彈性模量［11］。

模型參數(shù)經(jīng)過(guò)虛加剛性彈簧法處理后便可適用于膠結(jié)砂礫石料。 修正模型與經(jīng)典模型在表達(dá)式上一致，但不區(qū)分加載、卸載的參數(shù)差別，即切線彈性模量Et表達(dá)式為式（1）、泊松比μ表達(dá)式為式（3）、A的表達(dá)式為式（4）。

可以通過(guò)三軸固結(jié)排水試驗(yàn)得到2 種材料本構(gòu)模型參數(shù)（見(jiàn)表1）。

表1 本構(gòu)模型參數(shù)

1.3 ANSYS 二次開(kāi)發(fā)

1.3.1 建立鄧肯－張宏命令

創(chuàng)建本構(gòu)關(guān)系的宏命令文件后，可使各單元在每個(gè)荷載步后自動(dòng)根據(jù)應(yīng)力狀態(tài)更新彈性常數(shù)，簡(jiǎn)化計(jì)算流程［15］。

1.3.2 荷載逐級(jí)加載的模擬實(shí)現(xiàn)［9，16－17］

逐級(jí)加載與瞬間加載效果不同，非線性計(jì)算更需分步加載。 計(jì)算荷載施加分為9 步，模擬壩體填筑和蓄水過(guò)程，填筑過(guò)程考慮重力荷載，蓄水過(guò)程還考慮靜水壓力與揚(yáng)壓力，通過(guò)ANSYS 的生死單元命令實(shí)現(xiàn)該功能。 計(jì)算流程如圖4 所示。

圖4 ANSYS 計(jì)算流程

1.3.3 重啟動(dòng)技術(shù)

ANSYS 非線性計(jì)算在荷載步求解完畢后，都需將文件存儲(chǔ)目錄里的“rst”“rxxx”“l(fā)dhi”這3 個(gè)文件刪除，通過(guò)單點(diǎn)重啟動(dòng)保證后續(xù)計(jì)算正常運(yùn)行［18］。

1.3.4 約束荷載平衡法

在荷載步計(jì)算前讀入初始應(yīng)力文件，產(chǎn)生與原荷載大小相同而方向相反的附加位移，使結(jié)構(gòu)初始位移為零，消除基巖區(qū)因自重沉降而產(chǎn)生的影響［19－20］。

2 應(yīng)力與變形分析

2.1 材料力學(xué)法應(yīng)力校核

在工程設(shè)計(jì)初期，通過(guò)材料力學(xué)法進(jìn)行了壩體穩(wěn)定性校核，結(jié)果見(jiàn)表2，可知壩體抗滑穩(wěn)定是安全的，但壩踵處應(yīng)力超過(guò)了砂卵石地基的允許承載力，所以設(shè)計(jì)單位希望通過(guò)有限元分析進(jìn)行復(fù)核。 一方面可以通過(guò)有限元分析實(shí)現(xiàn)非線性本構(gòu)的定義，計(jì)算所得應(yīng)力結(jié)果更符合結(jié)構(gòu)真實(shí)應(yīng)力狀態(tài)，結(jié)構(gòu)各特征值極值也可為安全性判斷提供參考；另一方面可依據(jù)材料力學(xué)法與有限元分析結(jié)果，選取優(yōu)化設(shè)計(jì)目標(biāo)，通過(guò)調(diào)整斷面尺寸，改善結(jié)構(gòu)應(yīng)力狀態(tài)。

表2 材料力學(xué)法校核結(jié)果

2.2 有限元分析結(jié)果

由于單元材料間的彈性常數(shù)存在差異，因此ANSYS 呈現(xiàn)結(jié)果云圖時(shí)，非線性材料區(qū)呈網(wǎng)格狀，線彈性材料區(qū)則呈片狀。 防滲墻會(huì)對(duì)揚(yáng)壓力和單元彈性常數(shù)的連續(xù)性產(chǎn)生影響，有限元建模時(shí)對(duì)防滲墻的存在進(jìn)行了保留，結(jié)果分析仍只討論壩體和覆蓋層結(jié)構(gòu)。

選取竣工期與蓄水期結(jié)果進(jìn)行對(duì)比，關(guān)鍵特征值的極值見(jiàn)表3。 對(duì)于X向位移，結(jié)構(gòu)極值都不超過(guò)1 cm，本文不再分析該特征值。 對(duì)于Y方向沉降，壩體極值因揚(yáng)壓力作用更顯著而變小，覆蓋層極值因靜水壓力的傳遞效果而變大，但數(shù)值都不超過(guò)6.5 cm，不影響結(jié)構(gòu)安全。 對(duì)于最大拉應(yīng)力，壩體擠壓覆蓋層，在接觸面處產(chǎn)生了極值，該數(shù)值不影響結(jié)構(gòu)安全。 對(duì)于最大壓應(yīng)力，壩體極值為竣工期的2.78 MPa，小于壩體材料的抗壓強(qiáng)度（6 MPa），位置變化原因是上游水壓力改變了壩體一、二階結(jié)構(gòu)接觸點(diǎn)的擠壓效果，使上游接觸點(diǎn)擠壓減小而下游接觸點(diǎn)擠壓增大；覆蓋層極值為蓄水期時(shí)的1.09 MPa，極限點(diǎn)位置與材料力學(xué)法結(jié)果一致，都為壩踵處。 考慮到剪切破壞是土體主要破壞形式，提取蓄水期覆蓋層各節(jié)點(diǎn)抗剪強(qiáng)度與最大剪應(yīng)力差值相對(duì)于抗剪強(qiáng)度的比值（稱為抗剪強(qiáng)度相對(duì)值），該值最大為－0.038，可知覆蓋層存在剪切破壞點(diǎn)。

表3 特征值極值及位置

靜水壓力和揚(yáng)壓力對(duì)結(jié)構(gòu)特征值極值影響不大：一是壩體本身斷面較厚，且與地基接觸面寬，壩型穩(wěn)定；二是蓄水深度不大，且因防滲墻的存在而對(duì)揚(yáng)壓力進(jìn)行了折減，兩種荷載作用有限。

3 優(yōu)化設(shè)計(jì)

3.1 優(yōu)化設(shè)計(jì)要素

優(yōu)化設(shè)計(jì)涉及目標(biāo)函數(shù)、設(shè)計(jì)變量、約束條件3 個(gè)要素，需通過(guò)數(shù)學(xué)建模，將優(yōu)化問(wèn)題轉(zhuǎn)換為求函數(shù)極值問(wèn)題。

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)，但覆蓋層存在剪切破壞點(diǎn)，且其受壓極值在有限元法中并無(wú)具體控制值。 考慮蓄水期靜水壓力會(huì)加大地基受力，本優(yōu)化設(shè)計(jì)目標(biāo)定為蓄水期覆蓋層最大壓應(yīng)力絕對(duì)值、覆蓋層抗剪強(qiáng)度相對(duì)值，通過(guò)提取各節(jié)點(diǎn)的目標(biāo)數(shù)值進(jìn)行大小排序，對(duì)最劣點(diǎn)進(jìn)行優(yōu)化，實(shí)現(xiàn)地基應(yīng)力狀態(tài)的改良。

3.1.2 設(shè)計(jì)變量

設(shè)計(jì)變量為壩體一、二階的上下游坡比，即X＝［X1，X2，X3，X4］T，其中X1和X2為壩體二階結(jié)構(gòu)的下游和上游坡比，X3和X4為壩體一階結(jié)構(gòu)的下游和上游坡比。

3.1.3 約束條件

優(yōu)化設(shè)計(jì)的約束條件參考《膠結(jié)顆粒料筑壩導(dǎo)則》《混凝土重力壩設(shè)計(jì)規(guī)范》《水閘設(shè)計(jì)規(guī)范》［21－23］，結(jié)合材料本構(gòu)特性與設(shè)計(jì)需求，設(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)化不大于初始?jí)涡蛿嗝婷娣eS初始。

（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。

3.2 MATLAB 優(yōu)化計(jì)算與ANSYS 有限元分析的交互

APDL 本身具備優(yōu)化功能，但在進(jìn)程控制和邏輯框架編制上，操作難度較大［24］，且優(yōu)化算法只有零階和一階兩種，在優(yōu)化算法種類、問(wèn)題解決適用性上，不如具備豐富優(yōu)化工具箱的MATLAB。 因此，本優(yōu)化設(shè)計(jì)運(yùn)用了MATLAB 和ANSYS 的混編，即優(yōu)化主模塊在MATLAB 中進(jìn)行，將ANSYS 有限元分析視為計(jì)算子程序進(jìn)行調(diào)用。 軟件混編的實(shí)現(xiàn)，結(jié)合了MATLAB強(qiáng)大的優(yōu)化功能和ANSYS 良好的有限元分析性能，提供了一種有限元優(yōu)化設(shè)計(jì)的借鑒思路。 數(shù)據(jù)交互流程如圖5 所示。 MATLAB 優(yōu)化計(jì)算使用MATLAB 優(yōu)化工具箱中全局優(yōu)化工具Patten Searsh。

圖5 MATLAB 與ANSYS 數(shù)據(jù)交互流程

3.3 優(yōu)化結(jié)果分析

表4 列出了初始?jí)涡秃蛢?yōu)化壩型的特征值，優(yōu)化壩型一和優(yōu)化壩型二分別為覆蓋層抗壓、抗剪優(yōu)化壩型。 優(yōu)化壩型一和優(yōu)化壩型二斷面對(duì)比見(jiàn)圖6、主應(yīng)力云圖對(duì)比見(jiàn)圖7。

圖6 優(yōu)化壩型一和優(yōu)化壩型二斷面對(duì)比

圖7 優(yōu)化壩型一和優(yōu)化壩型二主應(yīng)力云圖對(duì)比（單位：Pa）

表4 各特征值優(yōu)化效果

3.4 優(yōu)化方案評(píng)估

優(yōu)化壩型一的4 個(gè)坡比都顯著增大，優(yōu)化原理是通過(guò)減小壩體體量，降低覆蓋層受壓效果。 優(yōu)化壩型二的二階坡比增大而一階坡比變化很小，優(yōu)化原理是通過(guò)地基接觸面寬度與斷面尺寸的平衡，實(shí)現(xiàn)了覆蓋層抗剪狀態(tài)的改善。

解決材料力學(xué)法結(jié)果體現(xiàn)的地基壓應(yīng)力過(guò)大問(wèn)題，是優(yōu)化設(shè)計(jì)的目的之一。 在覆蓋層壓應(yīng)力優(yōu)化上，優(yōu)化壩型一效果更好，且優(yōu)化壩型一、壩型二的極值點(diǎn)位置都在壩踵處，契合材料力學(xué)法所得結(jié)果，具有設(shè)計(jì)指導(dǎo)意義。 剪切破壞雖是土體主要破壞形式，但優(yōu)化壩型一抗剪狀態(tài)特征值為正值，可知實(shí)現(xiàn)了土體剪切破壞點(diǎn)數(shù)量的清零。 在其余特征值的對(duì)比上，優(yōu)化壩型一都比優(yōu)化壩型二更具優(yōu)勢(shì)。 綜上，優(yōu)化壩型一實(shí)現(xiàn)了斷面面積減少11.3%、覆蓋層最大壓應(yīng)力降低9.2%和抗剪破壞點(diǎn)清零，達(dá)到了材料用量減少和地基應(yīng)力狀態(tài)改良的目的，更符合設(shè)計(jì)需求。

4 結(jié)語(yǔ)

首先，本文從膠結(jié)砂礫石料和砂卵石土的材料本構(gòu)出發(fā)，介紹了鄧肯－張經(jīng)典與修正模型通過(guò)ANSYS實(shí)現(xiàn)的思路，并簡(jiǎn)要介紹了ANSYS 二次開(kāi)發(fā)的細(xì)節(jié)，基于此，對(duì)初設(shè)尺寸的壩型進(jìn)行了非線性有限元計(jì)算，獲得了竣工期與蓄水期的特征值，分析其應(yīng)力變形規(guī)律；其次，依據(jù)材料力學(xué)法和有限元分析結(jié)果，以覆蓋層抗壓、抗剪性能作為優(yōu)化設(shè)計(jì)目標(biāo)，運(yùn)用數(shù)據(jù)交互方式，實(shí)現(xiàn)ANSYS 與MATLAB 的混編，得到兩種優(yōu)化壩型；最后，將優(yōu)化壩型與初始?jí)涡偷母黜?xiàng)特征值進(jìn)行對(duì)比，評(píng)判優(yōu)化效果。

優(yōu)化目標(biāo)的設(shè)定應(yīng)符合工程設(shè)計(jì)需求，設(shè)計(jì)變量也不只局限于壩型尺寸，可擴(kuò)展到覆蓋層尺寸與結(jié)構(gòu)材料參數(shù)。 本文探究了ANSYS 有限元分析與MATLAB優(yōu)化計(jì)算結(jié)合的思路，所得優(yōu)化壩型與特征值可供相關(guān)設(shè)計(jì)單位借鑒，后期會(huì)針對(duì)不同的設(shè)計(jì)變量對(duì)各目標(biāo)函數(shù)的敏感度做進(jìn)一步的優(yōu)化設(shè)計(jì)。